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