lemon-project-template-glpk

annotate deps/glpk/src/glplux.c @ 9:33de93886c88

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