1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/glpluf.h Mon Dec 06 13:09:21 2010 +0100
1.3 @@ -0,0 +1,326 @@
1.4 +/* glpluf.h (LU-factorization) */
1.5 +
1.6 +/***********************************************************************
1.7 +* This code is part of GLPK (GNU Linear Programming Kit).
1.8 +*
1.9 +* Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
1.10 +* 2009, 2010 Andrew Makhorin, Department for Applied Informatics,
1.11 +* Moscow Aviation Institute, Moscow, Russia. All rights reserved.
1.12 +* E-mail: <mao@gnu.org>.
1.13 +*
1.14 +* GLPK is free software: you can redistribute it and/or modify it
1.15 +* under the terms of the GNU General Public License as published by
1.16 +* the Free Software Foundation, either version 3 of the License, or
1.17 +* (at your option) any later version.
1.18 +*
1.19 +* GLPK is distributed in the hope that it will be useful, but WITHOUT
1.20 +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
1.21 +* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
1.22 +* License for more details.
1.23 +*
1.24 +* You should have received a copy of the GNU General Public License
1.25 +* along with GLPK. If not, see <http://www.gnu.org/licenses/>.
1.26 +***********************************************************************/
1.27 +
1.28 +#ifndef GLPLUF_H
1.29 +#define GLPLUF_H
1.30 +
1.31 +/***********************************************************************
1.32 +* The structure LUF defines LU-factorization of a square matrix A and
1.33 +* is the following quartet:
1.34 +*
1.35 +* [A] = (F, V, P, Q), (1)
1.36 +*
1.37 +* where F and V are such matrices that
1.38 +*
1.39 +* A = F * V, (2)
1.40 +*
1.41 +* and P and Q are such permutation matrices that the matrix
1.42 +*
1.43 +* L = P * F * inv(P) (3)
1.44 +*
1.45 +* is lower triangular with unity diagonal, and the matrix
1.46 +*
1.47 +* U = P * V * Q (4)
1.48 +*
1.49 +* is upper triangular. All the matrices have the order n.
1.50 +*
1.51 +* Matrices F and V are stored in row- and column-wise sparse format
1.52 +* as row and column linked lists of non-zero elements. Unity elements
1.53 +* on the main diagonal of matrix F are not stored. Pivot elements of
1.54 +* matrix V (which correspond to diagonal elements of matrix U) are
1.55 +* stored separately in an ordinary array.
1.56 +*
1.57 +* Permutation matrices P and Q are stored in ordinary arrays in both
1.58 +* row- and column-like formats.
1.59 +*
1.60 +* Matrices L and U are completely defined by matrices F, V, P, and Q
1.61 +* and therefore not stored explicitly.
1.62 +*
1.63 +* The factorization (1)-(4) is a version of LU-factorization. Indeed,
1.64 +* from (3) and (4) it follows that:
1.65 +*
1.66 +* F = inv(P) * L * P,
1.67 +*
1.68 +* U = inv(P) * U * inv(Q),
1.69 +*
1.70 +* and substitution into (2) leads to:
1.71 +*
1.72 +* A = F * V = inv(P) * L * U * inv(Q).
1.73 +*
1.74 +* For more details see the program documentation. */
1.75 +
1.76 +typedef struct LUF LUF;
1.77 +
1.78 +struct LUF
1.79 +{ /* LU-factorization of a square matrix */
1.80 + int n_max;
1.81 + /* maximal value of n (increased automatically, if necessary) */
1.82 + int n;
1.83 + /* the order of matrices A, F, V, P, Q */
1.84 + int valid;
1.85 + /* the factorization is valid only if this flag is set */
1.86 + /*--------------------------------------------------------------*/
1.87 + /* matrix F in row-wise format */
1.88 + int *fr_ptr; /* int fr_ptr[1+n_max]; */
1.89 + /* fr_ptr[i], i = 1,...,n, is a pointer to the first element of
1.90 + i-th row in SVA */
1.91 + int *fr_len; /* int fr_len[1+n_max]; */
1.92 + /* fr_len[i], i = 1,...,n, is the number of elements in i-th row
1.93 + (except unity diagonal element) */
1.94 + /*--------------------------------------------------------------*/
1.95 + /* matrix F in column-wise format */
1.96 + int *fc_ptr; /* int fc_ptr[1+n_max]; */
1.97 + /* fc_ptr[j], j = 1,...,n, is a pointer to the first element of
1.98 + j-th column in SVA */
1.99 + int *fc_len; /* int fc_len[1+n_max]; */
1.100 + /* fc_len[j], j = 1,...,n, is the number of elements in j-th
1.101 + column (except unity diagonal element) */
1.102 + /*--------------------------------------------------------------*/
1.103 + /* matrix V in row-wise format */
1.104 + int *vr_ptr; /* int vr_ptr[1+n_max]; */
1.105 + /* vr_ptr[i], i = 1,...,n, is a pointer to the first element of
1.106 + i-th row in SVA */
1.107 + int *vr_len; /* int vr_len[1+n_max]; */
1.108 + /* vr_len[i], i = 1,...,n, is the number of elements in i-th row
1.109 + (except pivot element) */
1.110 + int *vr_cap; /* int vr_cap[1+n_max]; */
1.111 + /* vr_cap[i], i = 1,...,n, is the capacity of i-th row, i.e.
1.112 + maximal number of elements which can be stored in the row
1.113 + without relocating it, vr_cap[i] >= vr_len[i] */
1.114 + double *vr_piv; /* double vr_piv[1+n_max]; */
1.115 + /* vr_piv[p], p = 1,...,n, is the pivot element v[p,q] which
1.116 + corresponds to a diagonal element of matrix U = P*V*Q */
1.117 + /*--------------------------------------------------------------*/
1.118 + /* matrix V in column-wise format */
1.119 + int *vc_ptr; /* int vc_ptr[1+n_max]; */
1.120 + /* vc_ptr[j], j = 1,...,n, is a pointer to the first element of
1.121 + j-th column in SVA */
1.122 + int *vc_len; /* int vc_len[1+n_max]; */
1.123 + /* vc_len[j], j = 1,...,n, is the number of elements in j-th
1.124 + column (except pivot element) */
1.125 + int *vc_cap; /* int vc_cap[1+n_max]; */
1.126 + /* vc_cap[j], j = 1,...,n, is the capacity of j-th column, i.e.
1.127 + maximal number of elements which can be stored in the column
1.128 + without relocating it, vc_cap[j] >= vc_len[j] */
1.129 + /*--------------------------------------------------------------*/
1.130 + /* matrix P */
1.131 + int *pp_row; /* int pp_row[1+n_max]; */
1.132 + /* pp_row[i] = j means that P[i,j] = 1 */
1.133 + int *pp_col; /* int pp_col[1+n_max]; */
1.134 + /* pp_col[j] = i means that P[i,j] = 1 */
1.135 + /* if i-th row or column of matrix F is i'-th row or column of
1.136 + matrix L, or if i-th row of matrix V is i'-th row of matrix U,
1.137 + then pp_row[i'] = i and pp_col[i] = i' */
1.138 + /*--------------------------------------------------------------*/
1.139 + /* matrix Q */
1.140 + int *qq_row; /* int qq_row[1+n_max]; */
1.141 + /* qq_row[i] = j means that Q[i,j] = 1 */
1.142 + int *qq_col; /* int qq_col[1+n_max]; */
1.143 + /* qq_col[j] = i means that Q[i,j] = 1 */
1.144 + /* if j-th column of matrix V is j'-th column of matrix U, then
1.145 + qq_row[j] = j' and qq_col[j'] = j */
1.146 + /*--------------------------------------------------------------*/
1.147 + /* the Sparse Vector Area (SVA) is a set of locations used to
1.148 + store sparse vectors representing rows and columns of matrices
1.149 + F and V; each location is a doublet (ind, val), where ind is
1.150 + an index, and val is a numerical value of a sparse vector
1.151 + element; in the whole each sparse vector is a set of adjacent
1.152 + locations defined by a pointer to the first element and the
1.153 + number of elements; these pointer and number are stored in the
1.154 + corresponding matrix data structure (see above); the left part
1.155 + of SVA is used to store rows and columns of matrix V, and its
1.156 + right part is used to store rows and columns of matrix F; the
1.157 + middle part of SVA contains free (unused) locations */
1.158 + int sv_size;
1.159 + /* the size of SVA, in locations; all locations are numbered by
1.160 + integers 1, ..., n, and location 0 is not used; if necessary,
1.161 + the SVA size is automatically increased */
1.162 + int sv_beg, sv_end;
1.163 + /* SVA partitioning pointers:
1.164 + locations from 1 to sv_beg-1 belong to the left part
1.165 + locations from sv_beg to sv_end-1 belong to the middle part
1.166 + locations from sv_end to sv_size belong to the right part
1.167 + the size of the middle part is (sv_end - sv_beg) */
1.168 + int *sv_ind; /* sv_ind[1+sv_size]; */
1.169 + /* sv_ind[k], 1 <= k <= sv_size, is the index field of k-th
1.170 + location */
1.171 + double *sv_val; /* sv_val[1+sv_size]; */
1.172 + /* sv_val[k], 1 <= k <= sv_size, is the value field of k-th
1.173 + location */
1.174 + /*--------------------------------------------------------------*/
1.175 + /* in order to efficiently defragment the left part of SVA there
1.176 + is a doubly linked list of rows and columns of matrix V, where
1.177 + rows are numbered by 1, ..., n, while columns are numbered by
1.178 + n+1, ..., n+n, that allows uniquely identifying each row and
1.179 + column of V by only one integer; in this list rows and columns
1.180 + are ordered by ascending their pointers vr_ptr and vc_ptr */
1.181 + int sv_head;
1.182 + /* the number of leftmost row/column */
1.183 + int sv_tail;
1.184 + /* the number of rightmost row/column */
1.185 + int *sv_prev; /* int sv_prev[1+n_max+n_max]; */
1.186 + /* sv_prev[k], k = 1,...,n+n, is the number of a row/column which
1.187 + precedes k-th row/column */
1.188 + int *sv_next; /* int sv_next[1+n_max+n_max]; */
1.189 + /* sv_next[k], k = 1,...,n+n, is the number of a row/column which
1.190 + succedes k-th row/column */
1.191 + /*--------------------------------------------------------------*/
1.192 + /* working segment (used only during factorization) */
1.193 + double *vr_max; /* int vr_max[1+n_max]; */
1.194 + /* vr_max[i], 1 <= i <= n, is used only if i-th row of matrix V
1.195 + is active (i.e. belongs to the active submatrix), and is the
1.196 + largest magnitude of elements in i-th row; if vr_max[i] < 0,
1.197 + the largest magnitude is not known yet and should be computed
1.198 + by the pivoting routine */
1.199 + /*--------------------------------------------------------------*/
1.200 + /* in order to efficiently implement Markowitz strategy and Duff
1.201 + search technique there are two families {R[0], R[1], ..., R[n]}
1.202 + and {C[0], C[1], ..., C[n]}; member R[k] is the set of active
1.203 + rows of matrix V, which have k non-zeros, and member C[k] is
1.204 + the set of active columns of V, which have k non-zeros in the
1.205 + active submatrix (i.e. in the active rows); each set R[k] and
1.206 + C[k] is implemented as a separate doubly linked list */
1.207 + int *rs_head; /* int rs_head[1+n_max]; */
1.208 + /* rs_head[k], 0 <= k <= n, is the number of first active row,
1.209 + which has k non-zeros */
1.210 + int *rs_prev; /* int rs_prev[1+n_max]; */
1.211 + /* rs_prev[i], 1 <= i <= n, is the number of previous row, which
1.212 + has the same number of non-zeros as i-th row */
1.213 + int *rs_next; /* int rs_next[1+n_max]; */
1.214 + /* rs_next[i], 1 <= i <= n, is the number of next row, which has
1.215 + the same number of non-zeros as i-th row */
1.216 + int *cs_head; /* int cs_head[1+n_max]; */
1.217 + /* cs_head[k], 0 <= k <= n, is the number of first active column,
1.218 + which has k non-zeros (in the active rows) */
1.219 + int *cs_prev; /* int cs_prev[1+n_max]; */
1.220 + /* cs_prev[j], 1 <= j <= n, is the number of previous column,
1.221 + which has the same number of non-zeros (in the active rows) as
1.222 + j-th column */
1.223 + int *cs_next; /* int cs_next[1+n_max]; */
1.224 + /* cs_next[j], 1 <= j <= n, is the number of next column, which
1.225 + has the same number of non-zeros (in the active rows) as j-th
1.226 + column */
1.227 + /* (end of working segment) */
1.228 + /*--------------------------------------------------------------*/
1.229 + /* working arrays */
1.230 + int *flag; /* int flag[1+n_max]; */
1.231 + /* integer working array */
1.232 + double *work; /* double work[1+n_max]; */
1.233 + /* floating-point working array */
1.234 + /*--------------------------------------------------------------*/
1.235 + /* control parameters */
1.236 + int new_sva;
1.237 + /* new required size of the sparse vector area, in locations; set
1.238 + automatically by the factorizing routine */
1.239 + double piv_tol;
1.240 + /* threshold pivoting tolerance, 0 < piv_tol < 1; element v[i,j]
1.241 + of the active submatrix fits to be pivot if it satisfies to the
1.242 + stability criterion |v[i,j]| >= piv_tol * max |v[i,*]|, i.e. if
1.243 + it is not very small in the magnitude among other elements in
1.244 + the same row; decreasing this parameter gives better sparsity
1.245 + at the expense of numerical accuracy and vice versa */
1.246 + int piv_lim;
1.247 + /* maximal allowable number of pivot candidates to be considered;
1.248 + if piv_lim pivot candidates have been considered, the pivoting
1.249 + routine terminates the search with the best candidate found */
1.250 + int suhl;
1.251 + /* if this flag is set, the pivoting routine applies a heuristic
1.252 + proposed by Uwe Suhl: if a column of the active submatrix has
1.253 + no eligible pivot candidates (i.e. all its elements do not
1.254 + satisfy to the stability criterion), the routine excludes it
1.255 + from futher consideration until it becomes column singleton;
1.256 + in many cases this allows reducing the time needed for pivot
1.257 + searching */
1.258 + double eps_tol;
1.259 + /* epsilon tolerance; each element of the active submatrix, whose
1.260 + magnitude is less than eps_tol, is replaced by exact zero */
1.261 + double max_gro;
1.262 + /* maximal allowable growth of elements of matrix V during all
1.263 + the factorization process; if on some eliminaion step the ratio
1.264 + big_v / max_a (see below) becomes greater than max_gro, matrix
1.265 + A is considered as ill-conditioned (assuming that the pivoting
1.266 + tolerance piv_tol has an appropriate value) */
1.267 + /*--------------------------------------------------------------*/
1.268 + /* some statistics */
1.269 + int nnz_a;
1.270 + /* the number of non-zeros in matrix A */
1.271 + int nnz_f;
1.272 + /* the number of non-zeros in matrix F (except diagonal elements,
1.273 + which are not stored) */
1.274 + int nnz_v;
1.275 + /* the number of non-zeros in matrix V (except its pivot elements,
1.276 + which are stored in a separate array) */
1.277 + double max_a;
1.278 + /* the largest magnitude of elements of matrix A */
1.279 + double big_v;
1.280 + /* the largest magnitude of elements of matrix V appeared in the
1.281 + active submatrix during all the factorization process */
1.282 + int rank;
1.283 + /* estimated rank of matrix A */
1.284 +};
1.285 +
1.286 +/* return codes: */
1.287 +#define LUF_ESING 1 /* singular matrix */
1.288 +#define LUF_ECOND 2 /* ill-conditioned matrix */
1.289 +
1.290 +#define luf_create_it _glp_luf_create_it
1.291 +LUF *luf_create_it(void);
1.292 +/* create LU-factorization */
1.293 +
1.294 +#define luf_defrag_sva _glp_luf_defrag_sva
1.295 +void luf_defrag_sva(LUF *luf);
1.296 +/* defragment the sparse vector area */
1.297 +
1.298 +#define luf_enlarge_row _glp_luf_enlarge_row
1.299 +int luf_enlarge_row(LUF *luf, int i, int cap);
1.300 +/* enlarge row capacity */
1.301 +
1.302 +#define luf_enlarge_col _glp_luf_enlarge_col
1.303 +int luf_enlarge_col(LUF *luf, int j, int cap);
1.304 +/* enlarge column capacity */
1.305 +
1.306 +#define luf_factorize _glp_luf_factorize
1.307 +int luf_factorize(LUF *luf, int n, int (*col)(void *info, int j,
1.308 + int ind[], double val[]), void *info);
1.309 +/* compute LU-factorization */
1.310 +
1.311 +#define luf_f_solve _glp_luf_f_solve
1.312 +void luf_f_solve(LUF *luf, int tr, double x[]);
1.313 +/* solve system F*x = b or F'*x = b */
1.314 +
1.315 +#define luf_v_solve _glp_luf_v_solve
1.316 +void luf_v_solve(LUF *luf, int tr, double x[]);
1.317 +/* solve system V*x = b or V'*x = b */
1.318 +
1.319 +#define luf_a_solve _glp_luf_a_solve
1.320 +void luf_a_solve(LUF *luf, int tr, double x[]);
1.321 +/* solve system A*x = b or A'*x = b */
1.322 +
1.323 +#define luf_delete_it _glp_luf_delete_it
1.324 +void luf_delete_it(LUF *luf);
1.325 +/* delete LU-factorization */
1.326 +
1.327 +#endif
1.328 +
1.329 +/* eof */