lemon-project-template-glpk

comparison deps/glpk/src/glpluf.h @ 9:33de93886c88

find changesets by author, revision, files, or words in the commit message
Import GLPK 4.47
author Alpar Juttner <alpar@cs.elte.hu>
date Sun, 06 Nov 2011 20:59:10 +0100
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:af8ec7ccad64
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, 2011 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 */