1 /* glplux.h (LU-factorization, bignum arithmetic) */ 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 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 GLPLUX_H 26 #define GLPLUX_H 27 28 #include "glpdmp.h" 29 #include "glpgmp.h" 30 31 /*---------------------------------------------------------------------- 32 // The structure LUX defines LU-factorization of a square matrix A, 33 // which is the following quartet: 34 // 35 // [A] = (F, V, P, Q), (1) 36 // 37 // where F and V are such matrices that 38 // 39 // A = F * V, (2) 40 // 41 // and P and Q are such permutation matrices that the matrix 42 // 43 // L = P * F * inv(P) (3) 44 // 45 // is lower triangular with unity diagonal, and the matrix 46 // 47 // U = P * V * Q (4) 48 // 49 // is upper triangular. All the matrices have the order n. 50 // 51 // The matrices F and V are stored in row/column-wise sparse format as 52 // row and column linked lists of non-zero elements. Unity elements on 53 // the main diagonal of the matrix F are not stored. Pivot elements of 54 // the matrix V (that correspond to diagonal elements of the matrix U) 55 // are also missing from the row and column lists and stored separately 56 // in an ordinary array. 57 // 58 // The permutation matrices P and Q are stored as ordinary arrays using 59 // both row- and column-like formats. 60 // 61 // The matrices L and U being completely defined by the matrices F, V, 62 // P, and Q are not stored explicitly. 63 // 64 // It is easy to show that the factorization (1)-(3) is some version of 65 // LU-factorization. Indeed, from (3) and (4) it follows that: 66 // 67 // F = inv(P) * L * P, 68 // 69 // V = inv(P) * U * inv(Q), 70 // 71 // and substitution into (2) gives: 72 // 73 // A = F * V = inv(P) * L * U * inv(Q). 74 // 75 // For more details see the program documentation. */ 76 77 typedef struct LUX LUX; 78 typedef struct LUXELM LUXELM; 79 typedef struct LUXWKA LUXWKA; 80 81 struct LUX 82 { /* LU-factorization of a square matrix */ 83 int n; 84 /* the order of matrices A, F, V, P, Q */ 85 DMP *pool; 86 /* memory pool for elements of matrices F and V */ 87 LUXELM **F_row; /* LUXELM *F_row[1+n]; */ 88 /* F_row[0] is not used; 89 F_row[i], 1 <= i <= n, is a pointer to the list of elements in 90 i-th row of matrix F (diagonal elements are not stored) */ 91 LUXELM **F_col; /* LUXELM *F_col[1+n]; */ 92 /* F_col[0] is not used; 93 F_col[j], 1 <= j <= n, is a pointer to the list of elements in 94 j-th column of matrix F (diagonal elements are not stored) */ 95 mpq_t *V_piv; /* mpq_t V_piv[1+n]; */ 96 /* V_piv[0] is not used; 97 V_piv[p], 1 <= p <= n, is a pivot element v[p,q] corresponding 98 to a diagonal element u[k,k] of matrix U = P*V*Q (used on k-th 99 elimination step, k = 1, 2, ..., n) */ 100 LUXELM **V_row; /* LUXELM *V_row[1+n]; */ 101 /* V_row[0] is not used; 102 V_row[i], 1 <= i <= n, is a pointer to the list of elements in 103 i-th row of matrix V (except pivot elements) */ 104 LUXELM **V_col; /* LUXELM *V_col[1+n]; */ 105 /* V_col[0] is not used; 106 V_col[j], 1 <= j <= n, is a pointer to the list of elements in 107 j-th column of matrix V (except pivot elements) */ 108 int *P_row; /* int P_row[1+n]; */ 109 /* P_row[0] is not used; 110 P_row[i] = j means that p[i,j] = 1, where p[i,j] is an element 111 of permutation matrix P */ 112 int *P_col; /* int P_col[1+n]; */ 113 /* P_col[0] is not used; 114 P_col[j] = i means that p[i,j] = 1, where p[i,j] is an element 115 of permutation matrix P */ 116 /* if i-th row or column of matrix F is i'-th row or column of 117 matrix L = P*F*inv(P), or if i-th row of matrix V is i'-th row 118 of matrix U = P*V*Q, then P_row[i'] = i and P_col[i] = i' */ 119 int *Q_row; /* int Q_row[1+n]; */ 120 /* Q_row[0] is not used; 121 Q_row[i] = j means that q[i,j] = 1, where q[i,j] is an element 122 of permutation matrix Q */ 123 int *Q_col; /* int Q_col[1+n]; */ 124 /* Q_col[0] is not used; 125 Q_col[j] = i means that q[i,j] = 1, where q[i,j] is an element 126 of permutation matrix Q */ 127 /* if j-th column of matrix V is j'-th column of matrix U = P*V*Q, 128 then Q_row[j] = j' and Q_col[j'] = j */ 129 int rank; 130 /* the (exact) rank of matrices A and V */ 131 }; 132 133 struct LUXELM 134 { /* element of matrix F or V */ 135 int i; 136 /* row index, 1 <= i <= m */ 137 int j; 138 /* column index, 1 <= j <= n */ 139 mpq_t val; 140 /* numeric (non-zero) element value */ 141 LUXELM *r_prev; 142 /* pointer to previous element in the same row */ 143 LUXELM *r_next; 144 /* pointer to next element in the same row */ 145 LUXELM *c_prev; 146 /* pointer to previous element in the same column */ 147 LUXELM *c_next; 148 /* pointer to next element in the same column */ 149 }; 150 151 struct LUXWKA 152 { /* working area (used only during factorization) */ 153 /* in order to efficiently implement Markowitz strategy and Duff 154 search technique there are two families {R[0], R[1], ..., R[n]} 155 and {C[0], C[1], ..., C[n]}; member R[k] is a set of active 156 rows of matrix V having k non-zeros, and member C[k] is a set 157 of active columns of matrix V having k non-zeros (in the active 158 submatrix); each set R[k] and C[k] is implemented as a separate 159 doubly linked list */ 160 int *R_len; /* int R_len[1+n]; */ 161 /* R_len[0] is not used; 162 R_len[i], 1 <= i <= n, is the number of non-zero elements in 163 i-th row of matrix V (that is the length of i-th row) */ 164 int *R_head; /* int R_head[1+n]; */ 165 /* R_head[k], 0 <= k <= n, is the number of a first row, which is 166 active and whose length is k */ 167 int *R_prev; /* int R_prev[1+n]; */ 168 /* R_prev[0] is not used; 169 R_prev[i], 1 <= i <= n, is the number of a previous row, which 170 is active and has the same length as i-th row */ 171 int *R_next; /* int R_next[1+n]; */ 172 /* R_prev[0] is not used; 173 R_prev[i], 1 <= i <= n, is the number of a next row, which is 174 active and has the same length as i-th row */ 175 int *C_len; /* int C_len[1+n]; */ 176 /* C_len[0] is not used; 177 C_len[j], 1 <= j <= n, is the number of non-zero elements in 178 j-th column of the active submatrix of matrix V (that is the 179 length of j-th column in the active submatrix) */ 180 int *C_head; /* int C_head[1+n]; */ 181 /* C_head[k], 0 <= k <= n, is the number of a first column, which 182 is active and whose length is k */ 183 int *C_prev; /* int C_prev[1+n]; */ 184 /* C_prev[0] is not used; 185 C_prev[j], 1 <= j <= n, is the number of a previous column, 186 which is active and has the same length as j-th column */ 187 int *C_next; /* int C_next[1+n]; */ 188 /* C_next[0] is not used; 189 C_next[j], 1 <= j <= n, is the number of a next column, which 190 is active and has the same length as j-th column */ 191 }; 192 193 #define lux_create _glp_lux_create 194 #define lux_decomp _glp_lux_decomp 195 #define lux_f_solve _glp_lux_f_solve 196 #define lux_v_solve _glp_lux_v_solve 197 #define lux_solve _glp_lux_solve 198 #define lux_delete _glp_lux_delete 199 200 LUX *lux_create(int n); 201 /* create LU-factorization */ 202 203 int lux_decomp(LUX *lux, int (*col)(void *info, int j, int ind[], 204 mpq_t val[]), void *info); 205 /* compute LU-factorization */ 206 207 void lux_f_solve(LUX *lux, int tr, mpq_t x[]); 208 /* solve system F*x = b or F'*x = b */ 209 210 void lux_v_solve(LUX *lux, int tr, mpq_t x[]); 211 /* solve system V*x = b or V'*x = b */ 212 213 void lux_solve(LUX *lux, int tr, mpq_t x[]); 214 /* solve system A*x = b or A'*x = b */ 215 216 void lux_delete(LUX *lux); 217 /* delete LU-factorization */ 218 219 #endif 220 221 /* eof */ 222