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