|
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 */ |