/* glpluf.h */ /*---------------------------------------------------------------------- -- Copyright (C) 2000, 2001, 2002, 2003 Andrew Makhorin, Department -- for Applied Informatics, Moscow Aviation Institute, Moscow, Russia. -- All rights reserved. E-mail: . -- -- This file is part of GLPK (GNU Linear Programming Kit). -- -- GLPK is free software; you can redistribute it and/or modify it -- under the terms of the GNU General Public License as published by -- the Free Software Foundation; either version 2, or (at your option) -- any later version. -- -- GLPK is distributed in the hope that it will be useful, but WITHOUT -- ANY WARRANTY; without even the implied warranty of MERCHANTABILITY -- or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public -- License for more details. -- -- You should have received a copy of the GNU General Public License -- along with GLPK; see the file COPYING. If not, write to the Free -- Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA -- 02111-1307, USA. ----------------------------------------------------------------------*/ #ifndef _GLPLUF_H #define _GLPLUF_H #include "gnumeric-config.h" #include "gnumeric.h" #include "numbers.h" #define luf_create glp_luf_create #define luf_defrag_sva glp_luf_defrag_sva #define luf_enlarge_row glp_luf_enlarge_row #define luf_enlarge_col glp_luf_enlarge_col #define luf_alloc_wa glp_luf_alloc_wa #define luf_free_wa glp_luf_free_wa #define luf_decomp glp_luf_decomp #define luf_f_solve glp_luf_f_solve #define luf_v_solve glp_luf_v_solve #define luf_solve glp_luf_solve #define luf_delete glp_luf_delete /*---------------------------------------------------------------------- -- The structure LUF defines LU-factorization of a square matrix A, -- which is the following quartet: -- -- [A] = (F, V, P, Q), (1) -- -- where F and V are such matrices that -- -- A = F * V, (2) -- -- and P and Q are such permutation matrices that the matrix -- -- L = P * F * inv(P) (3) -- -- is lower triangular with unity diagonal, and the matrix -- -- U = P * V * Q (4) -- -- is upper triangular. All the matrices have the order n. -- -- The matrices F and V are stored in row/column-wise sparse format as -- row and column linked lists of non-zero elements. Unity elements on -- the main diagonal of the matrix F are not stored. Pivot elements of -- the matrix V (that correspond to diagonal elements of the matrix U) -- are also missing from the row and column lists and stored separately -- in an ordinary array. -- -- The permutation matrices P and Q are stored as ordinary arrays using -- both row- and column-like formats. -- -- The matrices L and U being completely defined by the matrices F, V, -- P, and Q are not stored explicitly. -- -- It can easily be shown that the factorization (1)-(3) is a version of -- LU-factorization. Indeed, from (3) and (4) it follows that -- -- F = inv(P) * L * P, -- -- V = inv(P) * U * inv(Q), -- -- and substitution into (2) gives -- -- A = F * V = inv(P) * L * U * inv(Q). -- -- For more details see the program documentation. */ typedef struct LUF LUF; typedef struct LUF_WA LUF_WA; struct LUF { /* LU-factorization of a square matrix */ int n; /* order of the matrices A, F, V, P, Q */ int valid; /* if this flag is not set, the factorization is invalid */ /*--------------------------------------------------------------*/ /* matrix F in row-wise format */ int *fr_ptr; /* int fr_ptr[1+n]; */ /* fr_ptr[0] is not used; fr_ptr[i], i = 1, ..., n, is a pointer to the first element of the i-th row in the sparse vector area */ int *fr_len; /* int fr_len[1+n]; */ /* fr_len[0] is not used; fr_len[i], i = 1, ..., n, is number of elements in the i-th row (except unity diagonal element) */ /*--------------------------------------------------------------*/ /* matrix F in column-wise format */ int *fc_ptr; /* int fc_ptr[1+n]; */ /* fc_ptr[0] is not used; fc_ptr[j], j = 1, ..., n, is a pointer to the first element of the j-th column in the sparse vector area */ int *fc_len; /* int fc_len[1+n]; */ /* fc_len[0] is not used; fc_len[j], j = 1, ..., n, is number of elements in the j-th column (except unity diagonal element) */ /*--------------------------------------------------------------*/ /* matrix V in row-wise format */ int *vr_ptr; /* int vr_ptr[1+n]; */ /* vr_ptr[0] is not used; vr_ptr[i], i = 1, ..., n, is a pointer to the first element of the i-th row in the sparse vector area */ int *vr_len; /* int vr_len[1+n]; */ /* vr_len[0] is not used; vr_len[i], i = 1, ..., n, is number of elements in the i-th row (except pivot element) */ int *vr_cap; /* int vr_cap[1+n]; */ /* vr_cap[0] is not used; vr_cap[i], i = 1, ..., n, is capacity of the i-th row, i.e. maximal number of elements, which can be stored there without relocating the row, vr_cap[i] >= vr_len[i] */ gnm_float *vr_piv; /* gnm_float vr_piv[1+n]; */ /* vr_piv[0] is not used; vr_piv[p], p = 1, ..., n, is the pivot element v[p,q], which corresponds to a diagonal element of the matrix U = P*V*Q */ /*--------------------------------------------------------------*/ /* matrix V in column-wise format */ int *vc_ptr; /* int vc_ptr[1+n]; */ /* vc_ptr[0] is not used; vc_ptr[j], j = 1, ..., n, is a pointer to the first element of the j-th column in the sparse vector area */ int *vc_len; /* int vc_len[1+n]; */ /* vc_len[0] is not used; vc_len[j], j = 1, ..., n, is number of elements in the j-th column (except pivot element) */ int *vc_cap; /* int vc_cap[1+n]; */ /* vc_cap[0] is not used; vc_cap[j], j = 1, ..., n, is capacity of the j-th column, i.e. maximal number of elements, which can be stored there without relocating the column, vc_cap[j] >= vc_len[j] */ /*--------------------------------------------------------------*/ /* matrix P */ int *pp_row; /* int pp_row[1+n]; */ /* pp_row[0] is not used; pp_row[i] = j means that p[i,j] = 1 */ int *pp_col; /* int pp_col[1+n]; */ /* pp_col[0] is not used; pp_col[j] = i means that p[i,j] = 1 */ /* if i-th row or column of the matrix F corresponds to i'-th row or column of the matrix L = P*F*inv(P), or if i-th row of the matrix V corresponds to i'-th row of the matrix U = P*V*Q, then pp_row[i'] = i and pp_col[i] = i' */ /*--------------------------------------------------------------*/ /* matrix Q */ int *qq_row; /* int qq_row[1+n]; */ /* qq_row[0] is not used; qq_row[i] = j means that q[i,j] = 1 */ int *qq_col; /* int qq_col[1+n]; */ /* qq_col[0] is not used; qq_col[j] = i means that q[i,j] = 1 */ /* if j-th column of the matrix V corresponds to j'-th column of the matrix U = P*V*Q, then qq_row[j] = j' and qq_col[j'] = j */ /*--------------------------------------------------------------*/ /* sparse vector area (SVA) is a set of locations intended to store sparse vectors that represent rows and columns of the matrices F and V; each location is the gnm_floatt (ndx, val), where ndx is an index and val is a numerical value of a sparse vector element; in the whole each sparse vector is a set of adjacent locations defined by a pointer to the first element and number of elements; these pointer and number are stored in the corresponding matrix data structure (see above); the left part of SVA is used to store rows and columns of the matrix V, the right part is used to store rows and columns of the matrix F; between the left and right parts there is the middle part, locations of which are free */ int sv_size; /* total size of the sparse vector area, in locations; locations are numbered by integers 1, 2, ..., sv_size, and location with the number 0 is not used; if it is necessary, the size of SVA is automatically increased */ int sv_beg, sv_end; /* SVA partitioning pointers: locations 1, ..., sv_beg-1 belong to the left part; locations sv_beg, ..., sv_end-1 belong to the middle part; locations sv_end, ..., sv_size belong to the right part; number of free locations, i.e. locations that belong to the middle part, is (sv_end - sv_beg) */ int *sv_ndx; /* int sv_ndx[1+sv_size]; */ /* sv_ndx[0] is not used; sv_ndx[k], 1 <= k <= sv_size, is the index field of the k-th location */ gnm_float *sv_val; /* gnm_float sv_val[1+sv_size]; */ /* sv_val[0] is not used; sv_val[k], 1 <= k <= sv_size, is the value field of the k-th location */ /* in order to efficiently defragment the left part of SVA there is a gnm_float linked list of rows and columns of the matrix V, where rows have numbers 1, ..., n, and columns have numbers n+1, ..., n+n, due to that each row and column can be uniquely identified by one integer; in this list rows and columns are ordered by ascending their pointers vr_ptr[i] and vc_ptr[j] */ int sv_head; /* the number of the leftmost row/column */ int sv_tail; /* the number of the rightmost row/column */ int *sv_prev; /* int sv_prev[1+n+n]; */ /* sv_prev[k], k = 1, ..., n+n, is the number of a row/column, which precedes the k-th row/column */ int *sv_next; /* int sv_next[1+n+n]; */ /* sv_next[k], k = 1, ..., n+n, is the number of a row/column, which succedes the k-th row/column */ /*--------------------------------------------------------------*/ /* working arrays */ int *flag; /* int flag[1+n]; */ /* integer working array */ gnm_float *work; /* gnm_float work[1+n]; */ /* floating-point working array */ /*--------------------------------------------------------------*/ /* control parameters */ int new_sva; /* new required size of the sparse vector area, in locations; set automatically by the factorizing routine */ gnm_float piv_tol; /* threshold pivoting tolerance, 0 < piv_tol < 1; element v[i,j] of the active submatrix fits to be pivot if it satisfies to the stability condition |v[i,j]| >= piv_tol * max|v[i,*]|, i.e. if this element is not very small (in absolute value) among other elements in the same row; decreasing this parameter involves better sparsity at the expense of numerical accuracy and vice versa */ int piv_lim; /* maximal allowable number of pivot candidates to be considered; if piv_lim pivot candidates have been considered, the pivoting routine terminates the search with the best candidate found */ int suhl; /* if this flag is set, the pivoting routine applies a heuristic rule proposed by Uwe Suhl: if a column of the active submatrix has no eligible pivot candidates (i.e. all its elements don't satisfy to the stability condition), the routine excludes such column from the futher consideration until it becomes a column singleton; in many cases this reduces a time needed for pivot searching */ gnm_float eps_tol; /* epsilon tolerance; each element of the matrix V with absolute value less than eps_tol is replaced by exact zero */ gnm_float max_gro; /* maximal allowable growth of elements of the matrix V during all the factorization process; if on some elimination step the ratio big_v / max_a (see below) becomes greater than max_gro, the matrix A is considered as ill-conditioned (it is assumed that the tolerance piv_tol has an adequate value) */ /*--------------------------------------------------------------*/ /* some statistics */ int nnz_a; /* number of non-zeros in the matrix A */ int nnz_f; /* number of non-zeros in the matrix F (except diagonal elements, which are always equal to one and therefore not stored) */ int nnz_v; /* number of non-zeros in the matrix V (except pivot elements, which correspond to diagonal elements of the matrix U = P*V*Q and which are stored separately in the array vr_piv) */ gnm_float max_a; /* largest of absolute values of elements of the matrix A */ gnm_float big_v; /* estimated largest of absolute values of elements appeared in the active submatrix during all the factorization process */ int rank; /* estimated rank of the matrix A */ }; struct LUF_WA { /* working area (used only during factorization) */ gnm_float *rs_max; /* gnm_float rs_max[1+n]; */ /* rs_max[0] is not used; rs_max[i], 1 <= i <= n, is used only if the i-th row of the matrix V belongs to the active submatrix and is the largest of absolute values of elements in this row; rs_max[i] < 0.0 means that the largest value is not known yet and should be determined by the pivoting routine */ /*--------------------------------------------------------------*/ /* in order to efficiently implement Markowitz strategy and Duff search technique there are two families {R[0], R[1], ..., R[n]} and {C[0], C[1], ..., C[n]}; member R[k] is a set of active rows of the matrix V, which have k non-zeros; similarly, member C[k] is a set of active columns of the matrix V, which have k non-zeros (in the active submatrix); each set R[k] and C[k] is implemented as a separate doubly linked list */ int *rs_head; /* int rs_head[1+n]; */ /* rs_head[k], 0 <= k <= n, is number of the first active row, which has k non-zeros */ int *rs_prev; /* int rs_prev[1+n]; */ /* rs_prev[0] is not used; rs_prev[i], 1 <= i <= n, is number of the previous active row, which has the same number of non-zeros as the i-th row */ int *rs_next; /* int rs_next[1+n]; */ /* rs_next[0] is not used; rs_next[i], 1 <= i <= n, is number of the next active row, which has the same number of non-zeros as the i-th row */ int *cs_head; /* int cs_head[1+n]; */ /* cs_head[k], 0 <= k <= n, is number of the first active column, which has k non-zeros (in the active submatrix) */ int *cs_prev; /* int cs_prev[1+n]; */ /* cs_prev[0] is not used; cs_prev[j], 1 <= j <= n, is number of the previous active column, which has the same number of non-zeros (in the active submatrix) as the j-th column */ int *cs_next; /* int cs_next[1+n]; */ /* cs_next[0] is not used; cs_next[j], 1 <= j <= n, is number of the next active column, which has the same number of non-zeros (in the active submatrix) as the j-th column */ }; LUF *luf_create(int n, int sv_size); /* create LU-factorization */ void luf_defrag_sva(LUF *luf); /* defragment the sparse vector area */ int luf_enlarge_row(LUF *luf, int i, int cap); /* enlarge row capacity */ int luf_enlarge_col(LUF *luf, int j, int cap); /* enlarge column capacity */ LUF_WA *luf_alloc_wa(LUF *luf); /* pre-allocate working area */ void luf_free_wa(LUF_WA *wa); /* free working area */ int luf_decomp(LUF *luf, void *info, int (*col)(void *info, int j, int rn[], gnm_float aj[]), LUF_WA *wa); /* compute LU-factorization */ void luf_f_solve(LUF *luf, int tr, gnm_float x[]); /* solve system F*x = b or F'*x = b */ void luf_v_solve(LUF *luf, int tr, gnm_float x[]); /* solve system V*x = b or V'*x = b */ void luf_solve(LUF *luf, int tr, gnm_float x[]); /* solve system A*x = b or A'*x = b */ void luf_delete(LUF *luf); /* delete LU-factorization */ #endif /* eof */