|
1 /* glpluf.h (LU-factorization) */ |
|
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 GLPLUF_H |
|
26 #define GLPLUF_H |
|
27 |
|
28 /*********************************************************************** |
|
29 * The structure LUF defines LU-factorization of a square matrix A and |
|
30 * is the following quartet: |
|
31 * |
|
32 * [A] = (F, V, P, Q), (1) |
|
33 * |
|
34 * where F and V are such matrices that |
|
35 * |
|
36 * A = F * V, (2) |
|
37 * |
|
38 * and P and Q are such permutation matrices that the matrix |
|
39 * |
|
40 * L = P * F * inv(P) (3) |
|
41 * |
|
42 * is lower triangular with unity diagonal, and the matrix |
|
43 * |
|
44 * U = P * V * Q (4) |
|
45 * |
|
46 * is upper triangular. All the matrices have the order n. |
|
47 * |
|
48 * Matrices F and V are stored in row- and column-wise sparse format |
|
49 * as row and column linked lists of non-zero elements. Unity elements |
|
50 * on the main diagonal of matrix F are not stored. Pivot elements of |
|
51 * matrix V (which correspond to diagonal elements of matrix U) are |
|
52 * stored separately in an ordinary array. |
|
53 * |
|
54 * Permutation matrices P and Q are stored in ordinary arrays in both |
|
55 * row- and column-like formats. |
|
56 * |
|
57 * Matrices L and U are completely defined by matrices F, V, P, and Q |
|
58 * and therefore not stored explicitly. |
|
59 * |
|
60 * The factorization (1)-(4) is a version of LU-factorization. Indeed, |
|
61 * from (3) and (4) it follows that: |
|
62 * |
|
63 * F = inv(P) * L * P, |
|
64 * |
|
65 * U = inv(P) * U * inv(Q), |
|
66 * |
|
67 * and substitution into (2) leads to: |
|
68 * |
|
69 * A = F * V = inv(P) * L * U * inv(Q). |
|
70 * |
|
71 * For more details see the program documentation. */ |
|
72 |
|
73 typedef struct LUF LUF; |
|
74 |
|
75 struct LUF |
|
76 { /* LU-factorization of a square matrix */ |
|
77 int n_max; |
|
78 /* maximal value of n (increased automatically, if necessary) */ |
|
79 int n; |
|
80 /* the order of matrices A, F, V, P, Q */ |
|
81 int valid; |
|
82 /* the factorization is valid only if this flag is set */ |
|
83 /*--------------------------------------------------------------*/ |
|
84 /* matrix F in row-wise format */ |
|
85 int *fr_ptr; /* int fr_ptr[1+n_max]; */ |
|
86 /* fr_ptr[i], i = 1,...,n, is a pointer to the first element of |
|
87 i-th row in SVA */ |
|
88 int *fr_len; /* int fr_len[1+n_max]; */ |
|
89 /* fr_len[i], i = 1,...,n, is the number of elements in i-th row |
|
90 (except unity diagonal element) */ |
|
91 /*--------------------------------------------------------------*/ |
|
92 /* matrix F in column-wise format */ |
|
93 int *fc_ptr; /* int fc_ptr[1+n_max]; */ |
|
94 /* fc_ptr[j], j = 1,...,n, is a pointer to the first element of |
|
95 j-th column in SVA */ |
|
96 int *fc_len; /* int fc_len[1+n_max]; */ |
|
97 /* fc_len[j], j = 1,...,n, is the number of elements in j-th |
|
98 column (except unity diagonal element) */ |
|
99 /*--------------------------------------------------------------*/ |
|
100 /* matrix V in row-wise format */ |
|
101 int *vr_ptr; /* int vr_ptr[1+n_max]; */ |
|
102 /* vr_ptr[i], i = 1,...,n, is a pointer to the first element of |
|
103 i-th row in SVA */ |
|
104 int *vr_len; /* int vr_len[1+n_max]; */ |
|
105 /* vr_len[i], i = 1,...,n, is the number of elements in i-th row |
|
106 (except pivot element) */ |
|
107 int *vr_cap; /* int vr_cap[1+n_max]; */ |
|
108 /* vr_cap[i], i = 1,...,n, is the capacity of i-th row, i.e. |
|
109 maximal number of elements which can be stored in the row |
|
110 without relocating it, vr_cap[i] >= vr_len[i] */ |
|
111 double *vr_piv; /* double vr_piv[1+n_max]; */ |
|
112 /* vr_piv[p], p = 1,...,n, is the pivot element v[p,q] which |
|
113 corresponds to a diagonal element of matrix U = P*V*Q */ |
|
114 /*--------------------------------------------------------------*/ |
|
115 /* matrix V in column-wise format */ |
|
116 int *vc_ptr; /* int vc_ptr[1+n_max]; */ |
|
117 /* vc_ptr[j], j = 1,...,n, is a pointer to the first element of |
|
118 j-th column in SVA */ |
|
119 int *vc_len; /* int vc_len[1+n_max]; */ |
|
120 /* vc_len[j], j = 1,...,n, is the number of elements in j-th |
|
121 column (except pivot element) */ |
|
122 int *vc_cap; /* int vc_cap[1+n_max]; */ |
|
123 /* vc_cap[j], j = 1,...,n, is the capacity of j-th column, i.e. |
|
124 maximal number of elements which can be stored in the column |
|
125 without relocating it, vc_cap[j] >= vc_len[j] */ |
|
126 /*--------------------------------------------------------------*/ |
|
127 /* matrix P */ |
|
128 int *pp_row; /* int pp_row[1+n_max]; */ |
|
129 /* pp_row[i] = j means that P[i,j] = 1 */ |
|
130 int *pp_col; /* int pp_col[1+n_max]; */ |
|
131 /* pp_col[j] = i means that P[i,j] = 1 */ |
|
132 /* if i-th row or column of matrix F is i'-th row or column of |
|
133 matrix L, or if i-th row of matrix V is i'-th row of matrix U, |
|
134 then pp_row[i'] = i and pp_col[i] = i' */ |
|
135 /*--------------------------------------------------------------*/ |
|
136 /* matrix Q */ |
|
137 int *qq_row; /* int qq_row[1+n_max]; */ |
|
138 /* qq_row[i] = j means that Q[i,j] = 1 */ |
|
139 int *qq_col; /* int qq_col[1+n_max]; */ |
|
140 /* qq_col[j] = i means that Q[i,j] = 1 */ |
|
141 /* if j-th column of matrix V is j'-th column of matrix U, then |
|
142 qq_row[j] = j' and qq_col[j'] = j */ |
|
143 /*--------------------------------------------------------------*/ |
|
144 /* the Sparse Vector Area (SVA) is a set of locations used to |
|
145 store sparse vectors representing rows and columns of matrices |
|
146 F and V; each location is a doublet (ind, val), where ind is |
|
147 an index, and val is a numerical value of a sparse vector |
|
148 element; in the whole each sparse vector is a set of adjacent |
|
149 locations defined by a pointer to the first element and the |
|
150 number of elements; these pointer and number are stored in the |
|
151 corresponding matrix data structure (see above); the left part |
|
152 of SVA is used to store rows and columns of matrix V, and its |
|
153 right part is used to store rows and columns of matrix F; the |
|
154 middle part of SVA contains free (unused) locations */ |
|
155 int sv_size; |
|
156 /* the size of SVA, in locations; all locations are numbered by |
|
157 integers 1, ..., n, and location 0 is not used; if necessary, |
|
158 the SVA size is automatically increased */ |
|
159 int sv_beg, sv_end; |
|
160 /* SVA partitioning pointers: |
|
161 locations from 1 to sv_beg-1 belong to the left part |
|
162 locations from sv_beg to sv_end-1 belong to the middle part |
|
163 locations from sv_end to sv_size belong to the right part |
|
164 the size of the middle part is (sv_end - sv_beg) */ |
|
165 int *sv_ind; /* sv_ind[1+sv_size]; */ |
|
166 /* sv_ind[k], 1 <= k <= sv_size, is the index field of k-th |
|
167 location */ |
|
168 double *sv_val; /* sv_val[1+sv_size]; */ |
|
169 /* sv_val[k], 1 <= k <= sv_size, is the value field of k-th |
|
170 location */ |
|
171 /*--------------------------------------------------------------*/ |
|
172 /* in order to efficiently defragment the left part of SVA there |
|
173 is a doubly linked list of rows and columns of matrix V, where |
|
174 rows are numbered by 1, ..., n, while columns are numbered by |
|
175 n+1, ..., n+n, that allows uniquely identifying each row and |
|
176 column of V by only one integer; in this list rows and columns |
|
177 are ordered by ascending their pointers vr_ptr and vc_ptr */ |
|
178 int sv_head; |
|
179 /* the number of leftmost row/column */ |
|
180 int sv_tail; |
|
181 /* the number of rightmost row/column */ |
|
182 int *sv_prev; /* int sv_prev[1+n_max+n_max]; */ |
|
183 /* sv_prev[k], k = 1,...,n+n, is the number of a row/column which |
|
184 precedes k-th row/column */ |
|
185 int *sv_next; /* int sv_next[1+n_max+n_max]; */ |
|
186 /* sv_next[k], k = 1,...,n+n, is the number of a row/column which |
|
187 succedes k-th row/column */ |
|
188 /*--------------------------------------------------------------*/ |
|
189 /* working segment (used only during factorization) */ |
|
190 double *vr_max; /* int vr_max[1+n_max]; */ |
|
191 /* vr_max[i], 1 <= i <= n, is used only if i-th row of matrix V |
|
192 is active (i.e. belongs to the active submatrix), and is the |
|
193 largest magnitude of elements in i-th row; if vr_max[i] < 0, |
|
194 the largest magnitude is not known yet and should be computed |
|
195 by the pivoting routine */ |
|
196 /*--------------------------------------------------------------*/ |
|
197 /* in order to efficiently implement Markowitz strategy and Duff |
|
198 search technique there are two families {R[0], R[1], ..., R[n]} |
|
199 and {C[0], C[1], ..., C[n]}; member R[k] is the set of active |
|
200 rows of matrix V, which have k non-zeros, and member C[k] is |
|
201 the set of active columns of V, which have k non-zeros in the |
|
202 active submatrix (i.e. in the active rows); each set R[k] and |
|
203 C[k] is implemented as a separate doubly linked list */ |
|
204 int *rs_head; /* int rs_head[1+n_max]; */ |
|
205 /* rs_head[k], 0 <= k <= n, is the number of first active row, |
|
206 which has k non-zeros */ |
|
207 int *rs_prev; /* int rs_prev[1+n_max]; */ |
|
208 /* rs_prev[i], 1 <= i <= n, is the number of previous row, which |
|
209 has the same number of non-zeros as i-th row */ |
|
210 int *rs_next; /* int rs_next[1+n_max]; */ |
|
211 /* rs_next[i], 1 <= i <= n, is the number of next row, which has |
|
212 the same number of non-zeros as i-th row */ |
|
213 int *cs_head; /* int cs_head[1+n_max]; */ |
|
214 /* cs_head[k], 0 <= k <= n, is the number of first active column, |
|
215 which has k non-zeros (in the active rows) */ |
|
216 int *cs_prev; /* int cs_prev[1+n_max]; */ |
|
217 /* cs_prev[j], 1 <= j <= n, is the number of previous column, |
|
218 which has the same number of non-zeros (in the active rows) as |
|
219 j-th column */ |
|
220 int *cs_next; /* int cs_next[1+n_max]; */ |
|
221 /* cs_next[j], 1 <= j <= n, is the number of next column, which |
|
222 has the same number of non-zeros (in the active rows) as j-th |
|
223 column */ |
|
224 /* (end of working segment) */ |
|
225 /*--------------------------------------------------------------*/ |
|
226 /* working arrays */ |
|
227 int *flag; /* int flag[1+n_max]; */ |
|
228 /* integer working array */ |
|
229 double *work; /* double work[1+n_max]; */ |
|
230 /* floating-point working array */ |
|
231 /*--------------------------------------------------------------*/ |
|
232 /* control parameters */ |
|
233 int new_sva; |
|
234 /* new required size of the sparse vector area, in locations; set |
|
235 automatically by the factorizing routine */ |
|
236 double piv_tol; |
|
237 /* threshold pivoting tolerance, 0 < piv_tol < 1; element v[i,j] |
|
238 of the active submatrix fits to be pivot if it satisfies to the |
|
239 stability criterion |v[i,j]| >= piv_tol * max |v[i,*]|, i.e. if |
|
240 it is not very small in the magnitude among other elements in |
|
241 the same row; decreasing this parameter gives better sparsity |
|
242 at the expense of numerical accuracy and vice versa */ |
|
243 int piv_lim; |
|
244 /* maximal allowable number of pivot candidates to be considered; |
|
245 if piv_lim pivot candidates have been considered, the pivoting |
|
246 routine terminates the search with the best candidate found */ |
|
247 int suhl; |
|
248 /* if this flag is set, the pivoting routine applies a heuristic |
|
249 proposed by Uwe Suhl: if a column of the active submatrix has |
|
250 no eligible pivot candidates (i.e. all its elements do not |
|
251 satisfy to the stability criterion), the routine excludes it |
|
252 from futher consideration until it becomes column singleton; |
|
253 in many cases this allows reducing the time needed for pivot |
|
254 searching */ |
|
255 double eps_tol; |
|
256 /* epsilon tolerance; each element of the active submatrix, whose |
|
257 magnitude is less than eps_tol, is replaced by exact zero */ |
|
258 double max_gro; |
|
259 /* maximal allowable growth of elements of matrix V during all |
|
260 the factorization process; if on some eliminaion step the ratio |
|
261 big_v / max_a (see below) becomes greater than max_gro, matrix |
|
262 A is considered as ill-conditioned (assuming that the pivoting |
|
263 tolerance piv_tol has an appropriate value) */ |
|
264 /*--------------------------------------------------------------*/ |
|
265 /* some statistics */ |
|
266 int nnz_a; |
|
267 /* the number of non-zeros in matrix A */ |
|
268 int nnz_f; |
|
269 /* the number of non-zeros in matrix F (except diagonal elements, |
|
270 which are not stored) */ |
|
271 int nnz_v; |
|
272 /* the number of non-zeros in matrix V (except its pivot elements, |
|
273 which are stored in a separate array) */ |
|
274 double max_a; |
|
275 /* the largest magnitude of elements of matrix A */ |
|
276 double big_v; |
|
277 /* the largest magnitude of elements of matrix V appeared in the |
|
278 active submatrix during all the factorization process */ |
|
279 int rank; |
|
280 /* estimated rank of matrix A */ |
|
281 }; |
|
282 |
|
283 /* return codes: */ |
|
284 #define LUF_ESING 1 /* singular matrix */ |
|
285 #define LUF_ECOND 2 /* ill-conditioned matrix */ |
|
286 |
|
287 #define luf_create_it _glp_luf_create_it |
|
288 LUF *luf_create_it(void); |
|
289 /* create LU-factorization */ |
|
290 |
|
291 #define luf_defrag_sva _glp_luf_defrag_sva |
|
292 void luf_defrag_sva(LUF *luf); |
|
293 /* defragment the sparse vector area */ |
|
294 |
|
295 #define luf_enlarge_row _glp_luf_enlarge_row |
|
296 int luf_enlarge_row(LUF *luf, int i, int cap); |
|
297 /* enlarge row capacity */ |
|
298 |
|
299 #define luf_enlarge_col _glp_luf_enlarge_col |
|
300 int luf_enlarge_col(LUF *luf, int j, int cap); |
|
301 /* enlarge column capacity */ |
|
302 |
|
303 #define luf_factorize _glp_luf_factorize |
|
304 int luf_factorize(LUF *luf, int n, int (*col)(void *info, int j, |
|
305 int ind[], double val[]), void *info); |
|
306 /* compute LU-factorization */ |
|
307 |
|
308 #define luf_f_solve _glp_luf_f_solve |
|
309 void luf_f_solve(LUF *luf, int tr, double x[]); |
|
310 /* solve system F*x = b or F'*x = b */ |
|
311 |
|
312 #define luf_v_solve _glp_luf_v_solve |
|
313 void luf_v_solve(LUF *luf, int tr, double x[]); |
|
314 /* solve system V*x = b or V'*x = b */ |
|
315 |
|
316 #define luf_a_solve _glp_luf_a_solve |
|
317 void luf_a_solve(LUF *luf, int tr, double x[]); |
|
318 /* solve system A*x = b or A'*x = b */ |
|
319 |
|
320 #define luf_delete_it _glp_luf_delete_it |
|
321 void luf_delete_it(LUF *luf); |
|
322 /* delete LU-factorization */ |
|
323 |
|
324 #endif |
|
325 |
|
326 /* eof */ |