[9] | 1 | /* glplpf.h (LP basis factorization, Schur complement version) */ |
---|
| 2 | |
---|
| 3 | /*********************************************************************** |
---|
| 4 | * This code is part of GLPK (GNU Linear Programming Kit). |
---|
| 5 | * |
---|
| 6 | * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, |
---|
| 7 | * 2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics, |
---|
| 8 | * Moscow Aviation Institute, Moscow, Russia. All rights reserved. |
---|
| 9 | * E-mail: <mao@gnu.org>. |
---|
| 10 | * |
---|
| 11 | * GLPK is free software: you can redistribute it and/or modify it |
---|
| 12 | * under the terms of the GNU General Public License as published by |
---|
| 13 | * the Free Software Foundation, either version 3 of the License, or |
---|
| 14 | * (at your option) any later version. |
---|
| 15 | * |
---|
| 16 | * GLPK is distributed in the hope that it will be useful, but WITHOUT |
---|
| 17 | * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY |
---|
| 18 | * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public |
---|
| 19 | * License for more details. |
---|
| 20 | * |
---|
| 21 | * You should have received a copy of the GNU General Public License |
---|
| 22 | * along with GLPK. If not, see <http://www.gnu.org/licenses/>. |
---|
| 23 | ***********************************************************************/ |
---|
| 24 | |
---|
| 25 | #ifndef GLPLPF_H |
---|
| 26 | #define GLPLPF_H |
---|
| 27 | |
---|
| 28 | #include "glpscf.h" |
---|
| 29 | #include "glpluf.h" |
---|
| 30 | |
---|
| 31 | /*********************************************************************** |
---|
| 32 | * The structure LPF defines the factorization of the basis mxm matrix |
---|
| 33 | * B, where m is the number of rows in corresponding problem instance. |
---|
| 34 | * |
---|
| 35 | * This factorization is the following septet: |
---|
| 36 | * |
---|
| 37 | * [B] = (L0, U0, R, S, C, P, Q), (1) |
---|
| 38 | * |
---|
| 39 | * and is based on the following main equality: |
---|
| 40 | * |
---|
| 41 | * ( B F^) ( B0 F ) ( L0 0 ) ( U0 R ) |
---|
| 42 | * ( ) = P ( ) Q = P ( ) ( ) Q, (2) |
---|
| 43 | * ( G^ H^) ( G H ) ( S I ) ( 0 C ) |
---|
| 44 | * |
---|
| 45 | * where: |
---|
| 46 | * |
---|
| 47 | * B is the current basis matrix (not stored); |
---|
| 48 | * |
---|
| 49 | * F^, G^, H^ are some additional matrices (not stored); |
---|
| 50 | * |
---|
| 51 | * B0 is some initial basis matrix (not stored); |
---|
| 52 | * |
---|
| 53 | * F, G, H are some additional matrices (not stored); |
---|
| 54 | * |
---|
| 55 | * P, Q are permutation matrices (stored in both row- and column-like |
---|
| 56 | * formats); |
---|
| 57 | * |
---|
| 58 | * L0, U0 are some matrices that defines a factorization of the initial |
---|
| 59 | * basis matrix B0 = L0 * U0 (stored in an invertable form); |
---|
| 60 | * |
---|
| 61 | * R is a matrix defined from L0 * R = F, so R = inv(L0) * F (stored in |
---|
| 62 | * a column-wise sparse format); |
---|
| 63 | * |
---|
| 64 | * S is a matrix defined from S * U0 = G, so S = G * inv(U0) (stored in |
---|
| 65 | * a row-wise sparse format); |
---|
| 66 | * |
---|
| 67 | * C is the Schur complement for matrix (B0 F G H). It is defined from |
---|
| 68 | * S * R + C = H, so C = H - S * R = H - G * inv(U0) * inv(L0) * F = |
---|
| 69 | * = H - G * inv(B0) * F. Matrix C is stored in an invertable form. |
---|
| 70 | * |
---|
| 71 | * REFERENCES |
---|
| 72 | * |
---|
| 73 | * 1. M.A.Saunders, "LUSOL: A basis package for constrained optimiza- |
---|
| 74 | * tion," SCCM, Stanford University, 2006. |
---|
| 75 | * |
---|
| 76 | * 2. M.A.Saunders, "Notes 5: Basis Updates," CME 318, Stanford Univer- |
---|
| 77 | * sity, Spring 2006. |
---|
| 78 | * |
---|
| 79 | * 3. M.A.Saunders, "Notes 6: LUSOL---a Basis Factorization Package," |
---|
| 80 | * ibid. */ |
---|
| 81 | |
---|
| 82 | typedef struct LPF LPF; |
---|
| 83 | |
---|
| 84 | struct LPF |
---|
| 85 | { /* LP basis factorization */ |
---|
| 86 | int valid; |
---|
| 87 | /* the factorization is valid only if this flag is set */ |
---|
| 88 | /*--------------------------------------------------------------*/ |
---|
| 89 | /* initial basis matrix B0 */ |
---|
| 90 | int m0_max; |
---|
| 91 | /* maximal value of m0 (increased automatically, if necessary) */ |
---|
| 92 | int m0; |
---|
| 93 | /* the order of B0 */ |
---|
| 94 | LUF *luf; |
---|
| 95 | /* LU-factorization of B0 */ |
---|
| 96 | /*--------------------------------------------------------------*/ |
---|
| 97 | /* current basis matrix B */ |
---|
| 98 | int m; |
---|
| 99 | /* the order of B */ |
---|
| 100 | double *B; /* double B[1+m*m]; */ |
---|
| 101 | /* B in dense format stored by rows and used only for debugging; |
---|
| 102 | normally this array is not allocated */ |
---|
| 103 | /*--------------------------------------------------------------*/ |
---|
| 104 | /* augmented matrix (B0 F G H) of the order m0+n */ |
---|
| 105 | int n_max; |
---|
| 106 | /* maximal number of additional rows and columns */ |
---|
| 107 | int n; |
---|
| 108 | /* current number of additional rows and columns */ |
---|
| 109 | /*--------------------------------------------------------------*/ |
---|
| 110 | /* m0xn matrix R in column-wise format */ |
---|
| 111 | int *R_ptr; /* int R_ptr[1+n_max]; */ |
---|
| 112 | /* R_ptr[j], 1 <= j <= n, is a pointer to j-th column */ |
---|
| 113 | int *R_len; /* int R_len[1+n_max]; */ |
---|
| 114 | /* R_len[j], 1 <= j <= n, is the length of j-th column */ |
---|
| 115 | /*--------------------------------------------------------------*/ |
---|
| 116 | /* nxm0 matrix S in row-wise format */ |
---|
| 117 | int *S_ptr; /* int S_ptr[1+n_max]; */ |
---|
| 118 | /* S_ptr[i], 1 <= i <= n, is a pointer to i-th row */ |
---|
| 119 | int *S_len; /* int S_len[1+n_max]; */ |
---|
| 120 | /* S_len[i], 1 <= i <= n, is the length of i-th row */ |
---|
| 121 | /*--------------------------------------------------------------*/ |
---|
| 122 | /* Schur complement C of the order n */ |
---|
| 123 | SCF *scf; /* SCF scf[1:n_max]; */ |
---|
| 124 | /* factorization of the Schur complement */ |
---|
| 125 | /*--------------------------------------------------------------*/ |
---|
| 126 | /* matrix P of the order m0+n */ |
---|
| 127 | int *P_row; /* int P_row[1+m0_max+n_max]; */ |
---|
| 128 | /* P_row[i] = j means that P[i,j] = 1 */ |
---|
| 129 | int *P_col; /* int P_col[1+m0_max+n_max]; */ |
---|
| 130 | /* P_col[j] = i means that P[i,j] = 1 */ |
---|
| 131 | /*--------------------------------------------------------------*/ |
---|
| 132 | /* matrix Q of the order m0+n */ |
---|
| 133 | int *Q_row; /* int Q_row[1+m0_max+n_max]; */ |
---|
| 134 | /* Q_row[i] = j means that Q[i,j] = 1 */ |
---|
| 135 | int *Q_col; /* int Q_col[1+m0_max+n_max]; */ |
---|
| 136 | /* Q_col[j] = i means that Q[i,j] = 1 */ |
---|
| 137 | /*--------------------------------------------------------------*/ |
---|
| 138 | /* Sparse Vector Area (SVA) is a set of locations intended to |
---|
| 139 | store sparse vectors which represent columns of matrix R and |
---|
| 140 | rows of matrix S; each location is a doublet (ind, val), where |
---|
| 141 | ind is an index, val is a numerical value of a sparse vector |
---|
| 142 | element; in the whole each sparse vector is a set of adjacent |
---|
| 143 | locations defined by a pointer to its first element and its |
---|
| 144 | length, i.e. the number of its elements */ |
---|
| 145 | int v_size; |
---|
| 146 | /* the SVA size, in locations; locations are numbered by integers |
---|
| 147 | 1, 2, ..., v_size, and location 0 is not used */ |
---|
| 148 | int v_ptr; |
---|
| 149 | /* pointer to the first available location */ |
---|
| 150 | int *v_ind; /* int v_ind[1+v_size]; */ |
---|
| 151 | /* v_ind[k], 1 <= k <= v_size, is the index field of location k */ |
---|
| 152 | double *v_val; /* double v_val[1+v_size]; */ |
---|
| 153 | /* v_val[k], 1 <= k <= v_size, is the value field of location k */ |
---|
| 154 | /*--------------------------------------------------------------*/ |
---|
| 155 | double *work1; /* double work1[1+m0+n_max]; */ |
---|
| 156 | /* working array */ |
---|
| 157 | double *work2; /* double work2[1+m0+n_max]; */ |
---|
| 158 | /* working array */ |
---|
| 159 | }; |
---|
| 160 | |
---|
| 161 | /* return codes: */ |
---|
| 162 | #define LPF_ESING 1 /* singular matrix */ |
---|
| 163 | #define LPF_ECOND 2 /* ill-conditioned matrix */ |
---|
| 164 | #define LPF_ELIMIT 3 /* update limit reached */ |
---|
| 165 | |
---|
| 166 | #define lpf_create_it _glp_lpf_create_it |
---|
| 167 | LPF *lpf_create_it(void); |
---|
| 168 | /* create LP basis factorization */ |
---|
| 169 | |
---|
| 170 | #define lpf_factorize _glp_lpf_factorize |
---|
| 171 | int lpf_factorize(LPF *lpf, int m, const int bh[], int (*col) |
---|
| 172 | (void *info, int j, int ind[], double val[]), void *info); |
---|
| 173 | /* compute LP basis factorization */ |
---|
| 174 | |
---|
| 175 | #define lpf_ftran _glp_lpf_ftran |
---|
| 176 | void lpf_ftran(LPF *lpf, double x[]); |
---|
| 177 | /* perform forward transformation (solve system B*x = b) */ |
---|
| 178 | |
---|
| 179 | #define lpf_btran _glp_lpf_btran |
---|
| 180 | void lpf_btran(LPF *lpf, double x[]); |
---|
| 181 | /* perform backward transformation (solve system B'*x = b) */ |
---|
| 182 | |
---|
| 183 | #define lpf_update_it _glp_lpf_update_it |
---|
| 184 | int lpf_update_it(LPF *lpf, int j, int bh, int len, const int ind[], |
---|
| 185 | const double val[]); |
---|
| 186 | /* update LP basis factorization */ |
---|
| 187 | |
---|
| 188 | #define lpf_delete_it _glp_lpf_delete_it |
---|
| 189 | void lpf_delete_it(LPF *lpf); |
---|
| 190 | /* delete LP basis factorization */ |
---|
| 191 | |
---|
| 192 | #endif |
---|
| 193 | |
---|
| 194 | /* eof */ |
---|