00001 /* glpluf.h */ 00002 00003 /*---------------------------------------------------------------------- 00004 -- Copyright (C) 2000, 2001, 2002 Andrew Makhorin <mao@mai2.rcnet.ru>, 00005 -- Department for Applied Informatics, Moscow Aviation 00006 -- Institute, Moscow, Russia. All rights reserved. 00007 -- 00008 -- This file is a part of GLPK (GNU Linear Programming Kit). 00009 -- 00010 -- Licensed under the Common Public License (CPL) by permission of the 00011 -- author for inclusion in the DyLP LP distribution. 00012 ----------------------------------------------------------------------*/ 00013 00014 /* 00015 @(#)glpluf.h 1.1 10/18/02 00016 svn/cvs: $Id: glpluf.h 158 2007-07-06 01:25:14Z lou $ 00017 */ 00018 00019 #ifndef _GLPLUF_H 00020 #define _GLPLUF_H 00021 00022 #define luf_create dy_glp_luf_create 00023 #define luf_defrag_sva dy_glp_luf_defrag_sva 00024 #define luf_enlarge_row dy_glp_luf_enlarge_row 00025 #define luf_enlarge_col dy_glp_luf_enlarge_col 00026 #define luf_alloc_wa dy_glp_luf_alloc_wa 00027 #define luf_free_wa dy_glp_luf_free_wa 00028 #define luf_decomp dy_glp_luf_decomp 00029 #define luf_f_solve dy_glp_luf_f_solve 00030 #define luf_v_solve dy_glp_luf_v_solve 00031 #define luf_solve dy_glp_luf_solve 00032 #define luf_delete dy_glp_luf_delete 00033 00034 /*---------------------------------------------------------------------- 00035 -- The structure LUF defines LU-factorization of a square matrix A, 00036 -- which is the following quartet: 00037 -- 00038 -- [A] = (F, V, P, Q), (1) 00039 -- 00040 -- where F and V are such matrices that 00041 -- 00042 -- A = F * V, (2) 00043 -- 00044 -- and P and Q are such permutation matrices that the matrix 00045 -- 00046 -- L = P * F * inv(P) (3) 00047 -- 00048 -- is lower triangular with unity diagonal, and the matrix 00049 -- 00050 -- U = P * V * Q (4) 00051 -- 00052 -- is upper triangular. All the matrices have the order n. 00053 -- 00054 -- The matrices F and V are stored in row/column-wise sparse format as 00055 -- row and column linked lists of non-zero elements. Unity elements on 00056 -- the main diagonal of the matrix F are not stored. Pivot elements of 00057 -- the matrix V (that correspond to diagonal elements of the matrix U) 00058 -- are also missing from the row and column lists and stored separately 00059 -- in an ordinary array. 00060 -- 00061 -- The permutation matrices P and Q are stored as ordinary arrays using 00062 -- both row- and column-like formats. 00063 -- 00064 -- The matrices L and U being completely defined by the matrices F, V, 00065 -- P, and Q are not stored explicitly. 00066 -- 00067 -- It can easily be shown that the factorization (1)-(3) is a version of 00068 -- LU-factorization. Indeed, from (3) and (4) it follows that 00069 -- 00070 -- F = inv(P) * L * P, 00071 -- 00072 -- V = inv(P) * U * inv(Q), 00073 -- 00074 -- and substitution into (2) gives 00075 -- 00076 -- A = F * V = inv(P) * L * U * inv(Q). 00077 -- 00078 -- For more details see the program documentation. */ 00079 00080 typedef struct LUF LUF; 00081 typedef struct LUF_WA LUF_WA; 00082 00083 struct LUF 00084 { /* LU-factorization of a square matrix */ 00085 int n; 00086 /* order of the matrices A, F, V, P, Q */ 00087 int valid; 00088 /* if this flag is not set, the factorization is invalid */ 00089 /*--------------------------------------------------------------*/ 00090 /* matrix F in row-wise format */ 00091 int *fr_ptr; /* int fr_ptr[1+n]; */ 00092 /* fr_ptr[0] is not used; 00093 fr_ptr[i], i = 1, ..., n, is a pointer to the first element of 00094 the i-th row in the sparse vector area */ 00095 int *fr_len; /* int fr_len[1+n]; */ 00096 /* fr_len[0] is not used; 00097 fr_len[i], i = 1, ..., n, is number of elements in the i-th 00098 row (except unity diagonal element) */ 00099 /*--------------------------------------------------------------*/ 00100 /* matrix F in column-wise format */ 00101 int *fc_ptr; /* int fc_ptr[1+n]; */ 00102 /* fc_ptr[0] is not used; 00103 fc_ptr[j], j = 1, ..., n, is a pointer to the first element of 00104 the j-th column in the sparse vector area */ 00105 int *fc_len; /* int fc_len[1+n]; */ 00106 /* fc_len[0] is not used; 00107 fc_len[j], j = 1, ..., n, is number of elements in the j-th 00108 column (except unity diagonal element) */ 00109 /*--------------------------------------------------------------*/ 00110 /* matrix V in row-wise format */ 00111 int *vr_ptr; /* int vr_ptr[1+n]; */ 00112 /* vr_ptr[0] is not used; 00113 vr_ptr[i], i = 1, ..., n, is a pointer to the first element of 00114 the i-th row in the sparse vector area */ 00115 int *vr_len; /* int vr_len[1+n]; */ 00116 /* vr_len[0] is not used; 00117 vr_len[i], i = 1, ..., n, is number of elements in the i-th 00118 row (except pivot element) */ 00119 int *vr_cap; /* int vr_cap[1+n]; */ 00120 /* vr_cap[0] is not used; 00121 vr_cap[i], i = 1, ..., n, is capacity of the i-th row, i.e. 00122 maximal number of elements, which can be stored there without 00123 relocating the row, vr_cap[i] >= vr_len[i] */ 00124 double *vr_piv; /* double vr_piv[1+n]; */ 00125 /* vr_piv[0] is not used; 00126 vr_piv[p], p = 1, ..., n, is the pivot element v[p,q], which 00127 corresponds to a diagonal element of the matrix U = P*V*Q */ 00128 /*--------------------------------------------------------------*/ 00129 /* matrix V in column-wise format */ 00130 int *vc_ptr; /* int vc_ptr[1+n]; */ 00131 /* vc_ptr[0] is not used; 00132 vc_ptr[j], j = 1, ..., n, is a pointer to the first element of 00133 the j-th column in the sparse vector area */ 00134 int *vc_len; /* int vc_len[1+n]; */ 00135 /* vc_len[0] is not used; 00136 vc_len[j], j = 1, ..., n, is number of elements in the j-th 00137 column (except pivot element) */ 00138 int *vc_cap; /* int vc_cap[1+n]; */ 00139 /* vc_cap[0] is not used; 00140 vc_cap[j], j = 1, ..., n, is capacity of the j-th column, i.e. 00141 maximal number of elements, which can be stored there without 00142 relocating the column, vc_cap[j] >= vc_len[j] */ 00143 /*--------------------------------------------------------------*/ 00144 /* matrix P */ 00145 int *pp_row; /* int pp_row[1+n]; */ 00146 /* pp_row[0] is not used; pp_row[i] = j means that p[i,j] = 1 */ 00147 int *pp_col; /* int pp_col[1+n]; */ 00148 /* pp_col[0] is not used; pp_col[j] = i means that p[i,j] = 1 */ 00149 /* if i-th row or column of the matrix F corresponds to i'-th row 00150 or column of the matrix L = P*F*inv(P), or if i-th row of the 00151 matrix V corresponds to i'-th row of the matrix U = P*V*Q, then 00152 pp_row[i'] = i and pp_col[i] = i' */ 00153 /*--------------------------------------------------------------*/ 00154 /* matrix Q */ 00155 int *qq_row; /* int qq_row[1+n]; */ 00156 /* qq_row[0] is not used; qq_row[i] = j means that q[i,j] = 1 */ 00157 int *qq_col; /* int qq_col[1+n]; */ 00158 /* qq_col[0] is not used; qq_col[j] = i means that q[i,j] = 1 */ 00159 /* if j-th column of the matrix V corresponds to j'-th column of 00160 the matrix U = P*V*Q, then qq_row[j] = j' and qq_col[j'] = j */ 00161 /*--------------------------------------------------------------*/ 00162 /* sparse vector area (SVA) is a set of locations intended to 00163 store sparse vectors that represent rows and columns of the 00164 matrices F and V; each location is the doublet (ndx, val), 00165 where ndx is an index and val is a numerical value of a sparse 00166 vector element; in the whole each sparse vector is a set of 00167 adjacent locations defined by a pointer to the first element 00168 and number of elements; these pointer and number are stored in 00169 the corresponding matrix data structure (see above); the left 00170 part of SVA is used to store rows and columns of the matrix V, 00171 the right part is used to store rows and columns of the matrix 00172 F; between the left and right parts there is the middle part, 00173 locations of which are free */ 00174 int sv_size; 00175 /* total size of the sparse vector area, in locations; locations 00176 are numbered by integers 1, 2, ..., sv_size, and location with 00177 the number 0 is not used; if it is necessary, the size of SVA 00178 is automatically increased */ 00179 int sv_beg, sv_end; 00180 /* SVA partitioning pointers: 00181 locations 1, ..., sv_beg-1 belong to the left part; 00182 locations sv_beg, ..., sv_end-1 belong to the middle part; 00183 locations sv_end, ..., sv_size belong to the right part; 00184 number of free locations, i.e. locations that belong to the 00185 middle part, is (sv_end - sv_beg) */ 00186 int *sv_ndx; /* int sv_ndx[1+sv_size]; */ 00187 /* sv_ndx[0] is not used; 00188 sv_ndx[k], 1 <= k <= sv_size, is the index field of the k-th 00189 location */ 00190 double *sv_val; /* double sv_val[1+sv_size]; */ 00191 /* sv_val[0] is not used; 00192 sv_val[k], 1 <= k <= sv_size, is the value field of the k-th 00193 location */ 00194 /* in order to efficiently defragment the left part of SVA there 00195 is a double linked list of rows and columns of the matrix V, 00196 where rows have numbers 1, ..., n, and columns have numbers 00197 n+1, ..., n+n, due to that each row and column can be uniquely 00198 identified by one integer; in this list rows and columns are 00199 ordered by ascending their pointers vr_ptr[i] and vc_ptr[j] */ 00200 int sv_head; 00201 /* the number of the leftmost row/column */ 00202 int sv_tail; 00203 /* the number of the rightmost row/column */ 00204 int *sv_prev; /* int sv_prev[1+n+n]; */ 00205 /* sv_prev[k], k = 1, ..., n+n, is the number of a row/column, 00206 which precedes the k-th row/column */ 00207 int *sv_next; /* int sv_next[1+n+n]; */ 00208 /* sv_next[k], k = 1, ..., n+n, is the number of a row/column, 00209 which succedes the k-th row/column */ 00210 /*--------------------------------------------------------------*/ 00211 /* working arrays */ 00212 int *flag; /* int flag[1+n]; */ 00213 /* integer working array */ 00214 double *work; /* double work[1+n]; */ 00215 /* floating-point working array */ 00216 /*--------------------------------------------------------------*/ 00217 /* control parameters */ 00218 int new_sva; 00219 /* new required size of the sparse vector area, in locations; set 00220 automatically by the factorizing routine */ 00221 double piv_tol; 00222 /* threshold pivoting tolerance, 0 < piv_tol < 1; element v[i,j] 00223 of the active submatrix fits to be pivot if it satisfies to the 00224 stability condition |v[i,j]| >= piv_tol * max|v[i,*]|, i.e. if 00225 this element is not very small (in absolute value) among other 00226 elements in the same row; decreasing this parameter involves 00227 better sparsity at the expense of numerical accuracy and vice 00228 versa */ 00229 int piv_lim; 00230 /* maximal allowable number of pivot candidates to be considered; 00231 if piv_lim pivot candidates have been considered, the pivoting 00232 routine terminates the search with the best candidate found */ 00233 int suhl; 00234 /* if this flag is set, the pivoting routine applies a heuristic 00235 rule proposed by Uwe Suhl: if a column of the active submatrix 00236 has no eligible pivot candidates (i.e. all its elements don't 00237 satisfy to the stability condition), the routine excludes such 00238 column from the futher consideration until it becomes a column 00239 singleton; in many cases this reduces a time needed for pivot 00240 searching */ 00241 double eps_tol; 00242 /* epsilon tolerance; each element of the matrix V with absolute 00243 value less than eps_tol is replaced by exact zero */ 00244 double max_gro; 00245 /* maximal allowable growth of elements of the matrix V during 00246 all the factorization process; if on some elimination step the 00247 ratio big_v / max_a (see below) becomes greater than max_gro, 00248 the matrix A is considered as ill-conditioned (it is assumed 00249 that the tolerance piv_tol has an adequate value) */ 00250 /*--------------------------------------------------------------*/ 00251 /* some statistics */ 00252 int nnz_a; 00253 /* number of non-zeros in the matrix A */ 00254 int nnz_f; 00255 /* number of non-zeros in the matrix F (except diagonal elements, 00256 which are always equal to one and therefore not stored) */ 00257 int nnz_v; 00258 /* number of non-zeros in the matrix V (except pivot elements, 00259 which correspond to diagonal elements of the matrix U = P*V*Q 00260 and which are stored separately in the array vr_piv) */ 00261 double max_a; 00262 /* largest of absolute values of elements of the matrix A */ 00263 double big_v; 00264 /* estimated largest of absolute values of elements appeared in 00265 the active submatrix during all the factorization process */ 00266 int rank; 00267 /* estimated rank of the matrix A */ 00268 }; 00269 00270 struct LUF_WA 00271 { /* working area (used only during factorization) */ 00272 double *rs_max; /* double rs_max[1+n]; */ 00273 /* rs_max[0] is not used; rs_max[i], 1 <= i <= n, is used only if 00274 the i-th row of the matrix V belongs to the active submatrix 00275 and is the largest of absolute values of elements in this row; 00276 rs_max[i] < 0.0 means that the largest value is not known yet 00277 and should be determined by the pivoting routine */ 00278 /*--------------------------------------------------------------*/ 00279 /* in order to efficiently implement Markowitz strategy and Duff 00280 search technique there are two families {R[0], R[1], ..., R[n]} 00281 and {C[0], C[1], ..., C[n]}; member R[k] is a set of active 00282 rows of the matrix V, which have k non-zeros; similarly, member 00283 C[k] is a set of active columns of the matrix V, which have k 00284 non-zeros (in the active submatrix); each set R[k] and C[k] is 00285 implemented as a separate doubly linked list */ 00286 int *rs_head; /* int rs_head[1+n]; */ 00287 /* rs_head[k], 0 <= k <= n, is number of the first active row, 00288 which has k non-zeros */ 00289 int *rs_prev; /* int rs_prev[1+n]; */ 00290 /* rs_prev[0] is not used; rs_prev[i], 1 <= i <= n, is number of 00291 the previous active row, which has the same number of non-zeros 00292 as the i-th row */ 00293 int *rs_next; /* int rs_next[1+n]; */ 00294 /* rs_next[0] is not used; rs_next[i], 1 <= i <= n, is number of 00295 the next active row, which has the same number of non-zeros as 00296 the i-th row */ 00297 int *cs_head; /* int cs_head[1+n]; */ 00298 /* cs_head[k], 0 <= k <= n, is number of the first active column, 00299 which has k non-zeros (in the active submatrix) */ 00300 int *cs_prev; /* int cs_prev[1+n]; */ 00301 /* cs_prev[0] is not used; cs_prev[j], 1 <= j <= n, is number of 00302 the previous active column, which has the same number of 00303 non-zeros (in the active submatrix) as the j-th column */ 00304 int *cs_next; /* int cs_next[1+n]; */ 00305 /* cs_next[0] is not used; cs_next[j], 1 <= j <= n, is number of 00306 the next active column, which has the same number of non-zeros 00307 (in the active submatrix) as the j-th column */ 00308 }; 00309 00310 LUF *luf_create(int n, int sv_size); 00311 /* create LU-factorization */ 00312 00313 void luf_defrag_sva(LUF *luf); 00314 /* defragment the sparse vector area */ 00315 00316 int luf_enlarge_row(LUF *luf, int i, int cap); 00317 /* enlarge row capacity */ 00318 00319 int luf_enlarge_col(LUF *luf, int j, int cap); 00320 /* enlarge column capacity */ 00321 00322 LUF_WA *luf_alloc_wa(LUF *luf); 00323 /* pre-allocate working area */ 00324 00325 void luf_free_wa(LUF_WA *wa); 00326 /* free working area */ 00327 00328 int luf_decomp(LUF *luf, 00329 void *info, int (*col)(void *info, int j, int rn[], double aj[]), 00330 LUF_WA *wa); 00331 /* compute LU-factorization */ 00332 00333 void luf_f_solve(LUF *luf, int tr, double x[]); 00334 /* solve system F*x = b or F'*x = b */ 00335 00336 void luf_v_solve(LUF *luf, int tr, double x[]); 00337 /* solve system V*x = b or V'*x = b */ 00338 00339 void luf_solve(LUF *luf, int tr, double x[]); 00340 /* solve system A*x = b or A'*x = b */ 00341 00342 void luf_delete(LUF *luf); 00343 /* delete LU-factorization */ 00344 00345 #endif 00346 00347 /* eof */