DyLP trunk
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$
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 Files Functions Variables Typedefs Enumerations Enumerator Friends Defines