DyLP 1.10.4
Loading...
Searching...
No Matches
glpluf.h
Go to the documentation of this file.
1/* glpluf.h */
2
3/*----------------------------------------------------------------------
4-- Copyright (C) 2000, 2001, 2002 Andrew Makhorin <mao@mai2.rcnet.ru>,
5-- Department for Applied Informatics, Moscow Aviation
6-- Institute, Moscow, Russia. All rights reserved.
7--
8-- This file is a part of GLPK (GNU Linear Programming Kit).
9--
10-- Licensed under the Eclipse Public License (EPL) by permission of the
11-- author for inclusion in the DyLP LP distribution.
12----------------------------------------------------------------------*/
13
14/*
15 @(#)glpluf.h 1.1 10/18/02
16 svn/cvs: $Id: glpluf.h 407 2010-12-31 20:48:48Z lou $
17*/
18
19#ifndef _GLPLUF_H
20#define _GLPLUF_H
21
22#define luf_create dy_glp_luf_create
23#define luf_defrag_sva dy_glp_luf_defrag_sva
24#define luf_enlarge_row dy_glp_luf_enlarge_row
25#define luf_enlarge_col dy_glp_luf_enlarge_col
26#define luf_alloc_wa dy_glp_luf_alloc_wa
27#define luf_free_wa dy_glp_luf_free_wa
28#define luf_decomp dy_glp_luf_decomp
29#define luf_f_solve dy_glp_luf_f_solve
30#define luf_v_solve dy_glp_luf_v_solve
31#define luf_solve dy_glp_luf_solve
32#define luf_delete dy_glp_luf_delete
33
34/*----------------------------------------------------------------------
35-- The structure LUF defines LU-factorization of a square matrix A,
36-- which is the following quartet:
37--
38-- [A] = (F, V, P, Q), (1)
39--
40-- where F and V are such matrices that
41--
42-- A = F * V, (2)
43--
44-- and P and Q are such permutation matrices that the matrix
45--
46-- L = P * F * inv(P) (3)
47--
48-- is lower triangular with unity diagonal, and the matrix
49--
50-- U = P * V * Q (4)
51--
52-- is upper triangular. All the matrices have the order n.
53--
54-- The matrices F and V are stored in row/column-wise sparse format as
55-- row and column linked lists of non-zero elements. Unity elements on
56-- the main diagonal of the matrix F are not stored. Pivot elements of
57-- the matrix V (that correspond to diagonal elements of the matrix U)
58-- are also missing from the row and column lists and stored separately
59-- in an ordinary array.
60--
61-- The permutation matrices P and Q are stored as ordinary arrays using
62-- both row- and column-like formats.
63--
64-- The matrices L and U being completely defined by the matrices F, V,
65-- P, and Q are not stored explicitly.
66--
67-- It can easily be shown that the factorization (1)-(3) is a version of
68-- LU-factorization. Indeed, from (3) and (4) it follows that
69--
70-- F = inv(P) * L * P,
71--
72-- V = inv(P) * U * inv(Q),
73--
74-- and substitution into (2) gives
75--
76-- A = F * V = inv(P) * L * U * inv(Q).
77--
78-- For more details see the program documentation. */
79
80typedef struct LUF LUF;
81typedef struct LUF_WA LUF_WA;
82
83struct LUF
84{ /* LU-factorization of a square matrix */
85 int n;
86 /* order of the matrices A, F, V, P, Q */
87 int valid;
88 /* if this flag is not set, the factorization is invalid */
89 /*--------------------------------------------------------------*/
90 /* matrix F in row-wise format */
91 int *fr_ptr; /* int fr_ptr[1+n]; */
92 /* fr_ptr[0] is not used;
93 fr_ptr[i], i = 1, ..., n, is a pointer to the first element of
94 the i-th row in the sparse vector area */
95 int *fr_len; /* int fr_len[1+n]; */
96 /* fr_len[0] is not used;
97 fr_len[i], i = 1, ..., n, is number of elements in the i-th
98 row (except unity diagonal element) */
99 /*--------------------------------------------------------------*/
100 /* matrix F in column-wise format */
101 int *fc_ptr; /* int fc_ptr[1+n]; */
102 /* fc_ptr[0] is not used;
103 fc_ptr[j], j = 1, ..., n, is a pointer to the first element of
104 the j-th column in the sparse vector area */
105 int *fc_len; /* int fc_len[1+n]; */
106 /* fc_len[0] is not used;
107 fc_len[j], j = 1, ..., n, is number of elements in the j-th
108 column (except unity diagonal element) */
109 /*--------------------------------------------------------------*/
110 /* matrix V in row-wise format */
111 int *vr_ptr; /* int vr_ptr[1+n]; */
112 /* vr_ptr[0] is not used;
113 vr_ptr[i], i = 1, ..., n, is a pointer to the first element of
114 the i-th row in the sparse vector area */
115 int *vr_len; /* int vr_len[1+n]; */
116 /* vr_len[0] is not used;
117 vr_len[i], i = 1, ..., n, is number of elements in the i-th
118 row (except pivot element) */
119 int *vr_cap; /* int vr_cap[1+n]; */
120 /* vr_cap[0] is not used;
121 vr_cap[i], i = 1, ..., n, is capacity of the i-th row, i.e.
122 maximal number of elements, which can be stored there without
123 relocating the row, vr_cap[i] >= vr_len[i] */
124 double *vr_piv; /* double vr_piv[1+n]; */
125 /* vr_piv[0] is not used;
126 vr_piv[p], p = 1, ..., n, is the pivot element v[p,q], which
127 corresponds to a diagonal element of the matrix U = P*V*Q */
128 /*--------------------------------------------------------------*/
129 /* matrix V in column-wise format */
130 int *vc_ptr; /* int vc_ptr[1+n]; */
131 /* vc_ptr[0] is not used;
132 vc_ptr[j], j = 1, ..., n, is a pointer to the first element of
133 the j-th column in the sparse vector area */
134 int *vc_len; /* int vc_len[1+n]; */
135 /* vc_len[0] is not used;
136 vc_len[j], j = 1, ..., n, is number of elements in the j-th
137 column (except pivot element) */
138 int *vc_cap; /* int vc_cap[1+n]; */
139 /* vc_cap[0] is not used;
140 vc_cap[j], j = 1, ..., n, is capacity of the j-th column, i.e.
141 maximal number of elements, which can be stored there without
142 relocating the column, vc_cap[j] >= vc_len[j] */
143 /*--------------------------------------------------------------*/
144 /* matrix P */
145 int *pp_row; /* int pp_row[1+n]; */
146 /* pp_row[0] is not used; pp_row[i] = j means that p[i,j] = 1 */
147 int *pp_col; /* int pp_col[1+n]; */
148 /* pp_col[0] is not used; pp_col[j] = i means that p[i,j] = 1 */
149 /* if i-th row or column of the matrix F corresponds to i'-th row
150 or column of the matrix L = P*F*inv(P), or if i-th row of the
151 matrix V corresponds to i'-th row of the matrix U = P*V*Q, then
152 pp_row[i'] = i and pp_col[i] = i' */
153 /*--------------------------------------------------------------*/
154 /* matrix Q */
155 int *qq_row; /* int qq_row[1+n]; */
156 /* qq_row[0] is not used; qq_row[i] = j means that q[i,j] = 1 */
157 int *qq_col; /* int qq_col[1+n]; */
158 /* qq_col[0] is not used; qq_col[j] = i means that q[i,j] = 1 */
159 /* if j-th column of the matrix V corresponds to j'-th column of
160 the matrix U = P*V*Q, then qq_row[j] = j' and qq_col[j'] = j */
161 /*--------------------------------------------------------------*/
162 /* sparse vector area (SVA) is a set of locations intended to
163 store sparse vectors that represent rows and columns of the
164 matrices F and V; each location is the doublet (ndx, val),
165 where ndx is an index and val is a numerical value of a sparse
166 vector element; in the whole each sparse vector is a set of
167 adjacent locations defined by a pointer to the first element
168 and number of elements; these pointer and number are stored in
169 the corresponding matrix data structure (see above); the left
170 part of SVA is used to store rows and columns of the matrix V,
171 the right part is used to store rows and columns of the matrix
172 F; between the left and right parts there is the middle part,
173 locations of which are free */
175 /* total size of the sparse vector area, in locations; locations
176 are numbered by integers 1, 2, ..., sv_size, and location with
177 the number 0 is not used; if it is necessary, the size of SVA
178 is automatically increased */
180 /* SVA partitioning pointers:
181 locations 1, ..., sv_beg-1 belong to the left part;
182 locations sv_beg, ..., sv_end-1 belong to the middle part;
183 locations sv_end, ..., sv_size belong to the right part;
184 number of free locations, i.e. locations that belong to the
185 middle part, is (sv_end - sv_beg) */
186 int *sv_ndx; /* int sv_ndx[1+sv_size]; */
187 /* sv_ndx[0] is not used;
188 sv_ndx[k], 1 <= k <= sv_size, is the index field of the k-th
189 location */
190 double *sv_val; /* double sv_val[1+sv_size]; */
191 /* sv_val[0] is not used;
192 sv_val[k], 1 <= k <= sv_size, is the value field of the k-th
193 location */
194 /* in order to efficiently defragment the left part of SVA there
195 is a double linked list of rows and columns of the matrix V,
196 where rows have numbers 1, ..., n, and columns have numbers
197 n+1, ..., n+n, due to that each row and column can be uniquely
198 identified by one integer; in this list rows and columns are
199 ordered by ascending their pointers vr_ptr[i] and vc_ptr[j] */
201 /* the number of the leftmost row/column */
203 /* the number of the rightmost row/column */
204 int *sv_prev; /* int sv_prev[1+n+n]; */
205 /* sv_prev[k], k = 1, ..., n+n, is the number of a row/column,
206 which precedes the k-th row/column */
207 int *sv_next; /* int sv_next[1+n+n]; */
208 /* sv_next[k], k = 1, ..., n+n, is the number of a row/column,
209 which succedes the k-th row/column */
210 /*--------------------------------------------------------------*/
211 /* working arrays */
212 int *flag; /* int flag[1+n]; */
213 /* integer working array */
214 double *work; /* double work[1+n]; */
215 /* floating-point working array */
216 /*--------------------------------------------------------------*/
217 /* control parameters */
219 /* new required size of the sparse vector area, in locations; set
220 automatically by the factorizing routine */
221 double piv_tol;
222 /* threshold pivoting tolerance, 0 < piv_tol < 1; element v[i,j]
223 of the active submatrix fits to be pivot if it satisfies to the
224 stability condition |v[i,j]| >= piv_tol * max|v[i,*]|, i.e. if
225 this element is not very small (in absolute value) among other
226 elements in the same row; decreasing this parameter involves
227 better sparsity at the expense of numerical accuracy and vice
228 versa */
230 /* maximal allowable number of pivot candidates to be considered;
231 if piv_lim pivot candidates have been considered, the pivoting
232 routine terminates the search with the best candidate found */
233 int suhl;
234 /* if this flag is set, the pivoting routine applies a heuristic
235 rule proposed by Uwe Suhl: if a column of the active submatrix
236 has no eligible pivot candidates (i.e. all its elements don't
237 satisfy to the stability condition), the routine excludes such
238 column from the futher consideration until it becomes a column
239 singleton; in many cases this reduces a time needed for pivot
240 searching */
241 double eps_tol;
242 /* epsilon tolerance; each element of the matrix V with absolute
243 value less than eps_tol is replaced by exact zero */
244 double max_gro;
245 /* maximal allowable growth of elements of the matrix V during
246 all the factorization process; if on some elimination step the
247 ratio big_v / max_a (see below) becomes greater than max_gro,
248 the matrix A is considered as ill-conditioned (it is assumed
249 that the tolerance piv_tol has an adequate value) */
250 /*--------------------------------------------------------------*/
251 /* some statistics */
252 int nnz_a;
253 /* number of non-zeros in the matrix A */
254 int nnz_f;
255 /* number of non-zeros in the matrix F (except diagonal elements,
256 which are always equal to one and therefore not stored) */
257 int nnz_v;
258 /* number of non-zeros in the matrix V (except pivot elements,
259 which correspond to diagonal elements of the matrix U = P*V*Q
260 and which are stored separately in the array vr_piv) */
261 double max_a;
262 /* largest of absolute values of elements of the matrix A */
263 double big_v;
264 /* estimated largest of absolute values of elements appeared in
265 the active submatrix during all the factorization process */
266 int rank;
267 /* estimated rank of the matrix A */
268};
269
270struct LUF_WA
271{ /* working area (used only during factorization) */
272 double *rs_max; /* double rs_max[1+n]; */
273 /* rs_max[0] is not used; rs_max[i], 1 <= i <= n, is used only if
274 the i-th row of the matrix V belongs to the active submatrix
275 and is the largest of absolute values of elements in this row;
276 rs_max[i] < 0.0 means that the largest value is not known yet
277 and should be determined by the pivoting routine */
278 /*--------------------------------------------------------------*/
279 /* in order to efficiently implement Markowitz strategy and Duff
280 search technique there are two families {R[0], R[1], ..., R[n]}
281 and {C[0], C[1], ..., C[n]}; member R[k] is a set of active
282 rows of the matrix V, which have k non-zeros; similarly, member
283 C[k] is a set of active columns of the matrix V, which have k
284 non-zeros (in the active submatrix); each set R[k] and C[k] is
285 implemented as a separate doubly linked list */
286 int *rs_head; /* int rs_head[1+n]; */
287 /* rs_head[k], 0 <= k <= n, is number of the first active row,
288 which has k non-zeros */
289 int *rs_prev; /* int rs_prev[1+n]; */
290 /* rs_prev[0] is not used; rs_prev[i], 1 <= i <= n, is number of
291 the previous active row, which has the same number of non-zeros
292 as the i-th row */
293 int *rs_next; /* int rs_next[1+n]; */
294 /* rs_next[0] is not used; rs_next[i], 1 <= i <= n, is number of
295 the next active row, which has the same number of non-zeros as
296 the i-th row */
297 int *cs_head; /* int cs_head[1+n]; */
298 /* cs_head[k], 0 <= k <= n, is number of the first active column,
299 which has k non-zeros (in the active submatrix) */
300 int *cs_prev; /* int cs_prev[1+n]; */
301 /* cs_prev[0] is not used; cs_prev[j], 1 <= j <= n, is number of
302 the previous active column, which has the same number of
303 non-zeros (in the active submatrix) as the j-th column */
304 int *cs_next; /* int cs_next[1+n]; */
305 /* cs_next[0] is not used; cs_next[j], 1 <= j <= n, is number of
306 the next active column, which has the same number of non-zeros
307 (in the active submatrix) as the j-th column */
308};
309
310LUF *luf_create(int n, int sv_size);
311/* create LU-factorization */
312
314/* defragment the sparse vector area */
315
316int luf_enlarge_row(LUF *luf, int i, int cap);
317/* enlarge row capacity */
318
319int luf_enlarge_col(LUF *luf, int j, int cap);
320/* enlarge column capacity */
321
323/* pre-allocate working area */
324
326/* free working area */
327
329 void *info, int (*col)(void *info, int j, int rn[], double aj[]),
330 LUF_WA *wa);
331/* compute LU-factorization */
332
333void luf_f_solve(LUF *luf, int tr, double x[]);
334/* solve system F*x = b or F'*x = b */
335
336void luf_v_solve(LUF *luf, int tr, double x[]);
337/* solve system V*x = b or V'*x = b */
338
339void luf_solve(LUF *luf, int tr, double x[]);
340/* solve system A*x = b or A'*x = b */
341
342void luf_delete(LUF *luf);
343/* delete LU-factorization */
344
345#endif
346
347/* eof */
#define luf_f_solve
Definition glpluf.h:29
#define luf_defrag_sva
Definition glpluf.h:23
#define luf_v_solve
Definition glpluf.h:30
#define luf_create
Definition glpluf.h:22
#define luf_delete
Definition glpluf.h:32
#define luf_solve
Definition glpluf.h:31
#define luf_enlarge_row
Definition glpluf.h:24
#define luf_alloc_wa
Definition glpluf.h:26
#define luf_decomp
Definition glpluf.h:28
#define luf_free_wa
Definition glpluf.h:27
#define luf_enlarge_col
Definition glpluf.h:25
int * cs_head
Definition glpluf.h:297
int * cs_next
Definition glpluf.h:304
int * rs_head
Definition glpluf.h:286
int * rs_prev
Definition glpluf.h:289
int * cs_prev
Definition glpluf.h:300
int * rs_next
Definition glpluf.h:293
double * rs_max
Definition glpluf.h:272
Definition glpluf.h:84
int sv_head
Definition glpluf.h:200
double big_v
Definition glpluf.h:263
int * sv_prev
Definition glpluf.h:204
int rank
Definition glpluf.h:266
int * vc_ptr
Definition glpluf.h:130
double piv_tol
Definition glpluf.h:221
int * fr_ptr
Definition glpluf.h:91
int * qq_row
Definition glpluf.h:155
int new_sva
Definition glpluf.h:218
int * pp_col
Definition glpluf.h:147
int * fc_len
Definition glpluf.h:105
double eps_tol
Definition glpluf.h:241
int sv_beg
Definition glpluf.h:179
int nnz_f
Definition glpluf.h:254
int suhl
Definition glpluf.h:233
double * work
Definition glpluf.h:214
int * vr_cap
Definition glpluf.h:119
double max_gro
Definition glpluf.h:244
int * vc_len
Definition glpluf.h:134
int piv_lim
Definition glpluf.h:229
int n
Definition glpluf.h:85
double * sv_val
Definition glpluf.h:190
int * pp_row
Definition glpluf.h:145
int nnz_v
Definition glpluf.h:257
int * vr_len
Definition glpluf.h:115
int * vc_cap
Definition glpluf.h:138
int * fr_len
Definition glpluf.h:95
int * qq_col
Definition glpluf.h:157
int nnz_a
Definition glpluf.h:252
int * flag
Definition glpluf.h:212
int * sv_next
Definition glpluf.h:207
double * vr_piv
Definition glpluf.h:124
int sv_size
Definition glpluf.h:174
int sv_end
Definition glpluf.h:179
int sv_tail
Definition glpluf.h:202
int valid
Definition glpluf.h:87
double max_a
Definition glpluf.h:261
int * vr_ptr
Definition glpluf.h:111
int * fc_ptr
Definition glpluf.h:101
int * sv_ndx
Definition glpluf.h:186