PageRenderTime 18ms CodeModel.GetById 8ms app.highlight 6ms RepoModel.GetById 1ms app.codeStats 0ms

/farmR/src/glpluf.h

https://code.google.com/p/javawfm/
C++ Header | 325 lines | 75 code | 17 blank | 233 comment | 0 complexity | 2234ded9d81391816113fbd647bfacac MD5 | raw file
  1/* glpluf.h (LU-factorization) */
  2
  3/***********************************************************************
  4*  This code is part of GLPK (GNU Linear Programming Kit).
  5*
  6*  Copyright (C) 2000,01,02,03,04,05,06,07,08,2009 Andrew Makhorin,
  7*  Department for Applied Informatics, Moscow Aviation Institute,
  8*  Moscow, Russia. All rights reserved. E-mail: <mao@mai2.rcnet.ru>.
  9*
 10*  GLPK is free software: you can redistribute it and/or modify it
 11*  under the terms of the GNU General Public License as published by
 12*  the Free Software Foundation, either version 3 of the License, or
 13*  (at your option) any later version.
 14*
 15*  GLPK is distributed in the hope that it will be useful, but WITHOUT
 16*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
 17*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
 18*  License for more details.
 19*
 20*  You should have received a copy of the GNU General Public License
 21*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
 22***********************************************************************/
 23
 24#ifndef _GLPLUF_H
 25#define _GLPLUF_H
 26
 27/***********************************************************************
 28*  The structure LUF defines LU-factorization of a square matrix A and
 29*  is the following quartet:
 30*
 31*     [A] = (F, V, P, Q),                                            (1)
 32*
 33*  where F and V are such matrices that
 34*
 35*     A = F * V,                                                     (2)
 36*
 37*  and P and Q are such permutation matrices that the matrix
 38*
 39*     L = P * F * inv(P)                                             (3)
 40*
 41*  is lower triangular with unity diagonal, and the matrix
 42*
 43*     U = P * V * Q                                                  (4)
 44*
 45*  is upper triangular. All the matrices have the order n.
 46*
 47*  Matrices F and V are stored in row- and column-wise sparse format
 48*  as row and column linked lists of non-zero elements. Unity elements
 49*  on the main diagonal of matrix F are not stored. Pivot elements of
 50*  matrix V (which correspond to diagonal elements of matrix U) are
 51*  stored separately in an ordinary array.
 52*
 53*  Permutation matrices P and Q are stored in ordinary arrays in both
 54*  row- and column-like formats.
 55*
 56*  Matrices L and U are completely defined by matrices F, V, P, and Q
 57*  and therefore not stored explicitly.
 58*
 59*  The factorization (1)-(4) is a version of LU-factorization. Indeed,
 60*  from (3) and (4) it follows that:
 61*
 62*     F = inv(P) * L * P,
 63*
 64*     U = inv(P) * U * inv(Q),
 65*
 66*  and substitution into (2) leads to:
 67*
 68*     A = F * V = inv(P) * L * U * inv(Q).
 69*
 70*  For more details see the program documentation. */
 71
 72typedef struct LUF LUF;
 73
 74struct LUF
 75{     /* LU-factorization of a square matrix */
 76      int n_max;
 77      /* maximal value of n (increased automatically, if necessary) */
 78      int n;
 79      /* the order of matrices A, F, V, P, Q */
 80      int valid;
 81      /* the factorization is valid only if this flag is set */
 82      /*--------------------------------------------------------------*/
 83      /* matrix F in row-wise format */
 84      int *fr_ptr; /* int fr_ptr[1+n_max]; */
 85      /* fr_ptr[i], i = 1,...,n, is a pointer to the first element of
 86         i-th row in SVA */
 87      int *fr_len; /* int fr_len[1+n_max]; */
 88      /* fr_len[i], i = 1,...,n, is the number of elements in i-th row
 89         (except unity diagonal element) */
 90      /*--------------------------------------------------------------*/
 91      /* matrix F in column-wise format */
 92      int *fc_ptr; /* int fc_ptr[1+n_max]; */
 93      /* fc_ptr[j], j = 1,...,n, is a pointer to the first element of
 94         j-th column in SVA */
 95      int *fc_len; /* int fc_len[1+n_max]; */
 96      /* fc_len[j], j = 1,...,n, is the number of elements in j-th
 97         column (except unity diagonal element) */
 98      /*--------------------------------------------------------------*/
 99      /* matrix V in row-wise format */
100      int *vr_ptr; /* int vr_ptr[1+n_max]; */
101      /* vr_ptr[i], i = 1,...,n, is a pointer to the first element of
102         i-th row in SVA */
103      int *vr_len; /* int vr_len[1+n_max]; */
104      /* vr_len[i], i = 1,...,n, is the number of elements in i-th row
105         (except pivot element) */
106      int *vr_cap; /* int vr_cap[1+n_max]; */
107      /* vr_cap[i], i = 1,...,n, is the capacity of i-th row, i.e.
108         maximal number of elements which can be stored in the row
109         without relocating it, vr_cap[i] >= vr_len[i] */
110      double *vr_piv; /* double vr_piv[1+n_max]; */
111      /* vr_piv[p], p = 1,...,n, is the pivot element v[p,q] which
112         corresponds to a diagonal element of matrix U = P*V*Q */
113      /*--------------------------------------------------------------*/
114      /* matrix V in column-wise format */
115      int *vc_ptr; /* int vc_ptr[1+n_max]; */
116      /* vc_ptr[j], j = 1,...,n, is a pointer to the first element of
117         j-th column in SVA */
118      int *vc_len; /* int vc_len[1+n_max]; */
119      /* vc_len[j], j = 1,...,n, is the number of elements in j-th
120         column (except pivot element) */
121      int *vc_cap; /* int vc_cap[1+n_max]; */
122      /* vc_cap[j], j = 1,...,n, is the capacity of j-th column, i.e.
123         maximal number of elements which can be stored in the column
124         without relocating it, vc_cap[j] >= vc_len[j] */
125      /*--------------------------------------------------------------*/
126      /* matrix P */
127      int *pp_row; /* int pp_row[1+n_max]; */
128      /* pp_row[i] = j means that P[i,j] = 1 */
129      int *pp_col; /* int pp_col[1+n_max]; */
130      /* pp_col[j] = i means that P[i,j] = 1 */
131      /* if i-th row or column of matrix F is i'-th row or column of
132         matrix L, or if i-th row of matrix V is i'-th row of matrix U,
133         then pp_row[i'] = i and pp_col[i] = i' */
134      /*--------------------------------------------------------------*/
135      /* matrix Q */
136      int *qq_row; /* int qq_row[1+n_max]; */
137      /* qq_row[i] = j means that Q[i,j] = 1 */
138      int *qq_col; /* int qq_col[1+n_max]; */
139      /* qq_col[j] = i means that Q[i,j] = 1 */
140      /* if j-th column of matrix V is j'-th column of matrix U, then
141         qq_row[j] = j' and qq_col[j'] = j */
142      /*--------------------------------------------------------------*/
143      /* the Sparse Vector Area (SVA) is a set of locations used to
144         store sparse vectors representing rows and columns of matrices
145         F and V; each location is a doublet (ind, val), where ind is
146         an index, and val is a numerical value of a sparse vector
147         element; in the whole each sparse vector is a set of adjacent
148         locations defined by a pointer to the first element and the
149         number of elements; these pointer and number are stored in the
150         corresponding matrix data structure (see above); the left part
151         of SVA is used to store rows and columns of matrix V, and its
152         right part is used to store rows and columns of matrix F; the
153         middle part of SVA contains free (unused) locations */
154      int sv_size;
155      /* the size of SVA, in locations; all locations are numbered by
156         integers 1, ..., n, and location 0 is not used; if necessary,
157         the SVA size is automatically increased */
158      int sv_beg, sv_end;
159      /* SVA partitioning pointers:
160         locations from 1 to sv_beg-1 belong to the left part
161         locations from sv_beg to sv_end-1 belong to the middle part
162         locations from sv_end to sv_size belong to the right part
163         the size of the middle part is (sv_end - sv_beg) */
164      int *sv_ind; /* sv_ind[1+sv_size]; */
165      /* sv_ind[k], 1 <= k <= sv_size, is the index field of k-th
166         location */
167      double *sv_val; /* sv_val[1+sv_size]; */
168      /* sv_val[k], 1 <= k <= sv_size, is the value field of k-th
169         location */
170      /*--------------------------------------------------------------*/
171      /* in order to efficiently defragment the left part of SVA there
172         is a doubly linked list of rows and columns of matrix V, where
173         rows are numbered by 1, ..., n, while columns are numbered by
174         n+1, ..., n+n, that allows uniquely identifying each row and
175         column of V by only one integer; in this list rows and columns
176         are ordered by ascending their pointers vr_ptr and vc_ptr */
177      int sv_head;
178      /* the number of leftmost row/column */
179      int sv_tail;
180      /* the number of rightmost row/column */
181      int *sv_prev; /* int sv_prev[1+n_max+n_max]; */
182      /* sv_prev[k], k = 1,...,n+n, is the number of a row/column which
183         precedes k-th row/column */
184      int *sv_next; /* int sv_next[1+n_max+n_max]; */
185      /* sv_next[k], k = 1,...,n+n, is the number of a row/column which
186         succedes k-th row/column */
187      /*--------------------------------------------------------------*/
188      /* working segment (used only during factorization) */
189      double *vr_max; /* int vr_max[1+n_max]; */
190      /* vr_max[i], 1 <= i <= n, is used only if i-th row of matrix V
191         is active (i.e. belongs to the active submatrix), and is the
192         largest magnitude of elements in i-th row; if vr_max[i] < 0,
193         the largest magnitude is not known yet and should be computed
194         by the pivoting routine */
195      /*--------------------------------------------------------------*/
196      /* in order to efficiently implement Markowitz strategy and Duff
197         search technique there are two families {R[0], R[1], ..., R[n]}
198         and {C[0], C[1], ..., C[n]}; member R[k] is the set of active
199         rows of matrix V, which have k non-zeros, and member C[k] is
200         the set of active columns of V, which have k non-zeros in the
201         active submatrix (i.e. in the active rows); each set R[k] and
202         C[k] is implemented as a separate doubly linked list */
203      int *rs_head; /* int rs_head[1+n_max]; */
204      /* rs_head[k], 0 <= k <= n, is the number of first active row,
205         which has k non-zeros */
206      int *rs_prev; /* int rs_prev[1+n_max]; */
207      /* rs_prev[i], 1 <= i <= n, is the number of previous row, which
208         has the same number of non-zeros as i-th row */
209      int *rs_next; /* int rs_next[1+n_max]; */
210      /* rs_next[i], 1 <= i <= n, is the number of next row, which has
211         the same number of non-zeros as i-th row */
212      int *cs_head; /* int cs_head[1+n_max]; */
213      /* cs_head[k], 0 <= k <= n, is the number of first active column,
214         which has k non-zeros (in the active rows) */
215      int *cs_prev; /* int cs_prev[1+n_max]; */
216      /* cs_prev[j], 1 <= j <= n, is the number of previous column,
217         which has the same number of non-zeros (in the active rows) as
218         j-th column */
219      int *cs_next; /* int cs_next[1+n_max]; */
220      /* cs_next[j], 1 <= j <= n, is the number of next column, which
221         has the same number of non-zeros (in the active rows) as j-th
222         column */
223      /* (end of working segment) */
224      /*--------------------------------------------------------------*/
225      /* working arrays */
226      int *flag; /* int flag[1+n_max]; */
227      /* integer working array */
228      double *work; /* double work[1+n_max]; */
229      /* floating-point working array */
230      /*--------------------------------------------------------------*/
231      /* control parameters */
232      int new_sva;
233      /* new required size of the sparse vector area, in locations; set
234         automatically by the factorizing routine */
235      double piv_tol;
236      /* threshold pivoting tolerance, 0 < piv_tol < 1; element v[i,j]
237         of the active submatrix fits to be pivot if it satisfies to the
238         stability criterion |v[i,j]| >= piv_tol * max |v[i,*]|, i.e. if
239         it is not very small in the magnitude among other elements in
240         the same row; decreasing this parameter gives better sparsity
241         at the expense of numerical accuracy and vice versa */
242      int piv_lim;
243      /* maximal allowable number of pivot candidates to be considered;
244         if piv_lim pivot candidates have been considered, the pivoting
245         routine terminates the search with the best candidate found */
246      int suhl;
247      /* if this flag is set, the pivoting routine applies a heuristic
248         proposed by Uwe Suhl: if a column of the active submatrix has
249         no eligible pivot candidates (i.e. all its elements do not
250         satisfy to the stability criterion), the routine excludes it
251         from futher consideration until it becomes column singleton;
252         in many cases this allows reducing the time needed for pivot
253         searching */
254      double eps_tol;
255      /* epsilon tolerance; each element of the active submatrix, whose
256         magnitude is less than eps_tol, is replaced by exact zero */
257      double max_gro;
258      /* maximal allowable growth of elements of matrix V during all
259         the factorization process; if on some eliminaion step the ratio
260         big_v / max_a (see below) becomes greater than max_gro, matrix
261         A is considered as ill-conditioned (assuming that the pivoting
262         tolerance piv_tol has an appropriate value) */
263      /*--------------------------------------------------------------*/
264      /* some statistics */
265      int nnz_a;
266      /* the number of non-zeros in matrix A */
267      int nnz_f;
268      /* the number of non-zeros in matrix F (except diagonal elements,
269         which are not stored) */
270      int nnz_v;
271      /* the number of non-zeros in matrix V (except its pivot elements,
272         which are stored in a separate array) */
273      double max_a;
274      /* the largest magnitude of elements of matrix A */
275      double big_v;
276      /* the largest magnitude of elements of matrix V appeared in the
277         active submatrix during all the factorization process */
278      int rank;
279      /* estimated rank of matrix A */
280};
281
282/* return codes: */
283#define LUF_ESING    1  /* singular matrix */
284#define LUF_ECOND    2  /* ill-conditioned matrix */
285
286#define luf_create_it _glp_luf_create_it
287LUF *luf_create_it(void);
288/* create LU-factorization */
289
290#define luf_defrag_sva _glp_luf_defrag_sva
291void luf_defrag_sva(LUF *luf);
292/* defragment the sparse vector area */
293
294#define luf_enlarge_row _glp_luf_enlarge_row
295int luf_enlarge_row(LUF *luf, int i, int cap);
296/* enlarge row capacity */
297
298#define luf_enlarge_col _glp_luf_enlarge_col
299int luf_enlarge_col(LUF *luf, int j, int cap);
300/* enlarge column capacity */
301
302#define luf_factorize _glp_luf_factorize
303int luf_factorize(LUF *luf, int n, int (*col)(void *info, int j,
304      int ind[], double val[]), void *info);
305/* compute LU-factorization */
306
307#define luf_f_solve _glp_luf_f_solve
308void luf_f_solve(LUF *luf, int tr, double x[]);
309/* solve system F*x = b or F'*x = b */
310
311#define luf_v_solve _glp_luf_v_solve
312void luf_v_solve(LUF *luf, int tr, double x[]);
313/* solve system V*x = b or V'*x = b */
314
315#define luf_a_solve _glp_luf_a_solve
316void luf_a_solve(LUF *luf, int tr, double x[]);
317/* solve system A*x = b or A'*x = b */
318
319#define luf_delete_it _glp_luf_delete_it
320void luf_delete_it(LUF *luf);
321/* delete LU-factorization */
322
323#endif
324
325/* eof */