[1] | 1 | /* glplux.c */ |
---|
| 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 | #include "glplux.h" |
---|
| 26 | #define xfault xerror |
---|
| 27 | #define dmp_create_poolx(size) dmp_create_pool() |
---|
| 28 | |
---|
| 29 | /*---------------------------------------------------------------------- |
---|
| 30 | // lux_create - create LU-factorization. |
---|
| 31 | // |
---|
| 32 | // SYNOPSIS |
---|
| 33 | // |
---|
| 34 | // #include "glplux.h" |
---|
| 35 | // LUX *lux_create(int n); |
---|
| 36 | // |
---|
| 37 | // DESCRIPTION |
---|
| 38 | // |
---|
| 39 | // The routine lux_create creates LU-factorization data structure for |
---|
| 40 | // a matrix of the order n. Initially the factorization corresponds to |
---|
| 41 | // the unity matrix (F = V = P = Q = I, so A = I). |
---|
| 42 | // |
---|
| 43 | // RETURNS |
---|
| 44 | // |
---|
| 45 | // The routine returns a pointer to the created LU-factorization data |
---|
| 46 | // structure, which represents the unity matrix of the order n. */ |
---|
| 47 | |
---|
| 48 | LUX *lux_create(int n) |
---|
| 49 | { LUX *lux; |
---|
| 50 | int k; |
---|
| 51 | if (n < 1) |
---|
| 52 | xfault("lux_create: n = %d; invalid parameter\n", n); |
---|
| 53 | lux = xmalloc(sizeof(LUX)); |
---|
| 54 | lux->n = n; |
---|
| 55 | lux->pool = dmp_create_poolx(sizeof(LUXELM)); |
---|
| 56 | lux->F_row = xcalloc(1+n, sizeof(LUXELM *)); |
---|
| 57 | lux->F_col = xcalloc(1+n, sizeof(LUXELM *)); |
---|
| 58 | lux->V_piv = xcalloc(1+n, sizeof(mpq_t)); |
---|
| 59 | lux->V_row = xcalloc(1+n, sizeof(LUXELM *)); |
---|
| 60 | lux->V_col = xcalloc(1+n, sizeof(LUXELM *)); |
---|
| 61 | lux->P_row = xcalloc(1+n, sizeof(int)); |
---|
| 62 | lux->P_col = xcalloc(1+n, sizeof(int)); |
---|
| 63 | lux->Q_row = xcalloc(1+n, sizeof(int)); |
---|
| 64 | lux->Q_col = xcalloc(1+n, sizeof(int)); |
---|
| 65 | for (k = 1; k <= n; k++) |
---|
| 66 | { lux->F_row[k] = lux->F_col[k] = NULL; |
---|
| 67 | mpq_init(lux->V_piv[k]); |
---|
| 68 | mpq_set_si(lux->V_piv[k], 1, 1); |
---|
| 69 | lux->V_row[k] = lux->V_col[k] = NULL; |
---|
| 70 | lux->P_row[k] = lux->P_col[k] = k; |
---|
| 71 | lux->Q_row[k] = lux->Q_col[k] = k; |
---|
| 72 | } |
---|
| 73 | lux->rank = n; |
---|
| 74 | return lux; |
---|
| 75 | } |
---|
| 76 | |
---|
| 77 | /*---------------------------------------------------------------------- |
---|
| 78 | // initialize - initialize LU-factorization data structures. |
---|
| 79 | // |
---|
| 80 | // This routine initializes data structures for subsequent computing |
---|
| 81 | // the LU-factorization of a given matrix A, which is specified by the |
---|
| 82 | // formal routine col. On exit V = A and F = P = Q = I, where I is the |
---|
| 83 | // unity matrix. */ |
---|
| 84 | |
---|
| 85 | static void initialize(LUX *lux, int (*col)(void *info, int j, |
---|
| 86 | int ind[], mpq_t val[]), void *info, LUXWKA *wka) |
---|
| 87 | { int n = lux->n; |
---|
| 88 | DMP *pool = lux->pool; |
---|
| 89 | LUXELM **F_row = lux->F_row; |
---|
| 90 | LUXELM **F_col = lux->F_col; |
---|
| 91 | mpq_t *V_piv = lux->V_piv; |
---|
| 92 | LUXELM **V_row = lux->V_row; |
---|
| 93 | LUXELM **V_col = lux->V_col; |
---|
| 94 | int *P_row = lux->P_row; |
---|
| 95 | int *P_col = lux->P_col; |
---|
| 96 | int *Q_row = lux->Q_row; |
---|
| 97 | int *Q_col = lux->Q_col; |
---|
| 98 | int *R_len = wka->R_len; |
---|
| 99 | int *R_head = wka->R_head; |
---|
| 100 | int *R_prev = wka->R_prev; |
---|
| 101 | int *R_next = wka->R_next; |
---|
| 102 | int *C_len = wka->C_len; |
---|
| 103 | int *C_head = wka->C_head; |
---|
| 104 | int *C_prev = wka->C_prev; |
---|
| 105 | int *C_next = wka->C_next; |
---|
| 106 | LUXELM *fij, *vij; |
---|
| 107 | int i, j, k, len, *ind; |
---|
| 108 | mpq_t *val; |
---|
| 109 | /* F := I */ |
---|
| 110 | for (i = 1; i <= n; i++) |
---|
| 111 | { while (F_row[i] != NULL) |
---|
| 112 | { fij = F_row[i], F_row[i] = fij->r_next; |
---|
| 113 | mpq_clear(fij->val); |
---|
| 114 | dmp_free_atom(pool, fij, sizeof(LUXELM)); |
---|
| 115 | } |
---|
| 116 | } |
---|
| 117 | for (j = 1; j <= n; j++) F_col[j] = NULL; |
---|
| 118 | /* V := 0 */ |
---|
| 119 | for (k = 1; k <= n; k++) mpq_set_si(V_piv[k], 0, 1); |
---|
| 120 | for (i = 1; i <= n; i++) |
---|
| 121 | { while (V_row[i] != NULL) |
---|
| 122 | { vij = V_row[i], V_row[i] = vij->r_next; |
---|
| 123 | mpq_clear(vij->val); |
---|
| 124 | dmp_free_atom(pool, vij, sizeof(LUXELM)); |
---|
| 125 | } |
---|
| 126 | } |
---|
| 127 | for (j = 1; j <= n; j++) V_col[j] = NULL; |
---|
| 128 | /* V := A */ |
---|
| 129 | ind = xcalloc(1+n, sizeof(int)); |
---|
| 130 | val = xcalloc(1+n, sizeof(mpq_t)); |
---|
| 131 | for (k = 1; k <= n; k++) mpq_init(val[k]); |
---|
| 132 | for (j = 1; j <= n; j++) |
---|
| 133 | { /* obtain j-th column of matrix A */ |
---|
| 134 | len = col(info, j, ind, val); |
---|
| 135 | if (!(0 <= len && len <= n)) |
---|
| 136 | xfault("lux_decomp: j = %d: len = %d; invalid column length" |
---|
| 137 | "\n", j, len); |
---|
| 138 | /* copy elements of j-th column to matrix V */ |
---|
| 139 | for (k = 1; k <= len; k++) |
---|
| 140 | { /* get row index of a[i,j] */ |
---|
| 141 | i = ind[k]; |
---|
| 142 | if (!(1 <= i && i <= n)) |
---|
| 143 | xfault("lux_decomp: j = %d: i = %d; row index out of ran" |
---|
| 144 | "ge\n", j, i); |
---|
| 145 | /* check for duplicate indices */ |
---|
| 146 | if (V_row[i] != NULL && V_row[i]->j == j) |
---|
| 147 | xfault("lux_decomp: j = %d: i = %d; duplicate row indice" |
---|
| 148 | "s not allowed\n", j, i); |
---|
| 149 | /* check for zero value */ |
---|
| 150 | if (mpq_sgn(val[k]) == 0) |
---|
| 151 | xfault("lux_decomp: j = %d: i = %d; zero elements not al" |
---|
| 152 | "lowed\n", j, i); |
---|
| 153 | /* add new element v[i,j] = a[i,j] to V */ |
---|
| 154 | vij = dmp_get_atom(pool, sizeof(LUXELM)); |
---|
| 155 | vij->i = i, vij->j = j; |
---|
| 156 | mpq_init(vij->val); |
---|
| 157 | mpq_set(vij->val, val[k]); |
---|
| 158 | vij->r_prev = NULL; |
---|
| 159 | vij->r_next = V_row[i]; |
---|
| 160 | vij->c_prev = NULL; |
---|
| 161 | vij->c_next = V_col[j]; |
---|
| 162 | if (vij->r_next != NULL) vij->r_next->r_prev = vij; |
---|
| 163 | if (vij->c_next != NULL) vij->c_next->c_prev = vij; |
---|
| 164 | V_row[i] = V_col[j] = vij; |
---|
| 165 | } |
---|
| 166 | } |
---|
| 167 | xfree(ind); |
---|
| 168 | for (k = 1; k <= n; k++) mpq_clear(val[k]); |
---|
| 169 | xfree(val); |
---|
| 170 | /* P := Q := I */ |
---|
| 171 | for (k = 1; k <= n; k++) |
---|
| 172 | P_row[k] = P_col[k] = Q_row[k] = Q_col[k] = k; |
---|
| 173 | /* the rank of A and V is not determined yet */ |
---|
| 174 | lux->rank = -1; |
---|
| 175 | /* initially the entire matrix V is active */ |
---|
| 176 | /* determine its row lengths */ |
---|
| 177 | for (i = 1; i <= n; i++) |
---|
| 178 | { len = 0; |
---|
| 179 | for (vij = V_row[i]; vij != NULL; vij = vij->r_next) len++; |
---|
| 180 | R_len[i] = len; |
---|
| 181 | } |
---|
| 182 | /* build linked lists of active rows */ |
---|
| 183 | for (len = 0; len <= n; len++) R_head[len] = 0; |
---|
| 184 | for (i = 1; i <= n; i++) |
---|
| 185 | { len = R_len[i]; |
---|
| 186 | R_prev[i] = 0; |
---|
| 187 | R_next[i] = R_head[len]; |
---|
| 188 | if (R_next[i] != 0) R_prev[R_next[i]] = i; |
---|
| 189 | R_head[len] = i; |
---|
| 190 | } |
---|
| 191 | /* determine its column lengths */ |
---|
| 192 | for (j = 1; j <= n; j++) |
---|
| 193 | { len = 0; |
---|
| 194 | for (vij = V_col[j]; vij != NULL; vij = vij->c_next) len++; |
---|
| 195 | C_len[j] = len; |
---|
| 196 | } |
---|
| 197 | /* build linked lists of active columns */ |
---|
| 198 | for (len = 0; len <= n; len++) C_head[len] = 0; |
---|
| 199 | for (j = 1; j <= n; j++) |
---|
| 200 | { len = C_len[j]; |
---|
| 201 | C_prev[j] = 0; |
---|
| 202 | C_next[j] = C_head[len]; |
---|
| 203 | if (C_next[j] != 0) C_prev[C_next[j]] = j; |
---|
| 204 | C_head[len] = j; |
---|
| 205 | } |
---|
| 206 | return; |
---|
| 207 | } |
---|
| 208 | |
---|
| 209 | /*---------------------------------------------------------------------- |
---|
| 210 | // find_pivot - choose a pivot element. |
---|
| 211 | // |
---|
| 212 | // This routine chooses a pivot element v[p,q] in the active submatrix |
---|
| 213 | // of matrix U = P*V*Q. |
---|
| 214 | // |
---|
| 215 | // It is assumed that on entry the matrix U has the following partially |
---|
| 216 | // triangularized form: |
---|
| 217 | // |
---|
| 218 | // 1 k n |
---|
| 219 | // 1 x x x x x x x x x x |
---|
| 220 | // . x x x x x x x x x |
---|
| 221 | // . . x x x x x x x x |
---|
| 222 | // . . . x x x x x x x |
---|
| 223 | // k . . . . * * * * * * |
---|
| 224 | // . . . . * * * * * * |
---|
| 225 | // . . . . * * * * * * |
---|
| 226 | // . . . . * * * * * * |
---|
| 227 | // . . . . * * * * * * |
---|
| 228 | // n . . . . * * * * * * |
---|
| 229 | // |
---|
| 230 | // where rows and columns k, k+1, ..., n belong to the active submatrix |
---|
| 231 | // (elements of the active submatrix are marked by '*'). |
---|
| 232 | // |
---|
| 233 | // Since the matrix U = P*V*Q is not stored, the routine works with the |
---|
| 234 | // matrix V. It is assumed that the row-wise representation corresponds |
---|
| 235 | // to the matrix V, but the column-wise representation corresponds to |
---|
| 236 | // the active submatrix of the matrix V, i.e. elements of the matrix V, |
---|
| 237 | // which does not belong to the active submatrix, are missing from the |
---|
| 238 | // column linked lists. It is also assumed that each active row of the |
---|
| 239 | // matrix V is in the set R[len], where len is number of non-zeros in |
---|
| 240 | // the row, and each active column of the matrix V is in the set C[len], |
---|
| 241 | // where len is number of non-zeros in the column (in the latter case |
---|
| 242 | // only elements of the active submatrix are counted; such elements are |
---|
| 243 | // marked by '*' on the figure above). |
---|
| 244 | // |
---|
| 245 | // Due to exact arithmetic any non-zero element of the active submatrix |
---|
| 246 | // can be chosen as a pivot. However, to keep sparsity of the matrix V |
---|
| 247 | // the routine uses Markowitz strategy, trying to choose such element |
---|
| 248 | // v[p,q], which has smallest Markowitz cost (nr[p]-1) * (nc[q]-1), |
---|
| 249 | // where nr[p] and nc[q] are the number of non-zero elements, resp., in |
---|
| 250 | // p-th row and in q-th column of the active submatrix. |
---|
| 251 | // |
---|
| 252 | // In order to reduce the search, i.e. not to walk through all elements |
---|
| 253 | // of the active submatrix, the routine exploits a technique proposed by |
---|
| 254 | // I.Duff. This technique is based on using the sets R[len] and C[len] |
---|
| 255 | // of active rows and columns. |
---|
| 256 | // |
---|
| 257 | // On exit the routine returns a pointer to a pivot v[p,q] chosen, or |
---|
| 258 | // NULL, if the active submatrix is empty. */ |
---|
| 259 | |
---|
| 260 | static LUXELM *find_pivot(LUX *lux, LUXWKA *wka) |
---|
| 261 | { int n = lux->n; |
---|
| 262 | LUXELM **V_row = lux->V_row; |
---|
| 263 | LUXELM **V_col = lux->V_col; |
---|
| 264 | int *R_len = wka->R_len; |
---|
| 265 | int *R_head = wka->R_head; |
---|
| 266 | int *R_next = wka->R_next; |
---|
| 267 | int *C_len = wka->C_len; |
---|
| 268 | int *C_head = wka->C_head; |
---|
| 269 | int *C_next = wka->C_next; |
---|
| 270 | LUXELM *piv, *some, *vij; |
---|
| 271 | int i, j, len, min_len, ncand, piv_lim = 5; |
---|
| 272 | double best, cost; |
---|
| 273 | /* nothing is chosen so far */ |
---|
| 274 | piv = NULL, best = DBL_MAX, ncand = 0; |
---|
| 275 | /* if in the active submatrix there is a column that has the only |
---|
| 276 | non-zero (column singleton), choose it as a pivot */ |
---|
| 277 | j = C_head[1]; |
---|
| 278 | if (j != 0) |
---|
| 279 | { xassert(C_len[j] == 1); |
---|
| 280 | piv = V_col[j]; |
---|
| 281 | xassert(piv != NULL && piv->c_next == NULL); |
---|
| 282 | goto done; |
---|
| 283 | } |
---|
| 284 | /* if in the active submatrix there is a row that has the only |
---|
| 285 | non-zero (row singleton), choose it as a pivot */ |
---|
| 286 | i = R_head[1]; |
---|
| 287 | if (i != 0) |
---|
| 288 | { xassert(R_len[i] == 1); |
---|
| 289 | piv = V_row[i]; |
---|
| 290 | xassert(piv != NULL && piv->r_next == NULL); |
---|
| 291 | goto done; |
---|
| 292 | } |
---|
| 293 | /* there are no singletons in the active submatrix; walk through |
---|
| 294 | other non-empty rows and columns */ |
---|
| 295 | for (len = 2; len <= n; len++) |
---|
| 296 | { /* consider active columns having len non-zeros */ |
---|
| 297 | for (j = C_head[len]; j != 0; j = C_next[j]) |
---|
| 298 | { /* j-th column has len non-zeros */ |
---|
| 299 | /* find an element in the row of minimal length */ |
---|
| 300 | some = NULL, min_len = INT_MAX; |
---|
| 301 | for (vij = V_col[j]; vij != NULL; vij = vij->c_next) |
---|
| 302 | { if (min_len > R_len[vij->i]) |
---|
| 303 | some = vij, min_len = R_len[vij->i]; |
---|
| 304 | /* if Markowitz cost of this element is not greater than |
---|
| 305 | (len-1)**2, it can be chosen right now; this heuristic |
---|
| 306 | reduces the search and works well in many cases */ |
---|
| 307 | if (min_len <= len) |
---|
| 308 | { piv = some; |
---|
| 309 | goto done; |
---|
| 310 | } |
---|
| 311 | } |
---|
| 312 | /* j-th column has been scanned */ |
---|
| 313 | /* the minimal element found is a next pivot candidate */ |
---|
| 314 | xassert(some != NULL); |
---|
| 315 | ncand++; |
---|
| 316 | /* compute its Markowitz cost */ |
---|
| 317 | cost = (double)(min_len - 1) * (double)(len - 1); |
---|
| 318 | /* choose between the current candidate and this element */ |
---|
| 319 | if (cost < best) piv = some, best = cost; |
---|
| 320 | /* if piv_lim candidates have been considered, there is a |
---|
| 321 | doubt that a much better candidate exists; therefore it |
---|
| 322 | is the time to terminate the search */ |
---|
| 323 | if (ncand == piv_lim) goto done; |
---|
| 324 | } |
---|
| 325 | /* now consider active rows having len non-zeros */ |
---|
| 326 | for (i = R_head[len]; i != 0; i = R_next[i]) |
---|
| 327 | { /* i-th row has len non-zeros */ |
---|
| 328 | /* find an element in the column of minimal length */ |
---|
| 329 | some = NULL, min_len = INT_MAX; |
---|
| 330 | for (vij = V_row[i]; vij != NULL; vij = vij->r_next) |
---|
| 331 | { if (min_len > C_len[vij->j]) |
---|
| 332 | some = vij, min_len = C_len[vij->j]; |
---|
| 333 | /* if Markowitz cost of this element is not greater than |
---|
| 334 | (len-1)**2, it can be chosen right now; this heuristic |
---|
| 335 | reduces the search and works well in many cases */ |
---|
| 336 | if (min_len <= len) |
---|
| 337 | { piv = some; |
---|
| 338 | goto done; |
---|
| 339 | } |
---|
| 340 | } |
---|
| 341 | /* i-th row has been scanned */ |
---|
| 342 | /* the minimal element found is a next pivot candidate */ |
---|
| 343 | xassert(some != NULL); |
---|
| 344 | ncand++; |
---|
| 345 | /* compute its Markowitz cost */ |
---|
| 346 | cost = (double)(len - 1) * (double)(min_len - 1); |
---|
| 347 | /* choose between the current candidate and this element */ |
---|
| 348 | if (cost < best) piv = some, best = cost; |
---|
| 349 | /* if piv_lim candidates have been considered, there is a |
---|
| 350 | doubt that a much better candidate exists; therefore it |
---|
| 351 | is the time to terminate the search */ |
---|
| 352 | if (ncand == piv_lim) goto done; |
---|
| 353 | } |
---|
| 354 | } |
---|
| 355 | done: /* bring the pivot v[p,q] to the factorizing routine */ |
---|
| 356 | return piv; |
---|
| 357 | } |
---|
| 358 | |
---|
| 359 | /*---------------------------------------------------------------------- |
---|
| 360 | // eliminate - perform gaussian elimination. |
---|
| 361 | // |
---|
| 362 | // This routine performs elementary gaussian transformations in order |
---|
| 363 | // to eliminate subdiagonal elements in the k-th column of the matrix |
---|
| 364 | // U = P*V*Q using the pivot element u[k,k], where k is the number of |
---|
| 365 | // the current elimination step. |
---|
| 366 | // |
---|
| 367 | // The parameter piv specifies the pivot element v[p,q] = u[k,k]. |
---|
| 368 | // |
---|
| 369 | // Each time when the routine applies the elementary transformation to |
---|
| 370 | // a non-pivot row of the matrix V, it stores the corresponding element |
---|
| 371 | // to the matrix F in order to keep the main equality A = F*V. |
---|
| 372 | // |
---|
| 373 | // The routine assumes that on entry the matrices L = P*F*inv(P) and |
---|
| 374 | // U = P*V*Q are the following: |
---|
| 375 | // |
---|
| 376 | // 1 k 1 k n |
---|
| 377 | // 1 1 . . . . . . . . . 1 x x x x x x x x x x |
---|
| 378 | // x 1 . . . . . . . . . x x x x x x x x x |
---|
| 379 | // x x 1 . . . . . . . . . x x x x x x x x |
---|
| 380 | // x x x 1 . . . . . . . . . x x x x x x x |
---|
| 381 | // k x x x x 1 . . . . . k . . . . * * * * * * |
---|
| 382 | // x x x x _ 1 . . . . . . . . # * * * * * |
---|
| 383 | // x x x x _ . 1 . . . . . . . # * * * * * |
---|
| 384 | // x x x x _ . . 1 . . . . . . # * * * * * |
---|
| 385 | // x x x x _ . . . 1 . . . . . # * * * * * |
---|
| 386 | // n x x x x _ . . . . 1 n . . . . # * * * * * |
---|
| 387 | // |
---|
| 388 | // matrix L matrix U |
---|
| 389 | // |
---|
| 390 | // where rows and columns of the matrix U with numbers k, k+1, ..., n |
---|
| 391 | // form the active submatrix (eliminated elements are marked by '#' and |
---|
| 392 | // other elements of the active submatrix are marked by '*'). Note that |
---|
| 393 | // each eliminated non-zero element u[i,k] of the matrix U gives the |
---|
| 394 | // corresponding element l[i,k] of the matrix L (marked by '_'). |
---|
| 395 | // |
---|
| 396 | // Actually all operations are performed on the matrix V. Should note |
---|
| 397 | // that the row-wise representation corresponds to the matrix V, but the |
---|
| 398 | // column-wise representation corresponds to the active submatrix of the |
---|
| 399 | // matrix V, i.e. elements of the matrix V, which doesn't belong to the |
---|
| 400 | // active submatrix, are missing from the column linked lists. |
---|
| 401 | // |
---|
| 402 | // Let u[k,k] = v[p,q] be the pivot. In order to eliminate subdiagonal |
---|
| 403 | // elements u[i',k] = v[i,q], i' = k+1, k+2, ..., n, the routine applies |
---|
| 404 | // the following elementary gaussian transformations: |
---|
| 405 | // |
---|
| 406 | // (i-th row of V) := (i-th row of V) - f[i,p] * (p-th row of V), |
---|
| 407 | // |
---|
| 408 | // where f[i,p] = v[i,q] / v[p,q] is a gaussian multiplier. |
---|
| 409 | // |
---|
| 410 | // Additionally, in order to keep the main equality A = F*V, each time |
---|
| 411 | // when the routine applies the transformation to i-th row of the matrix |
---|
| 412 | // V, it also adds f[i,p] as a new element to the matrix F. |
---|
| 413 | // |
---|
| 414 | // IMPORTANT: On entry the working arrays flag and work should contain |
---|
| 415 | // zeros. This status is provided by the routine on exit. */ |
---|
| 416 | |
---|
| 417 | static void eliminate(LUX *lux, LUXWKA *wka, LUXELM *piv, int flag[], |
---|
| 418 | mpq_t work[]) |
---|
| 419 | { DMP *pool = lux->pool; |
---|
| 420 | LUXELM **F_row = lux->F_row; |
---|
| 421 | LUXELM **F_col = lux->F_col; |
---|
| 422 | mpq_t *V_piv = lux->V_piv; |
---|
| 423 | LUXELM **V_row = lux->V_row; |
---|
| 424 | LUXELM **V_col = lux->V_col; |
---|
| 425 | int *R_len = wka->R_len; |
---|
| 426 | int *R_head = wka->R_head; |
---|
| 427 | int *R_prev = wka->R_prev; |
---|
| 428 | int *R_next = wka->R_next; |
---|
| 429 | int *C_len = wka->C_len; |
---|
| 430 | int *C_head = wka->C_head; |
---|
| 431 | int *C_prev = wka->C_prev; |
---|
| 432 | int *C_next = wka->C_next; |
---|
| 433 | LUXELM *fip, *vij, *vpj, *viq, *next; |
---|
| 434 | mpq_t temp; |
---|
| 435 | int i, j, p, q; |
---|
| 436 | mpq_init(temp); |
---|
| 437 | /* determine row and column indices of the pivot v[p,q] */ |
---|
| 438 | xassert(piv != NULL); |
---|
| 439 | p = piv->i, q = piv->j; |
---|
| 440 | /* remove p-th (pivot) row from the active set; it will never |
---|
| 441 | return there */ |
---|
| 442 | if (R_prev[p] == 0) |
---|
| 443 | R_head[R_len[p]] = R_next[p]; |
---|
| 444 | else |
---|
| 445 | R_next[R_prev[p]] = R_next[p]; |
---|
| 446 | if (R_next[p] == 0) |
---|
| 447 | ; |
---|
| 448 | else |
---|
| 449 | R_prev[R_next[p]] = R_prev[p]; |
---|
| 450 | /* remove q-th (pivot) column from the active set; it will never |
---|
| 451 | return there */ |
---|
| 452 | if (C_prev[q] == 0) |
---|
| 453 | C_head[C_len[q]] = C_next[q]; |
---|
| 454 | else |
---|
| 455 | C_next[C_prev[q]] = C_next[q]; |
---|
| 456 | if (C_next[q] == 0) |
---|
| 457 | ; |
---|
| 458 | else |
---|
| 459 | C_prev[C_next[q]] = C_prev[q]; |
---|
| 460 | /* store the pivot value in a separate array */ |
---|
| 461 | mpq_set(V_piv[p], piv->val); |
---|
| 462 | /* remove the pivot from p-th row */ |
---|
| 463 | if (piv->r_prev == NULL) |
---|
| 464 | V_row[p] = piv->r_next; |
---|
| 465 | else |
---|
| 466 | piv->r_prev->r_next = piv->r_next; |
---|
| 467 | if (piv->r_next == NULL) |
---|
| 468 | ; |
---|
| 469 | else |
---|
| 470 | piv->r_next->r_prev = piv->r_prev; |
---|
| 471 | R_len[p]--; |
---|
| 472 | /* remove the pivot from q-th column */ |
---|
| 473 | if (piv->c_prev == NULL) |
---|
| 474 | V_col[q] = piv->c_next; |
---|
| 475 | else |
---|
| 476 | piv->c_prev->c_next = piv->c_next; |
---|
| 477 | if (piv->c_next == NULL) |
---|
| 478 | ; |
---|
| 479 | else |
---|
| 480 | piv->c_next->c_prev = piv->c_prev; |
---|
| 481 | C_len[q]--; |
---|
| 482 | /* free the space occupied by the pivot */ |
---|
| 483 | mpq_clear(piv->val); |
---|
| 484 | dmp_free_atom(pool, piv, sizeof(LUXELM)); |
---|
| 485 | /* walk through p-th (pivot) row, which already does not contain |
---|
| 486 | the pivot v[p,q], and do the following... */ |
---|
| 487 | for (vpj = V_row[p]; vpj != NULL; vpj = vpj->r_next) |
---|
| 488 | { /* get column index of v[p,j] */ |
---|
| 489 | j = vpj->j; |
---|
| 490 | /* store v[p,j] in the working array */ |
---|
| 491 | flag[j] = 1; |
---|
| 492 | mpq_set(work[j], vpj->val); |
---|
| 493 | /* remove j-th column from the active set; it will return there |
---|
| 494 | later with a new length */ |
---|
| 495 | if (C_prev[j] == 0) |
---|
| 496 | C_head[C_len[j]] = C_next[j]; |
---|
| 497 | else |
---|
| 498 | C_next[C_prev[j]] = C_next[j]; |
---|
| 499 | if (C_next[j] == 0) |
---|
| 500 | ; |
---|
| 501 | else |
---|
| 502 | C_prev[C_next[j]] = C_prev[j]; |
---|
| 503 | /* v[p,j] leaves the active submatrix, so remove it from j-th |
---|
| 504 | column; however, v[p,j] is kept in p-th row */ |
---|
| 505 | if (vpj->c_prev == NULL) |
---|
| 506 | V_col[j] = vpj->c_next; |
---|
| 507 | else |
---|
| 508 | vpj->c_prev->c_next = vpj->c_next; |
---|
| 509 | if (vpj->c_next == NULL) |
---|
| 510 | ; |
---|
| 511 | else |
---|
| 512 | vpj->c_next->c_prev = vpj->c_prev; |
---|
| 513 | C_len[j]--; |
---|
| 514 | } |
---|
| 515 | /* now walk through q-th (pivot) column, which already does not |
---|
| 516 | contain the pivot v[p,q], and perform gaussian elimination */ |
---|
| 517 | while (V_col[q] != NULL) |
---|
| 518 | { /* element v[i,q] has to be eliminated */ |
---|
| 519 | viq = V_col[q]; |
---|
| 520 | /* get row index of v[i,q] */ |
---|
| 521 | i = viq->i; |
---|
| 522 | /* remove i-th row from the active set; later it will return |
---|
| 523 | there with a new length */ |
---|
| 524 | if (R_prev[i] == 0) |
---|
| 525 | R_head[R_len[i]] = R_next[i]; |
---|
| 526 | else |
---|
| 527 | R_next[R_prev[i]] = R_next[i]; |
---|
| 528 | if (R_next[i] == 0) |
---|
| 529 | ; |
---|
| 530 | else |
---|
| 531 | R_prev[R_next[i]] = R_prev[i]; |
---|
| 532 | /* compute gaussian multiplier f[i,p] = v[i,q] / v[p,q] and |
---|
| 533 | store it in the matrix F */ |
---|
| 534 | fip = dmp_get_atom(pool, sizeof(LUXELM)); |
---|
| 535 | fip->i = i, fip->j = p; |
---|
| 536 | mpq_init(fip->val); |
---|
| 537 | mpq_div(fip->val, viq->val, V_piv[p]); |
---|
| 538 | fip->r_prev = NULL; |
---|
| 539 | fip->r_next = F_row[i]; |
---|
| 540 | fip->c_prev = NULL; |
---|
| 541 | fip->c_next = F_col[p]; |
---|
| 542 | if (fip->r_next != NULL) fip->r_next->r_prev = fip; |
---|
| 543 | if (fip->c_next != NULL) fip->c_next->c_prev = fip; |
---|
| 544 | F_row[i] = F_col[p] = fip; |
---|
| 545 | /* v[i,q] has to be eliminated, so remove it from i-th row */ |
---|
| 546 | if (viq->r_prev == NULL) |
---|
| 547 | V_row[i] = viq->r_next; |
---|
| 548 | else |
---|
| 549 | viq->r_prev->r_next = viq->r_next; |
---|
| 550 | if (viq->r_next == NULL) |
---|
| 551 | ; |
---|
| 552 | else |
---|
| 553 | viq->r_next->r_prev = viq->r_prev; |
---|
| 554 | R_len[i]--; |
---|
| 555 | /* and also from q-th column */ |
---|
| 556 | V_col[q] = viq->c_next; |
---|
| 557 | C_len[q]--; |
---|
| 558 | /* free the space occupied by v[i,q] */ |
---|
| 559 | mpq_clear(viq->val); |
---|
| 560 | dmp_free_atom(pool, viq, sizeof(LUXELM)); |
---|
| 561 | /* perform gaussian transformation: |
---|
| 562 | (i-th row) := (i-th row) - f[i,p] * (p-th row) |
---|
| 563 | note that now p-th row, which is in the working array, |
---|
| 564 | does not contain the pivot v[p,q], and i-th row does not |
---|
| 565 | contain the element v[i,q] to be eliminated */ |
---|
| 566 | /* walk through i-th row and transform existing non-zero |
---|
| 567 | elements */ |
---|
| 568 | for (vij = V_row[i]; vij != NULL; vij = next) |
---|
| 569 | { next = vij->r_next; |
---|
| 570 | /* get column index of v[i,j] */ |
---|
| 571 | j = vij->j; |
---|
| 572 | /* v[i,j] := v[i,j] - f[i,p] * v[p,j] */ |
---|
| 573 | if (flag[j]) |
---|
| 574 | { /* v[p,j] != 0 */ |
---|
| 575 | flag[j] = 0; |
---|
| 576 | mpq_mul(temp, fip->val, work[j]); |
---|
| 577 | mpq_sub(vij->val, vij->val, temp); |
---|
| 578 | if (mpq_sgn(vij->val) == 0) |
---|
| 579 | { /* new v[i,j] is zero, so remove it from the active |
---|
| 580 | submatrix */ |
---|
| 581 | /* remove v[i,j] from i-th row */ |
---|
| 582 | if (vij->r_prev == NULL) |
---|
| 583 | V_row[i] = vij->r_next; |
---|
| 584 | else |
---|
| 585 | vij->r_prev->r_next = vij->r_next; |
---|
| 586 | if (vij->r_next == NULL) |
---|
| 587 | ; |
---|
| 588 | else |
---|
| 589 | vij->r_next->r_prev = vij->r_prev; |
---|
| 590 | R_len[i]--; |
---|
| 591 | /* remove v[i,j] from j-th column */ |
---|
| 592 | if (vij->c_prev == NULL) |
---|
| 593 | V_col[j] = vij->c_next; |
---|
| 594 | else |
---|
| 595 | vij->c_prev->c_next = vij->c_next; |
---|
| 596 | if (vij->c_next == NULL) |
---|
| 597 | ; |
---|
| 598 | else |
---|
| 599 | vij->c_next->c_prev = vij->c_prev; |
---|
| 600 | C_len[j]--; |
---|
| 601 | /* free the space occupied by v[i,j] */ |
---|
| 602 | mpq_clear(vij->val); |
---|
| 603 | dmp_free_atom(pool, vij, sizeof(LUXELM)); |
---|
| 604 | } |
---|
| 605 | } |
---|
| 606 | } |
---|
| 607 | /* now flag is the pattern of the set v[p,*] \ v[i,*] */ |
---|
| 608 | /* walk through p-th (pivot) row and create new elements in |
---|
| 609 | i-th row, which appear due to fill-in */ |
---|
| 610 | for (vpj = V_row[p]; vpj != NULL; vpj = vpj->r_next) |
---|
| 611 | { j = vpj->j; |
---|
| 612 | if (flag[j]) |
---|
| 613 | { /* create new non-zero v[i,j] = 0 - f[i,p] * v[p,j] and |
---|
| 614 | add it to i-th row and j-th column */ |
---|
| 615 | vij = dmp_get_atom(pool, sizeof(LUXELM)); |
---|
| 616 | vij->i = i, vij->j = j; |
---|
| 617 | mpq_init(vij->val); |
---|
| 618 | mpq_mul(vij->val, fip->val, work[j]); |
---|
| 619 | mpq_neg(vij->val, vij->val); |
---|
| 620 | vij->r_prev = NULL; |
---|
| 621 | vij->r_next = V_row[i]; |
---|
| 622 | vij->c_prev = NULL; |
---|
| 623 | vij->c_next = V_col[j]; |
---|
| 624 | if (vij->r_next != NULL) vij->r_next->r_prev = vij; |
---|
| 625 | if (vij->c_next != NULL) vij->c_next->c_prev = vij; |
---|
| 626 | V_row[i] = V_col[j] = vij; |
---|
| 627 | R_len[i]++, C_len[j]++; |
---|
| 628 | } |
---|
| 629 | else |
---|
| 630 | { /* there is no fill-in, because v[i,j] already exists in |
---|
| 631 | i-th row; restore the flag, which was reset before */ |
---|
| 632 | flag[j] = 1; |
---|
| 633 | } |
---|
| 634 | } |
---|
| 635 | /* now i-th row has been completely transformed and can return |
---|
| 636 | to the active set with a new length */ |
---|
| 637 | R_prev[i] = 0; |
---|
| 638 | R_next[i] = R_head[R_len[i]]; |
---|
| 639 | if (R_next[i] != 0) R_prev[R_next[i]] = i; |
---|
| 640 | R_head[R_len[i]] = i; |
---|
| 641 | } |
---|
| 642 | /* at this point q-th (pivot) column must be empty */ |
---|
| 643 | xassert(C_len[q] == 0); |
---|
| 644 | /* walk through p-th (pivot) row again and do the following... */ |
---|
| 645 | for (vpj = V_row[p]; vpj != NULL; vpj = vpj->r_next) |
---|
| 646 | { /* get column index of v[p,j] */ |
---|
| 647 | j = vpj->j; |
---|
| 648 | /* erase v[p,j] from the working array */ |
---|
| 649 | flag[j] = 0; |
---|
| 650 | mpq_set_si(work[j], 0, 1); |
---|
| 651 | /* now j-th column has been completely transformed, so it can |
---|
| 652 | return to the active list with a new length */ |
---|
| 653 | C_prev[j] = 0; |
---|
| 654 | C_next[j] = C_head[C_len[j]]; |
---|
| 655 | if (C_next[j] != 0) C_prev[C_next[j]] = j; |
---|
| 656 | C_head[C_len[j]] = j; |
---|
| 657 | } |
---|
| 658 | mpq_clear(temp); |
---|
| 659 | /* return to the factorizing routine */ |
---|
| 660 | return; |
---|
| 661 | } |
---|
| 662 | |
---|
| 663 | /*---------------------------------------------------------------------- |
---|
| 664 | // lux_decomp - compute LU-factorization. |
---|
| 665 | // |
---|
| 666 | // SYNOPSIS |
---|
| 667 | // |
---|
| 668 | // #include "glplux.h" |
---|
| 669 | // int lux_decomp(LUX *lux, int (*col)(void *info, int j, int ind[], |
---|
| 670 | // mpq_t val[]), void *info); |
---|
| 671 | // |
---|
| 672 | // DESCRIPTION |
---|
| 673 | // |
---|
| 674 | // The routine lux_decomp computes LU-factorization of a given square |
---|
| 675 | // matrix A. |
---|
| 676 | // |
---|
| 677 | // The parameter lux specifies LU-factorization data structure built by |
---|
| 678 | // means of the routine lux_create. |
---|
| 679 | // |
---|
| 680 | // The formal routine col specifies the original matrix A. In order to |
---|
| 681 | // obtain j-th column of the matrix A the routine lux_decomp calls the |
---|
| 682 | // routine col with the parameter j (1 <= j <= n, where n is the order |
---|
| 683 | // of A). In response the routine col should store row indices and |
---|
| 684 | // numerical values of non-zero elements of j-th column of A to the |
---|
| 685 | // locations ind[1], ..., ind[len] and val[1], ..., val[len], resp., |
---|
| 686 | // where len is the number of non-zeros in j-th column, which should be |
---|
| 687 | // returned on exit. Neiter zero nor duplicate elements are allowed. |
---|
| 688 | // |
---|
| 689 | // The parameter info is a transit pointer passed to the formal routine |
---|
| 690 | // col; it can be used for various purposes. |
---|
| 691 | // |
---|
| 692 | // RETURNS |
---|
| 693 | // |
---|
| 694 | // The routine lux_decomp returns the singularity flag. Zero flag means |
---|
| 695 | // that the original matrix A is non-singular while non-zero flag means |
---|
| 696 | // that A is (exactly!) singular. |
---|
| 697 | // |
---|
| 698 | // Note that LU-factorization is valid in both cases, however, in case |
---|
| 699 | // of singularity some rows of the matrix V (including pivot elements) |
---|
| 700 | // will be empty. |
---|
| 701 | // |
---|
| 702 | // REPAIRING SINGULAR MATRIX |
---|
| 703 | // |
---|
| 704 | // If the routine lux_decomp returns non-zero flag, it provides all |
---|
| 705 | // necessary information that can be used for "repairing" the matrix A, |
---|
| 706 | // where "repairing" means replacing linearly dependent columns of the |
---|
| 707 | // matrix A by appropriate columns of the unity matrix. This feature is |
---|
| 708 | // needed when the routine lux_decomp is used for reinverting the basis |
---|
| 709 | // matrix within the simplex method procedure. |
---|
| 710 | // |
---|
| 711 | // On exit linearly dependent columns of the matrix U have the numbers |
---|
| 712 | // rank+1, rank+2, ..., n, where rank is the exact rank of the matrix A |
---|
| 713 | // stored by the routine to the member lux->rank. The correspondence |
---|
| 714 | // between columns of A and U is the same as between columns of V and U. |
---|
| 715 | // Thus, linearly dependent columns of the matrix A have the numbers |
---|
| 716 | // Q_col[rank+1], Q_col[rank+2], ..., Q_col[n], where Q_col is an array |
---|
| 717 | // representing the permutation matrix Q in column-like format. It is |
---|
| 718 | // understood that each j-th linearly dependent column of the matrix U |
---|
| 719 | // should be replaced by the unity vector, where all elements are zero |
---|
| 720 | // except the unity diagonal element u[j,j]. On the other hand j-th row |
---|
| 721 | // of the matrix U corresponds to the row of the matrix V (and therefore |
---|
| 722 | // of the matrix A) with the number P_row[j], where P_row is an array |
---|
| 723 | // representing the permutation matrix P in row-like format. Thus, each |
---|
| 724 | // j-th linearly dependent column of the matrix U should be replaced by |
---|
| 725 | // a column of the unity matrix with the number P_row[j]. |
---|
| 726 | // |
---|
| 727 | // The code that repairs the matrix A may look like follows: |
---|
| 728 | // |
---|
| 729 | // for (j = rank+1; j <= n; j++) |
---|
| 730 | // { replace column Q_col[j] of the matrix A by column P_row[j] of |
---|
| 731 | // the unity matrix; |
---|
| 732 | // } |
---|
| 733 | // |
---|
| 734 | // where rank, P_row, and Q_col are members of the structure LUX. */ |
---|
| 735 | |
---|
| 736 | int lux_decomp(LUX *lux, int (*col)(void *info, int j, int ind[], |
---|
| 737 | mpq_t val[]), void *info) |
---|
| 738 | { int n = lux->n; |
---|
| 739 | LUXELM **V_row = lux->V_row; |
---|
| 740 | LUXELM **V_col = lux->V_col; |
---|
| 741 | int *P_row = lux->P_row; |
---|
| 742 | int *P_col = lux->P_col; |
---|
| 743 | int *Q_row = lux->Q_row; |
---|
| 744 | int *Q_col = lux->Q_col; |
---|
| 745 | LUXELM *piv, *vij; |
---|
| 746 | LUXWKA *wka; |
---|
| 747 | int i, j, k, p, q, t, *flag; |
---|
| 748 | mpq_t *work; |
---|
| 749 | /* allocate working area */ |
---|
| 750 | wka = xmalloc(sizeof(LUXWKA)); |
---|
| 751 | wka->R_len = xcalloc(1+n, sizeof(int)); |
---|
| 752 | wka->R_head = xcalloc(1+n, sizeof(int)); |
---|
| 753 | wka->R_prev = xcalloc(1+n, sizeof(int)); |
---|
| 754 | wka->R_next = xcalloc(1+n, sizeof(int)); |
---|
| 755 | wka->C_len = xcalloc(1+n, sizeof(int)); |
---|
| 756 | wka->C_head = xcalloc(1+n, sizeof(int)); |
---|
| 757 | wka->C_prev = xcalloc(1+n, sizeof(int)); |
---|
| 758 | wka->C_next = xcalloc(1+n, sizeof(int)); |
---|
| 759 | /* initialize LU-factorization data structures */ |
---|
| 760 | initialize(lux, col, info, wka); |
---|
| 761 | /* allocate working arrays */ |
---|
| 762 | flag = xcalloc(1+n, sizeof(int)); |
---|
| 763 | work = xcalloc(1+n, sizeof(mpq_t)); |
---|
| 764 | for (k = 1; k <= n; k++) |
---|
| 765 | { flag[k] = 0; |
---|
| 766 | mpq_init(work[k]); |
---|
| 767 | } |
---|
| 768 | /* main elimination loop */ |
---|
| 769 | for (k = 1; k <= n; k++) |
---|
| 770 | { /* choose a pivot element v[p,q] */ |
---|
| 771 | piv = find_pivot(lux, wka); |
---|
| 772 | if (piv == NULL) |
---|
| 773 | { /* no pivot can be chosen, because the active submatrix is |
---|
| 774 | empty */ |
---|
| 775 | break; |
---|
| 776 | } |
---|
| 777 | /* determine row and column indices of the pivot element */ |
---|
| 778 | p = piv->i, q = piv->j; |
---|
| 779 | /* let v[p,q] correspond to u[i',j']; permute k-th and i'-th |
---|
| 780 | rows and k-th and j'-th columns of the matrix U = P*V*Q to |
---|
| 781 | move the element u[i',j'] to the position u[k,k] */ |
---|
| 782 | i = P_col[p], j = Q_row[q]; |
---|
| 783 | xassert(k <= i && i <= n && k <= j && j <= n); |
---|
| 784 | /* permute k-th and i-th rows of the matrix U */ |
---|
| 785 | t = P_row[k]; |
---|
| 786 | P_row[i] = t, P_col[t] = i; |
---|
| 787 | P_row[k] = p, P_col[p] = k; |
---|
| 788 | /* permute k-th and j-th columns of the matrix U */ |
---|
| 789 | t = Q_col[k]; |
---|
| 790 | Q_col[j] = t, Q_row[t] = j; |
---|
| 791 | Q_col[k] = q, Q_row[q] = k; |
---|
| 792 | /* eliminate subdiagonal elements of k-th column of the matrix |
---|
| 793 | U = P*V*Q using the pivot element u[k,k] = v[p,q] */ |
---|
| 794 | eliminate(lux, wka, piv, flag, work); |
---|
| 795 | } |
---|
| 796 | /* determine the rank of A (and V) */ |
---|
| 797 | lux->rank = k - 1; |
---|
| 798 | /* free working arrays */ |
---|
| 799 | xfree(flag); |
---|
| 800 | for (k = 1; k <= n; k++) mpq_clear(work[k]); |
---|
| 801 | xfree(work); |
---|
| 802 | /* build column lists of the matrix V using its row lists */ |
---|
| 803 | for (j = 1; j <= n; j++) |
---|
| 804 | xassert(V_col[j] == NULL); |
---|
| 805 | for (i = 1; i <= n; i++) |
---|
| 806 | { for (vij = V_row[i]; vij != NULL; vij = vij->r_next) |
---|
| 807 | { j = vij->j; |
---|
| 808 | vij->c_prev = NULL; |
---|
| 809 | vij->c_next = V_col[j]; |
---|
| 810 | if (vij->c_next != NULL) vij->c_next->c_prev = vij; |
---|
| 811 | V_col[j] = vij; |
---|
| 812 | } |
---|
| 813 | } |
---|
| 814 | /* free working area */ |
---|
| 815 | xfree(wka->R_len); |
---|
| 816 | xfree(wka->R_head); |
---|
| 817 | xfree(wka->R_prev); |
---|
| 818 | xfree(wka->R_next); |
---|
| 819 | xfree(wka->C_len); |
---|
| 820 | xfree(wka->C_head); |
---|
| 821 | xfree(wka->C_prev); |
---|
| 822 | xfree(wka->C_next); |
---|
| 823 | xfree(wka); |
---|
| 824 | /* return to the calling program */ |
---|
| 825 | return (lux->rank < n); |
---|
| 826 | } |
---|
| 827 | |
---|
| 828 | /*---------------------------------------------------------------------- |
---|
| 829 | // lux_f_solve - solve system F*x = b or F'*x = b. |
---|
| 830 | // |
---|
| 831 | // SYNOPSIS |
---|
| 832 | // |
---|
| 833 | // #include "glplux.h" |
---|
| 834 | // void lux_f_solve(LUX *lux, int tr, mpq_t x[]); |
---|
| 835 | // |
---|
| 836 | // DESCRIPTION |
---|
| 837 | // |
---|
| 838 | // The routine lux_f_solve solves either the system F*x = b (if the |
---|
| 839 | // flag tr is zero) or the system F'*x = b (if the flag tr is non-zero), |
---|
| 840 | // where the matrix F is a component of LU-factorization specified by |
---|
| 841 | // the parameter lux, F' is a matrix transposed to F. |
---|
| 842 | // |
---|
| 843 | // On entry the array x should contain elements of the right-hand side |
---|
| 844 | // vector b in locations x[1], ..., x[n], where n is the order of the |
---|
| 845 | // matrix F. On exit this array will contain elements of the solution |
---|
| 846 | // vector x in the same locations. */ |
---|
| 847 | |
---|
| 848 | void lux_f_solve(LUX *lux, int tr, mpq_t x[]) |
---|
| 849 | { int n = lux->n; |
---|
| 850 | LUXELM **F_row = lux->F_row; |
---|
| 851 | LUXELM **F_col = lux->F_col; |
---|
| 852 | int *P_row = lux->P_row; |
---|
| 853 | LUXELM *fik, *fkj; |
---|
| 854 | int i, j, k; |
---|
| 855 | mpq_t temp; |
---|
| 856 | mpq_init(temp); |
---|
| 857 | if (!tr) |
---|
| 858 | { /* solve the system F*x = b */ |
---|
| 859 | for (j = 1; j <= n; j++) |
---|
| 860 | { k = P_row[j]; |
---|
| 861 | if (mpq_sgn(x[k]) != 0) |
---|
| 862 | { for (fik = F_col[k]; fik != NULL; fik = fik->c_next) |
---|
| 863 | { mpq_mul(temp, fik->val, x[k]); |
---|
| 864 | mpq_sub(x[fik->i], x[fik->i], temp); |
---|
| 865 | } |
---|
| 866 | } |
---|
| 867 | } |
---|
| 868 | } |
---|
| 869 | else |
---|
| 870 | { /* solve the system F'*x = b */ |
---|
| 871 | for (i = n; i >= 1; i--) |
---|
| 872 | { k = P_row[i]; |
---|
| 873 | if (mpq_sgn(x[k]) != 0) |
---|
| 874 | { for (fkj = F_row[k]; fkj != NULL; fkj = fkj->r_next) |
---|
| 875 | { mpq_mul(temp, fkj->val, x[k]); |
---|
| 876 | mpq_sub(x[fkj->j], x[fkj->j], temp); |
---|
| 877 | } |
---|
| 878 | } |
---|
| 879 | } |
---|
| 880 | } |
---|
| 881 | mpq_clear(temp); |
---|
| 882 | return; |
---|
| 883 | } |
---|
| 884 | |
---|
| 885 | /*---------------------------------------------------------------------- |
---|
| 886 | // lux_v_solve - solve system V*x = b or V'*x = b. |
---|
| 887 | // |
---|
| 888 | // SYNOPSIS |
---|
| 889 | // |
---|
| 890 | // #include "glplux.h" |
---|
| 891 | // void lux_v_solve(LUX *lux, int tr, double x[]); |
---|
| 892 | // |
---|
| 893 | // DESCRIPTION |
---|
| 894 | // |
---|
| 895 | // The routine lux_v_solve solves either the system V*x = b (if the |
---|
| 896 | // flag tr is zero) or the system V'*x = b (if the flag tr is non-zero), |
---|
| 897 | // where the matrix V is a component of LU-factorization specified by |
---|
| 898 | // the parameter lux, V' is a matrix transposed to V. |
---|
| 899 | // |
---|
| 900 | // On entry the array x should contain elements of the right-hand side |
---|
| 901 | // vector b in locations x[1], ..., x[n], where n is the order of the |
---|
| 902 | // matrix V. On exit this array will contain elements of the solution |
---|
| 903 | // vector x in the same locations. */ |
---|
| 904 | |
---|
| 905 | void lux_v_solve(LUX *lux, int tr, mpq_t x[]) |
---|
| 906 | { int n = lux->n; |
---|
| 907 | mpq_t *V_piv = lux->V_piv; |
---|
| 908 | LUXELM **V_row = lux->V_row; |
---|
| 909 | LUXELM **V_col = lux->V_col; |
---|
| 910 | int *P_row = lux->P_row; |
---|
| 911 | int *Q_col = lux->Q_col; |
---|
| 912 | LUXELM *vij; |
---|
| 913 | int i, j, k; |
---|
| 914 | mpq_t *b, temp; |
---|
| 915 | b = xcalloc(1+n, sizeof(mpq_t)); |
---|
| 916 | for (k = 1; k <= n; k++) |
---|
| 917 | mpq_init(b[k]), mpq_set(b[k], x[k]), mpq_set_si(x[k], 0, 1); |
---|
| 918 | mpq_init(temp); |
---|
| 919 | if (!tr) |
---|
| 920 | { /* solve the system V*x = b */ |
---|
| 921 | for (k = n; k >= 1; k--) |
---|
| 922 | { i = P_row[k], j = Q_col[k]; |
---|
| 923 | if (mpq_sgn(b[i]) != 0) |
---|
| 924 | { mpq_set(x[j], b[i]); |
---|
| 925 | mpq_div(x[j], x[j], V_piv[i]); |
---|
| 926 | for (vij = V_col[j]; vij != NULL; vij = vij->c_next) |
---|
| 927 | { mpq_mul(temp, vij->val, x[j]); |
---|
| 928 | mpq_sub(b[vij->i], b[vij->i], temp); |
---|
| 929 | } |
---|
| 930 | } |
---|
| 931 | } |
---|
| 932 | } |
---|
| 933 | else |
---|
| 934 | { /* solve the system V'*x = b */ |
---|
| 935 | for (k = 1; k <= n; k++) |
---|
| 936 | { i = P_row[k], j = Q_col[k]; |
---|
| 937 | if (mpq_sgn(b[j]) != 0) |
---|
| 938 | { mpq_set(x[i], b[j]); |
---|
| 939 | mpq_div(x[i], x[i], V_piv[i]); |
---|
| 940 | for (vij = V_row[i]; vij != NULL; vij = vij->r_next) |
---|
| 941 | { mpq_mul(temp, vij->val, x[i]); |
---|
| 942 | mpq_sub(b[vij->j], b[vij->j], temp); |
---|
| 943 | } |
---|
| 944 | } |
---|
| 945 | } |
---|
| 946 | } |
---|
| 947 | for (k = 1; k <= n; k++) mpq_clear(b[k]); |
---|
| 948 | mpq_clear(temp); |
---|
| 949 | xfree(b); |
---|
| 950 | return; |
---|
| 951 | } |
---|
| 952 | |
---|
| 953 | /*---------------------------------------------------------------------- |
---|
| 954 | // lux_solve - solve system A*x = b or A'*x = b. |
---|
| 955 | // |
---|
| 956 | // SYNOPSIS |
---|
| 957 | // |
---|
| 958 | // #include "glplux.h" |
---|
| 959 | // void lux_solve(LUX *lux, int tr, mpq_t x[]); |
---|
| 960 | // |
---|
| 961 | // DESCRIPTION |
---|
| 962 | // |
---|
| 963 | // The routine lux_solve solves either the system A*x = b (if the flag |
---|
| 964 | // tr is zero) or the system A'*x = b (if the flag tr is non-zero), |
---|
| 965 | // where the parameter lux specifies LU-factorization of the matrix A, |
---|
| 966 | // A' is a matrix transposed to A. |
---|
| 967 | // |
---|
| 968 | // On entry the array x should contain elements of the right-hand side |
---|
| 969 | // vector b in locations x[1], ..., x[n], where n is the order of the |
---|
| 970 | // matrix A. On exit this array will contain elements of the solution |
---|
| 971 | // vector x in the same locations. */ |
---|
| 972 | |
---|
| 973 | void lux_solve(LUX *lux, int tr, mpq_t x[]) |
---|
| 974 | { if (lux->rank < lux->n) |
---|
| 975 | xfault("lux_solve: LU-factorization has incomplete rank\n"); |
---|
| 976 | if (!tr) |
---|
| 977 | { /* A = F*V, therefore inv(A) = inv(V)*inv(F) */ |
---|
| 978 | lux_f_solve(lux, 0, x); |
---|
| 979 | lux_v_solve(lux, 0, x); |
---|
| 980 | } |
---|
| 981 | else |
---|
| 982 | { /* A' = V'*F', therefore inv(A') = inv(F')*inv(V') */ |
---|
| 983 | lux_v_solve(lux, 1, x); |
---|
| 984 | lux_f_solve(lux, 1, x); |
---|
| 985 | } |
---|
| 986 | return; |
---|
| 987 | } |
---|
| 988 | |
---|
| 989 | /*---------------------------------------------------------------------- |
---|
| 990 | // lux_delete - delete LU-factorization. |
---|
| 991 | // |
---|
| 992 | // SYNOPSIS |
---|
| 993 | // |
---|
| 994 | // #include "glplux.h" |
---|
| 995 | // void lux_delete(LUX *lux); |
---|
| 996 | // |
---|
| 997 | // DESCRIPTION |
---|
| 998 | // |
---|
| 999 | // The routine lux_delete deletes LU-factorization data structure, |
---|
| 1000 | // which the parameter lux points to, freeing all the memory allocated |
---|
| 1001 | // to this object. */ |
---|
| 1002 | |
---|
| 1003 | void lux_delete(LUX *lux) |
---|
| 1004 | { int n = lux->n; |
---|
| 1005 | LUXELM *fij, *vij; |
---|
| 1006 | int i; |
---|
| 1007 | for (i = 1; i <= n; i++) |
---|
| 1008 | { for (fij = lux->F_row[i]; fij != NULL; fij = fij->r_next) |
---|
| 1009 | mpq_clear(fij->val); |
---|
| 1010 | mpq_clear(lux->V_piv[i]); |
---|
| 1011 | for (vij = lux->V_row[i]; vij != NULL; vij = vij->r_next) |
---|
| 1012 | mpq_clear(vij->val); |
---|
| 1013 | } |
---|
| 1014 | dmp_delete_pool(lux->pool); |
---|
| 1015 | xfree(lux->F_row); |
---|
| 1016 | xfree(lux->F_col); |
---|
| 1017 | xfree(lux->V_piv); |
---|
| 1018 | xfree(lux->V_row); |
---|
| 1019 | xfree(lux->V_col); |
---|
| 1020 | xfree(lux->P_row); |
---|
| 1021 | xfree(lux->P_col); |
---|
| 1022 | xfree(lux->Q_row); |
---|
| 1023 | xfree(lux->Q_col); |
---|
| 1024 | xfree(lux); |
---|
| 1025 | return; |
---|
| 1026 | } |
---|
| 1027 | |
---|
| 1028 | /* eof */ |
---|