Loading...
Searching...
No Matches
glpluf.h
Go to the documentation of this file.
/* glpluf.h */
/*----------------------------------------------------------------------
-- Copyright (C) 2000, 2001, 2002 Andrew Makhorin <mao@mai2.rcnet.ru>,
-- Department for Applied Informatics, Moscow Aviation
-- Institute, Moscow, Russia. All rights reserved.
--
-- This file is a part of GLPK (GNU Linear Programming Kit).
--
-- Licensed under the Eclipse Public License (EPL) by permission of the
-- author for inclusion in the DyLP LP distribution.
----------------------------------------------------------------------*/
/*
@(#)glpluf.h 1.1 10/18/02
svn/cvs: $Id: glpluf.h 407 2010-12-31 20:48:48Z lou $
*/
#ifndef _GLPLUF_H
#define _GLPLUF_H
#define luf_create dy_glp_luf_create
#define luf_defrag_sva dy_glp_luf_defrag_sva
#define luf_enlarge_row dy_glp_luf_enlarge_row
#define luf_enlarge_col dy_glp_luf_enlarge_col
#define luf_alloc_wa dy_glp_luf_alloc_wa
#define luf_free_wa dy_glp_luf_free_wa
#define luf_decomp dy_glp_luf_decomp
#define luf_f_solve dy_glp_luf_f_solve
#define luf_v_solve dy_glp_luf_v_solve
#define luf_solve dy_glp_luf_solve
#define luf_delete dy_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] */
double *vr_piv; /* double 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 doublet (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 */
double *sv_val; /* double 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 double 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 */
double *work; /* double 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 */
double 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 */
double eps_tol;
/* epsilon tolerance; each element of the matrix V with absolute
value less than eps_tol is replaced by exact zero */
double 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) */
double max_a;
/* largest of absolute values of elements of the matrix A */
double 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) */
double *rs_max; /* double 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[], double aj[]),
LUF_WA *wa);
/* compute LU-factorization */
void luf_f_solve(LUF *luf, int tr, double x[]);
/* solve system F*x = b or F'*x = b */
void luf_v_solve(LUF *luf, int tr, double x[]);
/* solve system V*x = b or V'*x = b */
void luf_solve(LUF *luf, int tr, double x[]);
/* solve system A*x = b or A'*x = b */
void luf_delete(LUF *luf);
/* delete LU-factorization */
#endif
/* eof */
luf_f_solve
#define luf_f_solve
Definition glpluf.h:29
luf_defrag_sva
#define luf_defrag_sva
Definition glpluf.h:23
luf_v_solve
#define luf_v_solve
Definition glpluf.h:30
luf_create
#define luf_create
Definition glpluf.h:22
luf_delete
#define luf_delete
Definition glpluf.h:32
luf_solve
#define luf_solve
Definition glpluf.h:31
luf_enlarge_row
#define luf_enlarge_row
Definition glpluf.h:24
luf_alloc_wa
#define luf_alloc_wa
Definition glpluf.h:26
luf_decomp
#define luf_decomp
Definition glpluf.h:28
luf_free_wa
#define luf_free_wa
Definition glpluf.h:27
luf_enlarge_col
#define luf_enlarge_col
Definition glpluf.h:25
LUF_WA
Definition glpluf.h:271
LUF_WA::cs_head
int * cs_head
Definition glpluf.h:297
LUF_WA::cs_next
int * cs_next
Definition glpluf.h:304
LUF_WA::rs_head
int * rs_head
Definition glpluf.h:286
LUF_WA::rs_prev
int * rs_prev
Definition glpluf.h:289
LUF_WA::cs_prev
int * cs_prev
Definition glpluf.h:300
LUF_WA::rs_next
int * rs_next
Definition glpluf.h:293
LUF_WA::rs_max
double * rs_max
Definition glpluf.h:272
LUF
Definition glpluf.h:84
LUF::sv_head
int sv_head
Definition glpluf.h:200
LUF::big_v
double big_v
Definition glpluf.h:263
LUF::sv_prev
int * sv_prev
Definition glpluf.h:204
LUF::rank
int rank
Definition glpluf.h:266
LUF::vc_ptr
int * vc_ptr
Definition glpluf.h:130
LUF::piv_tol
double piv_tol
Definition glpluf.h:221
LUF::fr_ptr
int * fr_ptr
Definition glpluf.h:91
LUF::qq_row
int * qq_row
Definition glpluf.h:155
LUF::new_sva
int new_sva
Definition glpluf.h:218
LUF::pp_col
int * pp_col
Definition glpluf.h:147
LUF::fc_len
int * fc_len
Definition glpluf.h:105
LUF::eps_tol
double eps_tol
Definition glpluf.h:241
LUF::sv_beg
int sv_beg
Definition glpluf.h:179
LUF::nnz_f
int nnz_f
Definition glpluf.h:254
LUF::suhl
int suhl
Definition glpluf.h:233
LUF::work
double * work
Definition glpluf.h:214
LUF::vr_cap
int * vr_cap
Definition glpluf.h:119
LUF::max_gro
double max_gro
Definition glpluf.h:244
LUF::vc_len
int * vc_len
Definition glpluf.h:134
LUF::piv_lim
int piv_lim
Definition glpluf.h:229
LUF::n
int n
Definition glpluf.h:85
LUF::sv_val
double * sv_val
Definition glpluf.h:190
LUF::pp_row
int * pp_row
Definition glpluf.h:145
LUF::nnz_v
int nnz_v
Definition glpluf.h:257
LUF::vr_len
int * vr_len
Definition glpluf.h:115
LUF::vc_cap
int * vc_cap
Definition glpluf.h:138
LUF::fr_len
int * fr_len
Definition glpluf.h:95
LUF::qq_col
int * qq_col
Definition glpluf.h:157
LUF::nnz_a
int nnz_a
Definition glpluf.h:252
LUF::flag
int * flag
Definition glpluf.h:212
LUF::sv_next
int * sv_next
Definition glpluf.h:207
LUF::vr_piv
double * vr_piv
Definition glpluf.h:124
LUF::sv_size
int sv_size
Definition glpluf.h:174
LUF::sv_end
int sv_end
Definition glpluf.h:179
LUF::sv_tail
int sv_tail
Definition glpluf.h:202
LUF::valid
int valid
Definition glpluf.h:87
LUF::max_a
double max_a
Definition glpluf.h:261
LUF::vr_ptr
int * vr_ptr
Definition glpluf.h:111
LUF::fc_ptr
int * fc_ptr
Definition glpluf.h:101
LUF::sv_ndx
int * sv_ndx
Definition glpluf.h:186