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 */ 1029