diff options
Diffstat (limited to 'test/monniaux/glpk-4.65/src/bflib/scf.h')
| -rw-r--r-- | test/monniaux/glpk-4.65/src/bflib/scf.h | 211 |
1 files changed, 211 insertions, 0 deletions
diff --git a/test/monniaux/glpk-4.65/src/bflib/scf.h b/test/monniaux/glpk-4.65/src/bflib/scf.h new file mode 100644 index 00000000..69d8cfc2 --- /dev/null +++ b/test/monniaux/glpk-4.65/src/bflib/scf.h @@ -0,0 +1,211 @@ +/* scf.h (sparse updatable Schur-complement-based factorization) */ + +/*********************************************************************** +* This code is part of GLPK (GNU Linear Programming Kit). +* +* Copyright (C) 2013-2014 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 SCF_H +#define SCF_H + +#include "btf.h" +#include "ifu.h" +#include "luf.h" + +/*********************************************************************** +* The structure SCF describes sparse updatable factorization based on +* Schur complement. +* +* The SCF-factorization has the following format: +* +* ( A A1~ ) ( A0 A1 ) ( R0 ) ( S0 S ) +* ( ) = P ( ) Q = P ( ) ( ) Q, (1) +* ( A2~ A3~ ) ( A2 A3 ) ( R I ) ( C ) +* +* where: +* +* A is current (unsymmetric) square matrix (not stored); +* +* A1~, A2~, A3~ are some additional matrices (not stored); +* +* A0 is initial (unsymmetric) square matrix (not stored); +* +* A1, A2, A3 are some additional matrices (not stored); +* +* R0 and S0 are matrices that define factorization of the initial +* matrix A0 = R0 * S0 (stored in an invertable form); +* +* R is a matrix defined from R * S0 = A2, so R = A2 * inv(S0) (stored +* in row-wise sparse format); +* +* S is a matrix defined from R0 * S = A1, so S = inv(R0) * A1 (stored +* in column-wise sparse format); +* +* C is Schur complement (to matrix A0) defined from R * S + C = A3, +* so C = A3 - R * S = A3 - A2 * inv(A0) * A1 (stored in an invertable +* form). +* +* P, Q are permutation matrices (stored in both row- and column-like +* formats). */ + +typedef struct SCF SCF; + +struct SCF +{ /* Schur-complement-based factorization */ + int n; + /* order of current matrix A */ + /*--------------------------------------------------------------*/ + /* initial matrix A0 = R0 * S0 of order n0 in invertable form */ + int n0; + /* order of matrix A0 */ + int type; + /* type of factorization used: + * 1 - LU-factorization (R0 = F0, S0 = V0) + * 2 - BT-factorization (R0 = I, S0 = A0) */ + union + { LUF *luf; /* type = 1 */ + BTF *btf; /* type = 2 */ + } a0; + /* factorization of matrix A0 */ + /*--------------------------------------------------------------*/ + /* augmented matrix (A0, A1; A2, A3) of order n0+nn */ + int nn_max; + /* maximal number of additional rows and columns in the augmented + * matrix (this limits the number of updates) */ + int nn; + /* current number of additional rows and columns in the augmented + * matrix, 0 <= nn <= nn_max */ + SVA *sva; + /* associated sparse vector area (SVA) used to store rows of + * matrix R and columns of matrix S */ + /*--------------------------------------------------------------*/ + /* nn*n0-matrix R in row-wise format */ + int rr_ref; + /* reference number of sparse vector in SVA, which is the first + * row of matrix R */ +#if 0 + 0 + int *rr_ptr = &sva->ptr[rr_ref-1]; + /* rr_ptr[0] is not used; + * rr_ptr[i], 1 <= i <= nn, is pointer to i-th row in SVA; + * rr_ptr[nn+1,...,nn_max] are reserved locations */ + int *rr_len = &sva->len[rr_ref-1]; + /* rr_len[0] is not used; + * rr_len[i], 1 <= i <= nn, is length of i-th row; + * rr_len[nn+1,...,nn_max] are reserved locations */ +#endif + /*--------------------------------------------------------------*/ + /* n0*nn-matrix S in column-wise format */ + int ss_ref; + /* reference number of sparse vector in SVA, which is the first + * column of matrix S */ +#if 0 + 0 + int *ss_ptr = &sva->ptr[ss_ref-1]; + /* ss_ptr[0] is not used; + * ss_ptr[j], 1 <= j <= nn, is pointer to j-th column in SVA; + * ss_ptr[nn+1,...,nn_max] are reserved locations */ + int *ss_len = &sva->len[ss_ref-1]; + /* ss_len[0] is not used; + * ss_len[j], 1 <= j <= nn, is length of j-th column; + * ss_len[nn+1,...,nn_max] are reserved locations */ +#endif + /*--------------------------------------------------------------*/ + /* Schur complement C of order nn in invertable form */ + IFU ifu; + /* IFU-factorization of matrix C */ + /*--------------------------------------------------------------*/ + /* permutation matrix P of order n0+nn */ + int *pp_ind; /* int pp_ind[1+n0+nn_max]; */ + /* pp_ind[i] = j means that P[i,j] = 1 */ + int *pp_inv; /* int pp_inv[1+n0+nn_max]; */ + /* pp_inv[j] = i means that P[i,j] = 1 */ + /*--------------------------------------------------------------*/ + /* permutation matrix Q of order n0+nn */ + int *qq_ind; /* int qq_ind[1+n0+nn_max]; */ + /* qq_ind[i] = j means that Q[i,j] = 1 */ + int *qq_inv; /* int qq_inv[1+n0+nn_max]; */ + /* qq_inv[j] = i means that Q[i,j] = 1 */ +}; + +#define scf_swap_q_cols(j1, j2) \ + do \ + { int i1, i2; \ + i1 = qq_inv[j1], i2 = qq_inv[j2]; \ + qq_ind[i1] = j2, qq_inv[j2] = i1; \ + qq_ind[i2] = j1, qq_inv[j1] = i2; \ + } while (0) +/* swap columns j1 and j2 of permutation matrix Q */ + +#define scf_r0_solve _glp_scf_r0_solve +void scf_r0_solve(SCF *scf, int tr, double x[/*1+n0*/]); +/* solve system R0 * x = b or R0'* x = b */ + +#define scf_s0_solve _glp_scf_s0_solve +void scf_s0_solve(SCF *scf, int tr, double x[/*1+n0*/], + double w1[/*1+n0*/], double w2[/*1+n0*/], double w3[/*1+n0*/]); +/* solve system S0 * x = b or S0'* x = b */ + +#define scf_r_prod _glp_scf_r_prod +void scf_r_prod(SCF *scf, double y[/*1+nn*/], double a, const double + x[/*1+n0*/]); +/* compute product y := y + alpha * R * x */ + +#define scf_rt_prod _glp_scf_rt_prod +void scf_rt_prod(SCF *scf, double y[/*1+n0*/], double a, const double + x[/*1+nn*/]); +/* compute product y := y + alpha * R'* x */ + +#define scf_s_prod _glp_scf_s_prod +void scf_s_prod(SCF *scf, double y[/*1+n0*/], double a, const double + x[/*1+nn*/]); +/* compute product y := y + alpha * S * x */ + +#define scf_st_prod _glp_scf_st_prod +void scf_st_prod(SCF *scf, double y[/*1+nn*/], double a, const double + x[/*1+n0*/]); +/* compute product y := y + alpha * S'* x */ + +#define scf_a_solve _glp_scf_a_solve +void scf_a_solve(SCF *scf, double x[/*1+n*/], + double w[/*1+n0+nn*/], double work1[/*1+max(n0,nn)*/], + double work2[/*1+n*/], double work3[/*1+n*/]); +/* solve system A * x = b */ + +#define scf_at_solve _glp_scf_at_solve +void scf_at_solve(SCF *scf, double x[/*1+n*/], + double w[/*1+n0+nn*/], double work1[/*1+max(n0,nn)*/], + double work2[/*1+n*/], double work3[/*1+n*/]); +/* solve system A'* x = b */ + +#define scf_add_r_row _glp_scf_add_r_row +void scf_add_r_row(SCF *scf, const double w[/*1+n0*/]); +/* add new row to matrix R */ + +#define scf_add_s_col _glp_scf_add_s_col +void scf_add_s_col(SCF *scf, const double v[/*1+n0*/]); +/* add new column to matrix S */ + +#define scf_update_aug _glp_scf_update_aug +int scf_update_aug(SCF *scf, double b[/*1+n0*/], double d[/*1+n0*/], + double f[/*1+nn*/], double g[/*1+nn*/], double h, int upd, + double w1[/*1+n0*/], double w2[/*1+n0*/], double w3[/*1+n0*/]); +/* update factorization of augmented matrix */ + +#endif + +/* eof */ |