1 /* glplux.h (LU-factorization, bignum arithmetic) */
2 
3 /***********************************************************************
4 *  This code is part of GLPK (GNU Linear Programming Kit).
5 *
6 *  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
7 *  2009, 2010 Andrew Makhorin, Department for Applied Informatics,
8 *  Moscow Aviation Institute, Moscow, Russia. All rights reserved.
9 *  E-mail: <mao@gnu.org>.
10 *
11 *  GLPK is free software: you can redistribute it and/or modify it
12 *  under the terms of the GNU General Public License as published by
13 *  the Free Software Foundation, either version 3 of the License, or
14 *  (at your option) any later version.
15 *
16 *  GLPK is distributed in the hope that it will be useful, but WITHOUT
17 *  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
18 *  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
19 *  License for more details.
20 *
21 *  You should have received a copy of the GNU General Public License
22 *  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
23 ***********************************************************************/
24 
25 #ifndef GLPLUX_H
26 #define GLPLUX_H
27 
28 #include "glpdmp.h"
29 #include "glpgmp.h"
30 
31 /*----------------------------------------------------------------------
32 // The structure LUX defines LU-factorization of a square matrix A,
33 // which is the following quartet:
34 //
35 //    [A] = (F, V, P, Q),                                            (1)
36 //
37 // where F and V are such matrices that
38 //
39 //    A = F * V,                                                     (2)
40 //
41 // and P and Q are such permutation matrices that the matrix
42 //
43 //    L = P * F * inv(P)                                             (3)
44 //
45 // is lower triangular with unity diagonal, and the matrix
46 //
47 //    U = P * V * Q                                                  (4)
48 //
49 // is upper triangular. All the matrices have the order n.
50 //
51 // The matrices F and V are stored in row/column-wise sparse format as
52 // row and column linked lists of non-zero elements. Unity elements on
53 // the main diagonal of the matrix F are not stored. Pivot elements of
54 // the matrix V (that correspond to diagonal elements of the matrix U)
55 // are also missing from the row and column lists and stored separately
56 // in an ordinary array.
57 //
58 // The permutation matrices P and Q are stored as ordinary arrays using
59 // both row- and column-like formats.
60 //
61 // The matrices L and U being completely defined by the matrices F, V,
62 // P, and Q are not stored explicitly.
63 //
64 // It is easy to show that the factorization (1)-(3) is some version of
65 // LU-factorization. Indeed, from (3) and (4) it follows that:
66 //
67 //    F = inv(P) * L * P,
68 //
69 //    V = inv(P) * U * inv(Q),
70 //
71 // and substitution into (2) gives:
72 //
73 //    A = F * V = inv(P) * L * U * inv(Q).
74 //
75 // For more details see the program documentation. */
76 
77 typedef struct LUX LUX;
78 typedef struct LUXELM LUXELM;
79 typedef struct LUXWKA LUXWKA;
80 
81 struct LUX
82 {     /* LU-factorization of a square matrix */
83       int n;
84       /* the order of matrices A, F, V, P, Q */
85       DMP *pool;
86       /* memory pool for elements of matrices F and V */
87       LUXELM **F_row; /* LUXELM *F_row[1+n]; */
88       /* F_row[0] is not used;
89          F_row[i], 1 <= i <= n, is a pointer to the list of elements in
90          i-th row of matrix F (diagonal elements are not stored) */
91       LUXELM **F_col; /* LUXELM *F_col[1+n]; */
92       /* F_col[0] is not used;
93          F_col[j], 1 <= j <= n, is a pointer to the list of elements in
94          j-th column of matrix F (diagonal elements are not stored) */
95       mpq_t *V_piv; /* mpq_t V_piv[1+n]; */
96       /* V_piv[0] is not used;
97          V_piv[p], 1 <= p <= n, is a pivot element v[p,q] corresponding
98          to a diagonal element u[k,k] of matrix U = P*V*Q (used on k-th
99          elimination step, k = 1, 2, ..., n) */
100       LUXELM **V_row; /* LUXELM *V_row[1+n]; */
101       /* V_row[0] is not used;
102          V_row[i], 1 <= i <= n, is a pointer to the list of elements in
103          i-th row of matrix V (except pivot elements) */
104       LUXELM **V_col; /* LUXELM *V_col[1+n]; */
105       /* V_col[0] is not used;
106          V_col[j], 1 <= j <= n, is a pointer to the list of elements in
107          j-th column of matrix V (except pivot elements) */
108       int *P_row; /* int P_row[1+n]; */
109       /* P_row[0] is not used;
110          P_row[i] = j means that p[i,j] = 1, where p[i,j] is an element
111          of permutation matrix P */
112       int *P_col; /* int P_col[1+n]; */
113       /* P_col[0] is not used;
114          P_col[j] = i means that p[i,j] = 1, where p[i,j] is an element
115          of permutation matrix P */
116       /* if i-th row or column of matrix F is i'-th row or column of
117          matrix L = P*F*inv(P), or if i-th row of matrix V is i'-th row
118          of matrix U = P*V*Q, then P_row[i'] = i and P_col[i] = i' */
119       int *Q_row; /* int Q_row[1+n]; */
120       /* Q_row[0] is not used;
121          Q_row[i] = j means that q[i,j] = 1, where q[i,j] is an element
122          of permutation matrix Q */
123       int *Q_col; /* int Q_col[1+n]; */
124       /* Q_col[0] is not used;
125          Q_col[j] = i means that q[i,j] = 1, where q[i,j] is an element
126          of permutation matrix Q */
127       /* if j-th column of matrix V is j'-th column of matrix U = P*V*Q,
128          then Q_row[j] = j' and Q_col[j'] = j */
129       int rank;
130       /* the (exact) rank of matrices A and V */
131 };
132 
133 struct LUXELM
134 {     /* element of matrix F or V */
135       int i;
136       /* row index, 1 <= i <= m */
137       int j;
138       /* column index, 1 <= j <= n */
139       mpq_t val;
140       /* numeric (non-zero) element value */
141       LUXELM *r_prev;
142       /* pointer to previous element in the same row */
143       LUXELM *r_next;
144       /* pointer to next element in the same row */
145       LUXELM *c_prev;
146       /* pointer to previous element in the same column */
147       LUXELM *c_next;
148       /* pointer to next element in the same column */
149 };
150 
151 struct LUXWKA
152 {     /* working area (used only during factorization) */
153       /* in order to efficiently implement Markowitz strategy and Duff
154          search technique there are two families {R[0], R[1], ..., R[n]}
155          and {C[0], C[1], ..., C[n]}; member R[k] is a set of active
156          rows of matrix V having k non-zeros, and member C[k] is a set
157          of active columns of matrix V having k non-zeros (in the active
158          submatrix); each set R[k] and C[k] is implemented as a separate
159          doubly linked list */
160       int *R_len; /* int R_len[1+n]; */
161       /* R_len[0] is not used;
162          R_len[i], 1 <= i <= n, is the number of non-zero elements in
163          i-th row of matrix V (that is the length of i-th row) */
164       int *R_head; /* int R_head[1+n]; */
165       /* R_head[k], 0 <= k <= n, is the number of a first row, which is
166          active and whose length is k */
167       int *R_prev; /* int R_prev[1+n]; */
168       /* R_prev[0] is not used;
169          R_prev[i], 1 <= i <= n, is the number of a previous row, which
170          is active and has the same length as i-th row */
171       int *R_next; /* int R_next[1+n]; */
172       /* R_prev[0] is not used;
173          R_prev[i], 1 <= i <= n, is the number of a next row, which is
174          active and has the same length as i-th row */
175       int *C_len; /* int C_len[1+n]; */
176       /* C_len[0] is not used;
177          C_len[j], 1 <= j <= n, is the number of non-zero elements in
178          j-th column of the active submatrix of matrix V (that is the
179          length of j-th column in the active submatrix) */
180       int *C_head; /* int C_head[1+n]; */
181       /* C_head[k], 0 <= k <= n, is the number of a first column, which
182          is active and whose length is k */
183       int *C_prev; /* int C_prev[1+n]; */
184       /* C_prev[0] is not used;
185          C_prev[j], 1 <= j <= n, is the number of a previous column,
186          which is active and has the same length as j-th column */
187       int *C_next; /* int C_next[1+n]; */
188       /* C_next[0] is not used;
189          C_next[j], 1 <= j <= n, is the number of a next column, which
190          is active and has the same length as j-th column */
191 };
192 
193 #define lux_create            _glp_lux_create
194 #define lux_decomp            _glp_lux_decomp
195 #define lux_f_solve           _glp_lux_f_solve
196 #define lux_v_solve           _glp_lux_v_solve
197 #define lux_solve             _glp_lux_solve
198 #define lux_delete            _glp_lux_delete
199 
200 LUX *lux_create(int n);
201 /* create LU-factorization */
202 
203 int lux_decomp(LUX *lux, int (*col)(void *info, int j, int ind[],
204       mpq_t val[]), void *info);
205 /* compute LU-factorization */
206 
207 void lux_f_solve(LUX *lux, int tr, mpq_t x[]);
208 /* solve system F*x = b or F'*x = b */
209 
210 void lux_v_solve(LUX *lux, int tr, mpq_t x[]);
211 /* solve system V*x = b or V'*x = b */
212 
213 void lux_solve(LUX *lux, int tr, mpq_t x[]);
214 /* solve system A*x = b or A'*x = b */
215 
216 void lux_delete(LUX *lux);
217 /* delete LU-factorization */
218 
219 #endif
220 
221 /* eof */
222