# glpluf.h

Go to the documentation of this file.

    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 Eclipse Public License (EPL) 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 407 2010-12-31 20:48:48Z 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 */

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines