1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/deps/glpk/src/glpluf.h	Sun Nov 06 20:59:10 2011 +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, 2011 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 */