You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
203 lines
8.2 KiB
203 lines
8.2 KiB
/* sgf.h (sparse Gaussian factorizer) */
|
|
|
|
/***********************************************************************
|
|
* This code is part of GLPK (GNU Linear Programming Kit).
|
|
*
|
|
* Copyright (C) 2012, 2013 Andrew Makhorin, Department for Applied
|
|
* Informatics, Moscow Aviation Institute, Moscow, Russia. All rights
|
|
* reserved. E-mail: <mao@gnu.org>.
|
|
*
|
|
* 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 3 of the License, 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. If not, see <http://www.gnu.org/licenses/>.
|
|
***********************************************************************/
|
|
|
|
#ifndef SGF_H
|
|
#define SGF_H
|
|
|
|
#include "luf.h"
|
|
|
|
typedef struct SGF SGF;
|
|
|
|
struct SGF
|
|
{ /* sparse Gaussian factorizer workspace */
|
|
LUF *luf;
|
|
/* LU-factorization being computed */
|
|
/*--------------------------------------------------------------*/
|
|
/* to efficiently choose pivot elements according to Markowitz
|
|
* strategy, the search technique proposed by Iain Duff is used;
|
|
* it is based on using two families of sets {R[0], ..., R[n]}
|
|
* and {C[0], ..., C[n]}, where R[k] and C[k], 0 <= k <= n, are,
|
|
* respectively, sets of rows and columns of the active submatrix
|
|
* of matrix V having k non-zeros (i.e. whose length is k); each
|
|
* set R[k] and C[k] is implemented as a doubly linked list */
|
|
int *rs_head; /* int rs_head[1+n]; */
|
|
/* rs_head[k], 0 <= k <= n, is the number of first row, which
|
|
* has k non-zeros in the active submatrix */
|
|
int *rs_prev; /* int rs_prev[1+n]; */
|
|
/* rs_prev[0] is not used;
|
|
* rs_prev[i], 1 <= i <= n, is the number of previous row, which
|
|
* has the same number of non-zeros as i-th row;
|
|
* rs_prev[i] < 0 means that i-th row is inactive */
|
|
int *rs_next; /* int rs_next[1+n]; */
|
|
/* rs_next[0] is not used;
|
|
* rs_next[i], 1 <= i <= n, is the number of next row, which has
|
|
* the same number of non-zeros as i-th row;
|
|
* rs_next[i] < 0 means that i-th row is inactive */
|
|
int *cs_head; /* int cs_head[1+n]; */
|
|
/* cs_head[k], 0 <= k <= n, is the number of first 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 the number of previous column,
|
|
* which has the same number of non-zeros as j-th column;
|
|
* cs_prev[j] < 0 means that j-th column is inactive */
|
|
int *cs_next; /* int cs_next[1+n]; */
|
|
/* cs_next[0] is not used;
|
|
* cs_next[j], 1 <= j <= n, is the number of next column, which
|
|
* has the same number of non-zeros as j-th column;
|
|
* cs_next[j] < 0 means that j-th column is inactive */
|
|
/* NOTE: cs_prev[j] = cs_next[j] = j means that j-th column was
|
|
* temporarily removed from corresponding set C[k] by the
|
|
* pivoting routine according to Uwe Suhl's heuristic */
|
|
/*--------------------------------------------------------------*/
|
|
/* working arrays */
|
|
double *vr_max; /* int vr_max[1+n]; */
|
|
/* vr_max[0] is not used;
|
|
* vr_max[i], 1 <= i <= n, is used only if i-th row of matrix V
|
|
* is active (i.e. belongs to the active submatrix), and is the
|
|
* largest magnitude of elements in that row; if vr_max[i] < 0,
|
|
* the largest magnitude is unknown yet */
|
|
char *flag; /* char flag[1+n]; */
|
|
/* boolean working array */
|
|
double *work; /* double work[1+n]; */
|
|
/* floating-point working array */
|
|
/*--------------------------------------------------------------*/
|
|
/* control parameters */
|
|
int updat;
|
|
/* if this flag is set, the matrix V is assumed to be updatable;
|
|
* in this case factorized (non-active) part of V is stored in
|
|
* the left part of SVA rather than in its right part */
|
|
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 criterion |v[i,j]| >= piv_tol * max |v[i,*]|,
|
|
* i.e. if it is not very small in the magnitude among other
|
|
* elements in the same row; decreasing this parameter gives
|
|
* 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
|
|
* proposed by Uwe Suhl: if a column of the active submatrix has
|
|
* no eligible pivot candidates (i.e. all its elements do not
|
|
* satisfy to the stability criterion), the routine excludes it
|
|
* from futher consideration until it becomes column singleton;
|
|
* in many cases this allows reducing the time needed to choose
|
|
* the pivot */
|
|
double eps_tol;
|
|
/* epsilon tolerance; each element of the active submatrix, whose
|
|
* magnitude is less than eps_tol, is replaced by exact zero */
|
|
#if 0 /* FIXME */
|
|
double den_lim;
|
|
/* density limit; if the density of the active submatrix reaches
|
|
* this limit, the factorization routine switches from sparse to
|
|
* dense mode */
|
|
#endif
|
|
};
|
|
|
|
#define sgf_activate_row(i) \
|
|
do \
|
|
{ int len = vr_len[i]; \
|
|
rs_prev[i] = 0; \
|
|
rs_next[i] = rs_head[len]; \
|
|
if (rs_next[i] != 0) \
|
|
rs_prev[rs_next[i]] = i; \
|
|
rs_head[len] = i; \
|
|
} while (0)
|
|
/* include i-th row of matrix V in active set R[len] */
|
|
|
|
#define sgf_deactivate_row(i) \
|
|
do \
|
|
{ if (rs_prev[i] == 0) \
|
|
rs_head[vr_len[i]] = rs_next[i]; \
|
|
else \
|
|
rs_next[rs_prev[i]] = rs_next[i]; \
|
|
if (rs_next[i] == 0) \
|
|
; \
|
|
else \
|
|
rs_prev[rs_next[i]] = rs_prev[i]; \
|
|
rs_prev[i] = rs_next[i] = -1; \
|
|
} while (0)
|
|
/* remove i-th row of matrix V from active set R[len] */
|
|
|
|
#define sgf_activate_col(j) \
|
|
do \
|
|
{ int len = vc_len[j]; \
|
|
cs_prev[j] = 0; \
|
|
cs_next[j] = cs_head[len]; \
|
|
if (cs_next[j] != 0) \
|
|
cs_prev[cs_next[j]] = j; \
|
|
cs_head[len] = j; \
|
|
} while (0)
|
|
/* include j-th column of matrix V in active set C[len] */
|
|
|
|
#define sgf_deactivate_col(j) \
|
|
do \
|
|
{ if (cs_prev[j] == 0) \
|
|
cs_head[vc_len[j]] = cs_next[j]; \
|
|
else \
|
|
cs_next[cs_prev[j]] = cs_next[j]; \
|
|
if (cs_next[j] == 0) \
|
|
; \
|
|
else \
|
|
cs_prev[cs_next[j]] = cs_prev[j]; \
|
|
cs_prev[j] = cs_next[j] = -1; \
|
|
} while (0)
|
|
/* remove j-th column of matrix V from active set C[len] */
|
|
|
|
#define sgf_reduce_nuc _glp_sgf_reduce_nuc
|
|
void sgf_reduce_nuc(LUF *luf, int *k1, int *k2, int cnt[/*1+n*/],
|
|
int list[/*1+n*/]);
|
|
/* initial reordering to minimize nucleus size */
|
|
|
|
#define sgf_singl_phase _glp_sgf_singl_phase
|
|
int sgf_singl_phase(LUF *luf, int k1, int k2, int updat,
|
|
int ind[/*1+n*/], double val[/*1+n*/]);
|
|
/* compute LU-factorization (singleton phase) */
|
|
|
|
#define sgf_choose_pivot _glp_sgf_choose_pivot
|
|
int sgf_choose_pivot(SGF *sgf, int *p, int *q);
|
|
/* choose pivot element v[p,q] */
|
|
|
|
#define sgf_eliminate _glp_sgf_eliminate
|
|
int sgf_eliminate(SGF *sgf, int p, int q);
|
|
/* perform gaussian elimination */
|
|
|
|
#define sgf_dense_lu _glp_sgf_dense_lu
|
|
int sgf_dense_lu(int n, double a[], int r[], int c[], double eps);
|
|
/* compute dense LU-factorization with full pivoting */
|
|
|
|
#define sgf_dense_phase _glp_sgf_dense_phase
|
|
int sgf_dense_phase(LUF *luf, int k, int updat);
|
|
/* compute LU-factorization (dense phase) */
|
|
|
|
#define sgf_factorize _glp_sgf_factorize
|
|
int sgf_factorize(SGF *sgf, int singl);
|
|
/* compute LU-factorization (main routine) */
|
|
|
|
#endif
|
|
|
|
/* eof */
|