alpar@1: /* glpluf.h (LU-factorization) */
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  This code is part of GLPK (GNU Linear Programming Kit).
alpar@1: *
alpar@1: *  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
alpar@1: *  2009, 2010 Andrew Makhorin, Department for Applied Informatics,
alpar@1: *  Moscow Aviation Institute, Moscow, Russia. All rights reserved.
alpar@1: *  E-mail: <mao@gnu.org>.
alpar@1: *
alpar@1: *  GLPK is free software: you can redistribute it and/or modify it
alpar@1: *  under the terms of the GNU General Public License as published by
alpar@1: *  the Free Software Foundation, either version 3 of the License, or
alpar@1: *  (at your option) any later version.
alpar@1: *
alpar@1: *  GLPK is distributed in the hope that it will be useful, but WITHOUT
alpar@1: *  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
alpar@1: *  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
alpar@1: *  License for more details.
alpar@1: *
alpar@1: *  You should have received a copy of the GNU General Public License
alpar@1: *  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
alpar@1: ***********************************************************************/
alpar@1: 
alpar@1: #ifndef GLPLUF_H
alpar@1: #define GLPLUF_H
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  The structure LUF defines LU-factorization of a square matrix A and
alpar@1: *  is the following quartet:
alpar@1: *
alpar@1: *     [A] = (F, V, P, Q),                                            (1)
alpar@1: *
alpar@1: *  where F and V are such matrices that
alpar@1: *
alpar@1: *     A = F * V,                                                     (2)
alpar@1: *
alpar@1: *  and P and Q are such permutation matrices that the matrix
alpar@1: *
alpar@1: *     L = P * F * inv(P)                                             (3)
alpar@1: *
alpar@1: *  is lower triangular with unity diagonal, and the matrix
alpar@1: *
alpar@1: *     U = P * V * Q                                                  (4)
alpar@1: *
alpar@1: *  is upper triangular. All the matrices have the order n.
alpar@1: *
alpar@1: *  Matrices F and V are stored in row- and column-wise sparse format
alpar@1: *  as row and column linked lists of non-zero elements. Unity elements
alpar@1: *  on the main diagonal of matrix F are not stored. Pivot elements of
alpar@1: *  matrix V (which correspond to diagonal elements of matrix U) are
alpar@1: *  stored separately in an ordinary array.
alpar@1: *
alpar@1: *  Permutation matrices P and Q are stored in ordinary arrays in both
alpar@1: *  row- and column-like formats.
alpar@1: *
alpar@1: *  Matrices L and U are completely defined by matrices F, V, P, and Q
alpar@1: *  and therefore not stored explicitly.
alpar@1: *
alpar@1: *  The factorization (1)-(4) is a version of LU-factorization. Indeed,
alpar@1: *  from (3) and (4) it follows that:
alpar@1: *
alpar@1: *     F = inv(P) * L * P,
alpar@1: *
alpar@1: *     U = inv(P) * U * inv(Q),
alpar@1: *
alpar@1: *  and substitution into (2) leads to:
alpar@1: *
alpar@1: *     A = F * V = inv(P) * L * U * inv(Q).
alpar@1: *
alpar@1: *  For more details see the program documentation. */
alpar@1: 
alpar@1: typedef struct LUF LUF;
alpar@1: 
alpar@1: struct LUF
alpar@1: {     /* LU-factorization of a square matrix */
alpar@1:       int n_max;
alpar@1:       /* maximal value of n (increased automatically, if necessary) */
alpar@1:       int n;
alpar@1:       /* the order of matrices A, F, V, P, Q */
alpar@1:       int valid;
alpar@1:       /* the factorization is valid only if this flag is set */
alpar@1:       /*--------------------------------------------------------------*/
alpar@1:       /* matrix F in row-wise format */
alpar@1:       int *fr_ptr; /* int fr_ptr[1+n_max]; */
alpar@1:       /* fr_ptr[i], i = 1,...,n, is a pointer to the first element of
alpar@1:          i-th row in SVA */
alpar@1:       int *fr_len; /* int fr_len[1+n_max]; */
alpar@1:       /* fr_len[i], i = 1,...,n, is the number of elements in i-th row
alpar@1:          (except unity diagonal element) */
alpar@1:       /*--------------------------------------------------------------*/
alpar@1:       /* matrix F in column-wise format */
alpar@1:       int *fc_ptr; /* int fc_ptr[1+n_max]; */
alpar@1:       /* fc_ptr[j], j = 1,...,n, is a pointer to the first element of
alpar@1:          j-th column in SVA */
alpar@1:       int *fc_len; /* int fc_len[1+n_max]; */
alpar@1:       /* fc_len[j], j = 1,...,n, is the number of elements in j-th
alpar@1:          column (except unity diagonal element) */
alpar@1:       /*--------------------------------------------------------------*/
alpar@1:       /* matrix V in row-wise format */
alpar@1:       int *vr_ptr; /* int vr_ptr[1+n_max]; */
alpar@1:       /* vr_ptr[i], i = 1,...,n, is a pointer to the first element of
alpar@1:          i-th row in SVA */
alpar@1:       int *vr_len; /* int vr_len[1+n_max]; */
alpar@1:       /* vr_len[i], i = 1,...,n, is the number of elements in i-th row
alpar@1:          (except pivot element) */
alpar@1:       int *vr_cap; /* int vr_cap[1+n_max]; */
alpar@1:       /* vr_cap[i], i = 1,...,n, is the capacity of i-th row, i.e.
alpar@1:          maximal number of elements which can be stored in the row
alpar@1:          without relocating it, vr_cap[i] >= vr_len[i] */
alpar@1:       double *vr_piv; /* double vr_piv[1+n_max]; */
alpar@1:       /* vr_piv[p], p = 1,...,n, is the pivot element v[p,q] which
alpar@1:          corresponds to a diagonal element of matrix U = P*V*Q */
alpar@1:       /*--------------------------------------------------------------*/
alpar@1:       /* matrix V in column-wise format */
alpar@1:       int *vc_ptr; /* int vc_ptr[1+n_max]; */
alpar@1:       /* vc_ptr[j], j = 1,...,n, is a pointer to the first element of
alpar@1:          j-th column in SVA */
alpar@1:       int *vc_len; /* int vc_len[1+n_max]; */
alpar@1:       /* vc_len[j], j = 1,...,n, is the number of elements in j-th
alpar@1:          column (except pivot element) */
alpar@1:       int *vc_cap; /* int vc_cap[1+n_max]; */
alpar@1:       /* vc_cap[j], j = 1,...,n, is the capacity of j-th column, i.e.
alpar@1:          maximal number of elements which can be stored in the column
alpar@1:          without relocating it, vc_cap[j] >= vc_len[j] */
alpar@1:       /*--------------------------------------------------------------*/
alpar@1:       /* matrix P */
alpar@1:       int *pp_row; /* int pp_row[1+n_max]; */
alpar@1:       /* pp_row[i] = j means that P[i,j] = 1 */
alpar@1:       int *pp_col; /* int pp_col[1+n_max]; */
alpar@1:       /* pp_col[j] = i means that P[i,j] = 1 */
alpar@1:       /* if i-th row or column of matrix F is i'-th row or column of
alpar@1:          matrix L, or if i-th row of matrix V is i'-th row of matrix U,
alpar@1:          then pp_row[i'] = i and pp_col[i] = i' */
alpar@1:       /*--------------------------------------------------------------*/
alpar@1:       /* matrix Q */
alpar@1:       int *qq_row; /* int qq_row[1+n_max]; */
alpar@1:       /* qq_row[i] = j means that Q[i,j] = 1 */
alpar@1:       int *qq_col; /* int qq_col[1+n_max]; */
alpar@1:       /* qq_col[j] = i means that Q[i,j] = 1 */
alpar@1:       /* if j-th column of matrix V is j'-th column of matrix U, then
alpar@1:          qq_row[j] = j' and qq_col[j'] = j */
alpar@1:       /*--------------------------------------------------------------*/
alpar@1:       /* the Sparse Vector Area (SVA) is a set of locations used to
alpar@1:          store sparse vectors representing rows and columns of matrices
alpar@1:          F and V; each location is a doublet (ind, val), where ind is
alpar@1:          an index, and val is a numerical value of a sparse vector
alpar@1:          element; in the whole each sparse vector is a set of adjacent
alpar@1:          locations defined by a pointer to the first element and the
alpar@1:          number of elements; these pointer and number are stored in the
alpar@1:          corresponding matrix data structure (see above); the left part
alpar@1:          of SVA is used to store rows and columns of matrix V, and its
alpar@1:          right part is used to store rows and columns of matrix F; the
alpar@1:          middle part of SVA contains free (unused) locations */
alpar@1:       int sv_size;
alpar@1:       /* the size of SVA, in locations; all locations are numbered by
alpar@1:          integers 1, ..., n, and location 0 is not used; if necessary,
alpar@1:          the SVA size is automatically increased */
alpar@1:       int sv_beg, sv_end;
alpar@1:       /* SVA partitioning pointers:
alpar@1:          locations from 1 to sv_beg-1 belong to the left part
alpar@1:          locations from sv_beg to sv_end-1 belong to the middle part
alpar@1:          locations from sv_end to sv_size belong to the right part
alpar@1:          the size of the middle part is (sv_end - sv_beg) */
alpar@1:       int *sv_ind; /* sv_ind[1+sv_size]; */
alpar@1:       /* sv_ind[k], 1 <= k <= sv_size, is the index field of k-th
alpar@1:          location */
alpar@1:       double *sv_val; /* sv_val[1+sv_size]; */
alpar@1:       /* sv_val[k], 1 <= k <= sv_size, is the value field of k-th
alpar@1:          location */
alpar@1:       /*--------------------------------------------------------------*/
alpar@1:       /* in order to efficiently defragment the left part of SVA there
alpar@1:          is a doubly linked list of rows and columns of matrix V, where
alpar@1:          rows are numbered by 1, ..., n, while columns are numbered by
alpar@1:          n+1, ..., n+n, that allows uniquely identifying each row and
alpar@1:          column of V by only one integer; in this list rows and columns
alpar@1:          are ordered by ascending their pointers vr_ptr and vc_ptr */
alpar@1:       int sv_head;
alpar@1:       /* the number of leftmost row/column */
alpar@1:       int sv_tail;
alpar@1:       /* the number of rightmost row/column */
alpar@1:       int *sv_prev; /* int sv_prev[1+n_max+n_max]; */
alpar@1:       /* sv_prev[k], k = 1,...,n+n, is the number of a row/column which
alpar@1:          precedes k-th row/column */
alpar@1:       int *sv_next; /* int sv_next[1+n_max+n_max]; */
alpar@1:       /* sv_next[k], k = 1,...,n+n, is the number of a row/column which
alpar@1:          succedes k-th row/column */
alpar@1:       /*--------------------------------------------------------------*/
alpar@1:       /* working segment (used only during factorization) */
alpar@1:       double *vr_max; /* int vr_max[1+n_max]; */
alpar@1:       /* vr_max[i], 1 <= i <= n, is used only if i-th row of matrix V
alpar@1:          is active (i.e. belongs to the active submatrix), and is the
alpar@1:          largest magnitude of elements in i-th row; if vr_max[i] < 0,
alpar@1:          the largest magnitude is not known yet and should be computed
alpar@1:          by the pivoting routine */
alpar@1:       /*--------------------------------------------------------------*/
alpar@1:       /* in order to efficiently implement Markowitz strategy and Duff
alpar@1:          search technique there are two families {R[0], R[1], ..., R[n]}
alpar@1:          and {C[0], C[1], ..., C[n]}; member R[k] is the set of active
alpar@1:          rows of matrix V, which have k non-zeros, and member C[k] is
alpar@1:          the set of active columns of V, which have k non-zeros in the
alpar@1:          active submatrix (i.e. in the active rows); each set R[k] and
alpar@1:          C[k] is implemented as a separate doubly linked list */
alpar@1:       int *rs_head; /* int rs_head[1+n_max]; */
alpar@1:       /* rs_head[k], 0 <= k <= n, is the number of first active row,
alpar@1:          which has k non-zeros */
alpar@1:       int *rs_prev; /* int rs_prev[1+n_max]; */
alpar@1:       /* rs_prev[i], 1 <= i <= n, is the number of previous row, which
alpar@1:          has the same number of non-zeros as i-th row */
alpar@1:       int *rs_next; /* int rs_next[1+n_max]; */
alpar@1:       /* rs_next[i], 1 <= i <= n, is the number of next row, which has
alpar@1:          the same number of non-zeros as i-th row */
alpar@1:       int *cs_head; /* int cs_head[1+n_max]; */
alpar@1:       /* cs_head[k], 0 <= k <= n, is the number of first active column,
alpar@1:          which has k non-zeros (in the active rows) */
alpar@1:       int *cs_prev; /* int cs_prev[1+n_max]; */
alpar@1:       /* cs_prev[j], 1 <= j <= n, is the number of previous column,
alpar@1:          which has the same number of non-zeros (in the active rows) as
alpar@1:          j-th column */
alpar@1:       int *cs_next; /* int cs_next[1+n_max]; */
alpar@1:       /* cs_next[j], 1 <= j <= n, is the number of next column, which
alpar@1:          has the same number of non-zeros (in the active rows) as j-th
alpar@1:          column */
alpar@1:       /* (end of working segment) */
alpar@1:       /*--------------------------------------------------------------*/
alpar@1:       /* working arrays */
alpar@1:       int *flag; /* int flag[1+n_max]; */
alpar@1:       /* integer working array */
alpar@1:       double *work; /* double work[1+n_max]; */
alpar@1:       /* floating-point working array */
alpar@1:       /*--------------------------------------------------------------*/
alpar@1:       /* control parameters */
alpar@1:       int new_sva;
alpar@1:       /* new required size of the sparse vector area, in locations; set
alpar@1:          automatically by the factorizing routine */
alpar@1:       double piv_tol;
alpar@1:       /* threshold pivoting tolerance, 0 < piv_tol < 1; element v[i,j]
alpar@1:          of the active submatrix fits to be pivot if it satisfies to the
alpar@1:          stability criterion |v[i,j]| >= piv_tol * max |v[i,*]|, i.e. if
alpar@1:          it is not very small in the magnitude among other elements in
alpar@1:          the same row; decreasing this parameter gives better sparsity
alpar@1:          at the expense of numerical accuracy and vice versa */
alpar@1:       int piv_lim;
alpar@1:       /* maximal allowable number of pivot candidates to be considered;
alpar@1:          if piv_lim pivot candidates have been considered, the pivoting
alpar@1:          routine terminates the search with the best candidate found */
alpar@1:       int suhl;
alpar@1:       /* if this flag is set, the pivoting routine applies a heuristic
alpar@1:          proposed by Uwe Suhl: if a column of the active submatrix has
alpar@1:          no eligible pivot candidates (i.e. all its elements do not
alpar@1:          satisfy to the stability criterion), the routine excludes it
alpar@1:          from futher consideration until it becomes column singleton;
alpar@1:          in many cases this allows reducing the time needed for pivot
alpar@1:          searching */
alpar@1:       double eps_tol;
alpar@1:       /* epsilon tolerance; each element of the active submatrix, whose
alpar@1:          magnitude is less than eps_tol, is replaced by exact zero */
alpar@1:       double max_gro;
alpar@1:       /* maximal allowable growth of elements of matrix V during all
alpar@1:          the factorization process; if on some eliminaion step the ratio
alpar@1:          big_v / max_a (see below) becomes greater than max_gro, matrix
alpar@1:          A is considered as ill-conditioned (assuming that the pivoting
alpar@1:          tolerance piv_tol has an appropriate value) */
alpar@1:       /*--------------------------------------------------------------*/
alpar@1:       /* some statistics */
alpar@1:       int nnz_a;
alpar@1:       /* the number of non-zeros in matrix A */
alpar@1:       int nnz_f;
alpar@1:       /* the number of non-zeros in matrix F (except diagonal elements,
alpar@1:          which are not stored) */
alpar@1:       int nnz_v;
alpar@1:       /* the number of non-zeros in matrix V (except its pivot elements,
alpar@1:          which are stored in a separate array) */
alpar@1:       double max_a;
alpar@1:       /* the largest magnitude of elements of matrix A */
alpar@1:       double big_v;
alpar@1:       /* the largest magnitude of elements of matrix V appeared in the
alpar@1:          active submatrix during all the factorization process */
alpar@1:       int rank;
alpar@1:       /* estimated rank of matrix A */
alpar@1: };
alpar@1: 
alpar@1: /* return codes: */
alpar@1: #define LUF_ESING    1  /* singular matrix */
alpar@1: #define LUF_ECOND    2  /* ill-conditioned matrix */
alpar@1: 
alpar@1: #define luf_create_it _glp_luf_create_it
alpar@1: LUF *luf_create_it(void);
alpar@1: /* create LU-factorization */
alpar@1: 
alpar@1: #define luf_defrag_sva _glp_luf_defrag_sva
alpar@1: void luf_defrag_sva(LUF *luf);
alpar@1: /* defragment the sparse vector area */
alpar@1: 
alpar@1: #define luf_enlarge_row _glp_luf_enlarge_row
alpar@1: int luf_enlarge_row(LUF *luf, int i, int cap);
alpar@1: /* enlarge row capacity */
alpar@1: 
alpar@1: #define luf_enlarge_col _glp_luf_enlarge_col
alpar@1: int luf_enlarge_col(LUF *luf, int j, int cap);
alpar@1: /* enlarge column capacity */
alpar@1: 
alpar@1: #define luf_factorize _glp_luf_factorize
alpar@1: int luf_factorize(LUF *luf, int n, int (*col)(void *info, int j,
alpar@1:       int ind[], double val[]), void *info);
alpar@1: /* compute LU-factorization */
alpar@1: 
alpar@1: #define luf_f_solve _glp_luf_f_solve
alpar@1: void luf_f_solve(LUF *luf, int tr, double x[]);
alpar@1: /* solve system F*x = b or F'*x = b */
alpar@1: 
alpar@1: #define luf_v_solve _glp_luf_v_solve
alpar@1: void luf_v_solve(LUF *luf, int tr, double x[]);
alpar@1: /* solve system V*x = b or V'*x = b */
alpar@1: 
alpar@1: #define luf_a_solve _glp_luf_a_solve
alpar@1: void luf_a_solve(LUF *luf, int tr, double x[]);
alpar@1: /* solve system A*x = b or A'*x = b */
alpar@1: 
alpar@1: #define luf_delete_it _glp_luf_delete_it
alpar@1: void luf_delete_it(LUF *luf);
alpar@1: /* delete LU-factorization */
alpar@1: 
alpar@1: #endif
alpar@1: 
alpar@1: /* eof */