1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/glplux.c Mon Dec 06 13:09:21 2010 +0100
1.3 @@ -0,0 +1,1028 @@
1.4 +/* glplux.c */
1.5 +
1.6 +/***********************************************************************
1.7 +* This code is part of GLPK (GNU Linear Programming Kit).
1.8 +*
1.9 +* Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
1.10 +* 2009, 2010 Andrew Makhorin, Department for Applied Informatics,
1.11 +* Moscow Aviation Institute, Moscow, Russia. All rights reserved.
1.12 +* E-mail: <mao@gnu.org>.
1.13 +*
1.14 +* GLPK is free software: you can redistribute it and/or modify it
1.15 +* under the terms of the GNU General Public License as published by
1.16 +* the Free Software Foundation, either version 3 of the License, or
1.17 +* (at your option) any later version.
1.18 +*
1.19 +* GLPK is distributed in the hope that it will be useful, but WITHOUT
1.20 +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
1.21 +* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
1.22 +* License for more details.
1.23 +*
1.24 +* You should have received a copy of the GNU General Public License
1.25 +* along with GLPK. If not, see <http://www.gnu.org/licenses/>.
1.26 +***********************************************************************/
1.27 +
1.28 +#include "glplux.h"
1.29 +#define xfault xerror
1.30 +#define dmp_create_poolx(size) dmp_create_pool()
1.31 +
1.32 +/*----------------------------------------------------------------------
1.33 +// lux_create - create LU-factorization.
1.34 +//
1.35 +// SYNOPSIS
1.36 +//
1.37 +// #include "glplux.h"
1.38 +// LUX *lux_create(int n);
1.39 +//
1.40 +// DESCRIPTION
1.41 +//
1.42 +// The routine lux_create creates LU-factorization data structure for
1.43 +// a matrix of the order n. Initially the factorization corresponds to
1.44 +// the unity matrix (F = V = P = Q = I, so A = I).
1.45 +//
1.46 +// RETURNS
1.47 +//
1.48 +// The routine returns a pointer to the created LU-factorization data
1.49 +// structure, which represents the unity matrix of the order n. */
1.50 +
1.51 +LUX *lux_create(int n)
1.52 +{ LUX *lux;
1.53 + int k;
1.54 + if (n < 1)
1.55 + xfault("lux_create: n = %d; invalid parameter\n", n);
1.56 + lux = xmalloc(sizeof(LUX));
1.57 + lux->n = n;
1.58 + lux->pool = dmp_create_poolx(sizeof(LUXELM));
1.59 + lux->F_row = xcalloc(1+n, sizeof(LUXELM *));
1.60 + lux->F_col = xcalloc(1+n, sizeof(LUXELM *));
1.61 + lux->V_piv = xcalloc(1+n, sizeof(mpq_t));
1.62 + lux->V_row = xcalloc(1+n, sizeof(LUXELM *));
1.63 + lux->V_col = xcalloc(1+n, sizeof(LUXELM *));
1.64 + lux->P_row = xcalloc(1+n, sizeof(int));
1.65 + lux->P_col = xcalloc(1+n, sizeof(int));
1.66 + lux->Q_row = xcalloc(1+n, sizeof(int));
1.67 + lux->Q_col = xcalloc(1+n, sizeof(int));
1.68 + for (k = 1; k <= n; k++)
1.69 + { lux->F_row[k] = lux->F_col[k] = NULL;
1.70 + mpq_init(lux->V_piv[k]);
1.71 + mpq_set_si(lux->V_piv[k], 1, 1);
1.72 + lux->V_row[k] = lux->V_col[k] = NULL;
1.73 + lux->P_row[k] = lux->P_col[k] = k;
1.74 + lux->Q_row[k] = lux->Q_col[k] = k;
1.75 + }
1.76 + lux->rank = n;
1.77 + return lux;
1.78 +}
1.79 +
1.80 +/*----------------------------------------------------------------------
1.81 +// initialize - initialize LU-factorization data structures.
1.82 +//
1.83 +// This routine initializes data structures for subsequent computing
1.84 +// the LU-factorization of a given matrix A, which is specified by the
1.85 +// formal routine col. On exit V = A and F = P = Q = I, where I is the
1.86 +// unity matrix. */
1.87 +
1.88 +static void initialize(LUX *lux, int (*col)(void *info, int j,
1.89 + int ind[], mpq_t val[]), void *info, LUXWKA *wka)
1.90 +{ int n = lux->n;
1.91 + DMP *pool = lux->pool;
1.92 + LUXELM **F_row = lux->F_row;
1.93 + LUXELM **F_col = lux->F_col;
1.94 + mpq_t *V_piv = lux->V_piv;
1.95 + LUXELM **V_row = lux->V_row;
1.96 + LUXELM **V_col = lux->V_col;
1.97 + int *P_row = lux->P_row;
1.98 + int *P_col = lux->P_col;
1.99 + int *Q_row = lux->Q_row;
1.100 + int *Q_col = lux->Q_col;
1.101 + int *R_len = wka->R_len;
1.102 + int *R_head = wka->R_head;
1.103 + int *R_prev = wka->R_prev;
1.104 + int *R_next = wka->R_next;
1.105 + int *C_len = wka->C_len;
1.106 + int *C_head = wka->C_head;
1.107 + int *C_prev = wka->C_prev;
1.108 + int *C_next = wka->C_next;
1.109 + LUXELM *fij, *vij;
1.110 + int i, j, k, len, *ind;
1.111 + mpq_t *val;
1.112 + /* F := I */
1.113 + for (i = 1; i <= n; i++)
1.114 + { while (F_row[i] != NULL)
1.115 + { fij = F_row[i], F_row[i] = fij->r_next;
1.116 + mpq_clear(fij->val);
1.117 + dmp_free_atom(pool, fij, sizeof(LUXELM));
1.118 + }
1.119 + }
1.120 + for (j = 1; j <= n; j++) F_col[j] = NULL;
1.121 + /* V := 0 */
1.122 + for (k = 1; k <= n; k++) mpq_set_si(V_piv[k], 0, 1);
1.123 + for (i = 1; i <= n; i++)
1.124 + { while (V_row[i] != NULL)
1.125 + { vij = V_row[i], V_row[i] = vij->r_next;
1.126 + mpq_clear(vij->val);
1.127 + dmp_free_atom(pool, vij, sizeof(LUXELM));
1.128 + }
1.129 + }
1.130 + for (j = 1; j <= n; j++) V_col[j] = NULL;
1.131 + /* V := A */
1.132 + ind = xcalloc(1+n, sizeof(int));
1.133 + val = xcalloc(1+n, sizeof(mpq_t));
1.134 + for (k = 1; k <= n; k++) mpq_init(val[k]);
1.135 + for (j = 1; j <= n; j++)
1.136 + { /* obtain j-th column of matrix A */
1.137 + len = col(info, j, ind, val);
1.138 + if (!(0 <= len && len <= n))
1.139 + xfault("lux_decomp: j = %d: len = %d; invalid column length"
1.140 + "\n", j, len);
1.141 + /* copy elements of j-th column to matrix V */
1.142 + for (k = 1; k <= len; k++)
1.143 + { /* get row index of a[i,j] */
1.144 + i = ind[k];
1.145 + if (!(1 <= i && i <= n))
1.146 + xfault("lux_decomp: j = %d: i = %d; row index out of ran"
1.147 + "ge\n", j, i);
1.148 + /* check for duplicate indices */
1.149 + if (V_row[i] != NULL && V_row[i]->j == j)
1.150 + xfault("lux_decomp: j = %d: i = %d; duplicate row indice"
1.151 + "s not allowed\n", j, i);
1.152 + /* check for zero value */
1.153 + if (mpq_sgn(val[k]) == 0)
1.154 + xfault("lux_decomp: j = %d: i = %d; zero elements not al"
1.155 + "lowed\n", j, i);
1.156 + /* add new element v[i,j] = a[i,j] to V */
1.157 + vij = dmp_get_atom(pool, sizeof(LUXELM));
1.158 + vij->i = i, vij->j = j;
1.159 + mpq_init(vij->val);
1.160 + mpq_set(vij->val, val[k]);
1.161 + vij->r_prev = NULL;
1.162 + vij->r_next = V_row[i];
1.163 + vij->c_prev = NULL;
1.164 + vij->c_next = V_col[j];
1.165 + if (vij->r_next != NULL) vij->r_next->r_prev = vij;
1.166 + if (vij->c_next != NULL) vij->c_next->c_prev = vij;
1.167 + V_row[i] = V_col[j] = vij;
1.168 + }
1.169 + }
1.170 + xfree(ind);
1.171 + for (k = 1; k <= n; k++) mpq_clear(val[k]);
1.172 + xfree(val);
1.173 + /* P := Q := I */
1.174 + for (k = 1; k <= n; k++)
1.175 + P_row[k] = P_col[k] = Q_row[k] = Q_col[k] = k;
1.176 + /* the rank of A and V is not determined yet */
1.177 + lux->rank = -1;
1.178 + /* initially the entire matrix V is active */
1.179 + /* determine its row lengths */
1.180 + for (i = 1; i <= n; i++)
1.181 + { len = 0;
1.182 + for (vij = V_row[i]; vij != NULL; vij = vij->r_next) len++;
1.183 + R_len[i] = len;
1.184 + }
1.185 + /* build linked lists of active rows */
1.186 + for (len = 0; len <= n; len++) R_head[len] = 0;
1.187 + for (i = 1; i <= n; i++)
1.188 + { len = R_len[i];
1.189 + R_prev[i] = 0;
1.190 + R_next[i] = R_head[len];
1.191 + if (R_next[i] != 0) R_prev[R_next[i]] = i;
1.192 + R_head[len] = i;
1.193 + }
1.194 + /* determine its column lengths */
1.195 + for (j = 1; j <= n; j++)
1.196 + { len = 0;
1.197 + for (vij = V_col[j]; vij != NULL; vij = vij->c_next) len++;
1.198 + C_len[j] = len;
1.199 + }
1.200 + /* build linked lists of active columns */
1.201 + for (len = 0; len <= n; len++) C_head[len] = 0;
1.202 + for (j = 1; j <= n; j++)
1.203 + { len = C_len[j];
1.204 + C_prev[j] = 0;
1.205 + C_next[j] = C_head[len];
1.206 + if (C_next[j] != 0) C_prev[C_next[j]] = j;
1.207 + C_head[len] = j;
1.208 + }
1.209 + return;
1.210 +}
1.211 +
1.212 +/*----------------------------------------------------------------------
1.213 +// find_pivot - choose a pivot element.
1.214 +//
1.215 +// This routine chooses a pivot element v[p,q] in the active submatrix
1.216 +// of matrix U = P*V*Q.
1.217 +//
1.218 +// It is assumed that on entry the matrix U has the following partially
1.219 +// triangularized form:
1.220 +//
1.221 +// 1 k n
1.222 +// 1 x x x x x x x x x x
1.223 +// . x x x x x x x x x
1.224 +// . . x x x x x x x x
1.225 +// . . . x x x x x x x
1.226 +// k . . . . * * * * * *
1.227 +// . . . . * * * * * *
1.228 +// . . . . * * * * * *
1.229 +// . . . . * * * * * *
1.230 +// . . . . * * * * * *
1.231 +// n . . . . * * * * * *
1.232 +//
1.233 +// where rows and columns k, k+1, ..., n belong to the active submatrix
1.234 +// (elements of the active submatrix are marked by '*').
1.235 +//
1.236 +// Since the matrix U = P*V*Q is not stored, the routine works with the
1.237 +// matrix V. It is assumed that the row-wise representation corresponds
1.238 +// to the matrix V, but the column-wise representation corresponds to
1.239 +// the active submatrix of the matrix V, i.e. elements of the matrix V,
1.240 +// which does not belong to the active submatrix, are missing from the
1.241 +// column linked lists. It is also assumed that each active row of the
1.242 +// matrix V is in the set R[len], where len is number of non-zeros in
1.243 +// the row, and each active column of the matrix V is in the set C[len],
1.244 +// where len is number of non-zeros in the column (in the latter case
1.245 +// only elements of the active submatrix are counted; such elements are
1.246 +// marked by '*' on the figure above).
1.247 +//
1.248 +// Due to exact arithmetic any non-zero element of the active submatrix
1.249 +// can be chosen as a pivot. However, to keep sparsity of the matrix V
1.250 +// the routine uses Markowitz strategy, trying to choose such element
1.251 +// v[p,q], which has smallest Markowitz cost (nr[p]-1) * (nc[q]-1),
1.252 +// where nr[p] and nc[q] are the number of non-zero elements, resp., in
1.253 +// p-th row and in q-th column of the active submatrix.
1.254 +//
1.255 +// In order to reduce the search, i.e. not to walk through all elements
1.256 +// of the active submatrix, the routine exploits a technique proposed by
1.257 +// I.Duff. This technique is based on using the sets R[len] and C[len]
1.258 +// of active rows and columns.
1.259 +//
1.260 +// On exit the routine returns a pointer to a pivot v[p,q] chosen, or
1.261 +// NULL, if the active submatrix is empty. */
1.262 +
1.263 +static LUXELM *find_pivot(LUX *lux, LUXWKA *wka)
1.264 +{ int n = lux->n;
1.265 + LUXELM **V_row = lux->V_row;
1.266 + LUXELM **V_col = lux->V_col;
1.267 + int *R_len = wka->R_len;
1.268 + int *R_head = wka->R_head;
1.269 + int *R_next = wka->R_next;
1.270 + int *C_len = wka->C_len;
1.271 + int *C_head = wka->C_head;
1.272 + int *C_next = wka->C_next;
1.273 + LUXELM *piv, *some, *vij;
1.274 + int i, j, len, min_len, ncand, piv_lim = 5;
1.275 + double best, cost;
1.276 + /* nothing is chosen so far */
1.277 + piv = NULL, best = DBL_MAX, ncand = 0;
1.278 + /* if in the active submatrix there is a column that has the only
1.279 + non-zero (column singleton), choose it as a pivot */
1.280 + j = C_head[1];
1.281 + if (j != 0)
1.282 + { xassert(C_len[j] == 1);
1.283 + piv = V_col[j];
1.284 + xassert(piv != NULL && piv->c_next == NULL);
1.285 + goto done;
1.286 + }
1.287 + /* if in the active submatrix there is a row that has the only
1.288 + non-zero (row singleton), choose it as a pivot */
1.289 + i = R_head[1];
1.290 + if (i != 0)
1.291 + { xassert(R_len[i] == 1);
1.292 + piv = V_row[i];
1.293 + xassert(piv != NULL && piv->r_next == NULL);
1.294 + goto done;
1.295 + }
1.296 + /* there are no singletons in the active submatrix; walk through
1.297 + other non-empty rows and columns */
1.298 + for (len = 2; len <= n; len++)
1.299 + { /* consider active columns having len non-zeros */
1.300 + for (j = C_head[len]; j != 0; j = C_next[j])
1.301 + { /* j-th column has len non-zeros */
1.302 + /* find an element in the row of minimal length */
1.303 + some = NULL, min_len = INT_MAX;
1.304 + for (vij = V_col[j]; vij != NULL; vij = vij->c_next)
1.305 + { if (min_len > R_len[vij->i])
1.306 + some = vij, min_len = R_len[vij->i];
1.307 + /* if Markowitz cost of this element is not greater than
1.308 + (len-1)**2, it can be chosen right now; this heuristic
1.309 + reduces the search and works well in many cases */
1.310 + if (min_len <= len)
1.311 + { piv = some;
1.312 + goto done;
1.313 + }
1.314 + }
1.315 + /* j-th column has been scanned */
1.316 + /* the minimal element found is a next pivot candidate */
1.317 + xassert(some != NULL);
1.318 + ncand++;
1.319 + /* compute its Markowitz cost */
1.320 + cost = (double)(min_len - 1) * (double)(len - 1);
1.321 + /* choose between the current candidate and this element */
1.322 + if (cost < best) piv = some, best = cost;
1.323 + /* if piv_lim candidates have been considered, there is a
1.324 + doubt that a much better candidate exists; therefore it
1.325 + is the time to terminate the search */
1.326 + if (ncand == piv_lim) goto done;
1.327 + }
1.328 + /* now consider active rows having len non-zeros */
1.329 + for (i = R_head[len]; i != 0; i = R_next[i])
1.330 + { /* i-th row has len non-zeros */
1.331 + /* find an element in the column of minimal length */
1.332 + some = NULL, min_len = INT_MAX;
1.333 + for (vij = V_row[i]; vij != NULL; vij = vij->r_next)
1.334 + { if (min_len > C_len[vij->j])
1.335 + some = vij, min_len = C_len[vij->j];
1.336 + /* if Markowitz cost of this element is not greater than
1.337 + (len-1)**2, it can be chosen right now; this heuristic
1.338 + reduces the search and works well in many cases */
1.339 + if (min_len <= len)
1.340 + { piv = some;
1.341 + goto done;
1.342 + }
1.343 + }
1.344 + /* i-th row has been scanned */
1.345 + /* the minimal element found is a next pivot candidate */
1.346 + xassert(some != NULL);
1.347 + ncand++;
1.348 + /* compute its Markowitz cost */
1.349 + cost = (double)(len - 1) * (double)(min_len - 1);
1.350 + /* choose between the current candidate and this element */
1.351 + if (cost < best) piv = some, best = cost;
1.352 + /* if piv_lim candidates have been considered, there is a
1.353 + doubt that a much better candidate exists; therefore it
1.354 + is the time to terminate the search */
1.355 + if (ncand == piv_lim) goto done;
1.356 + }
1.357 + }
1.358 +done: /* bring the pivot v[p,q] to the factorizing routine */
1.359 + return piv;
1.360 +}
1.361 +
1.362 +/*----------------------------------------------------------------------
1.363 +// eliminate - perform gaussian elimination.
1.364 +//
1.365 +// This routine performs elementary gaussian transformations in order
1.366 +// to eliminate subdiagonal elements in the k-th column of the matrix
1.367 +// U = P*V*Q using the pivot element u[k,k], where k is the number of
1.368 +// the current elimination step.
1.369 +//
1.370 +// The parameter piv specifies the pivot element v[p,q] = u[k,k].
1.371 +//
1.372 +// Each time when the routine applies the elementary transformation to
1.373 +// a non-pivot row of the matrix V, it stores the corresponding element
1.374 +// to the matrix F in order to keep the main equality A = F*V.
1.375 +//
1.376 +// The routine assumes that on entry the matrices L = P*F*inv(P) and
1.377 +// U = P*V*Q are the following:
1.378 +//
1.379 +// 1 k 1 k n
1.380 +// 1 1 . . . . . . . . . 1 x x x x x x x x x x
1.381 +// x 1 . . . . . . . . . x x x x x x x x x
1.382 +// x x 1 . . . . . . . . . x x x x x x x x
1.383 +// x x x 1 . . . . . . . . . x x x x x x x
1.384 +// k x x x x 1 . . . . . k . . . . * * * * * *
1.385 +// x x x x _ 1 . . . . . . . . # * * * * *
1.386 +// x x x x _ . 1 . . . . . . . # * * * * *
1.387 +// x x x x _ . . 1 . . . . . . # * * * * *
1.388 +// x x x x _ . . . 1 . . . . . # * * * * *
1.389 +// n x x x x _ . . . . 1 n . . . . # * * * * *
1.390 +//
1.391 +// matrix L matrix U
1.392 +//
1.393 +// where rows and columns of the matrix U with numbers k, k+1, ..., n
1.394 +// form the active submatrix (eliminated elements are marked by '#' and
1.395 +// other elements of the active submatrix are marked by '*'). Note that
1.396 +// each eliminated non-zero element u[i,k] of the matrix U gives the
1.397 +// corresponding element l[i,k] of the matrix L (marked by '_').
1.398 +//
1.399 +// Actually all operations are performed on the matrix V. Should note
1.400 +// that the row-wise representation corresponds to the matrix V, but the
1.401 +// column-wise representation corresponds to the active submatrix of the
1.402 +// matrix V, i.e. elements of the matrix V, which doesn't belong to the
1.403 +// active submatrix, are missing from the column linked lists.
1.404 +//
1.405 +// Let u[k,k] = v[p,q] be the pivot. In order to eliminate subdiagonal
1.406 +// elements u[i',k] = v[i,q], i' = k+1, k+2, ..., n, the routine applies
1.407 +// the following elementary gaussian transformations:
1.408 +//
1.409 +// (i-th row of V) := (i-th row of V) - f[i,p] * (p-th row of V),
1.410 +//
1.411 +// where f[i,p] = v[i,q] / v[p,q] is a gaussian multiplier.
1.412 +//
1.413 +// Additionally, in order to keep the main equality A = F*V, each time
1.414 +// when the routine applies the transformation to i-th row of the matrix
1.415 +// V, it also adds f[i,p] as a new element to the matrix F.
1.416 +//
1.417 +// IMPORTANT: On entry the working arrays flag and work should contain
1.418 +// zeros. This status is provided by the routine on exit. */
1.419 +
1.420 +static void eliminate(LUX *lux, LUXWKA *wka, LUXELM *piv, int flag[],
1.421 + mpq_t work[])
1.422 +{ DMP *pool = lux->pool;
1.423 + LUXELM **F_row = lux->F_row;
1.424 + LUXELM **F_col = lux->F_col;
1.425 + mpq_t *V_piv = lux->V_piv;
1.426 + LUXELM **V_row = lux->V_row;
1.427 + LUXELM **V_col = lux->V_col;
1.428 + int *R_len = wka->R_len;
1.429 + int *R_head = wka->R_head;
1.430 + int *R_prev = wka->R_prev;
1.431 + int *R_next = wka->R_next;
1.432 + int *C_len = wka->C_len;
1.433 + int *C_head = wka->C_head;
1.434 + int *C_prev = wka->C_prev;
1.435 + int *C_next = wka->C_next;
1.436 + LUXELM *fip, *vij, *vpj, *viq, *next;
1.437 + mpq_t temp;
1.438 + int i, j, p, q;
1.439 + mpq_init(temp);
1.440 + /* determine row and column indices of the pivot v[p,q] */
1.441 + xassert(piv != NULL);
1.442 + p = piv->i, q = piv->j;
1.443 + /* remove p-th (pivot) row from the active set; it will never
1.444 + return there */
1.445 + if (R_prev[p] == 0)
1.446 + R_head[R_len[p]] = R_next[p];
1.447 + else
1.448 + R_next[R_prev[p]] = R_next[p];
1.449 + if (R_next[p] == 0)
1.450 + ;
1.451 + else
1.452 + R_prev[R_next[p]] = R_prev[p];
1.453 + /* remove q-th (pivot) column from the active set; it will never
1.454 + return there */
1.455 + if (C_prev[q] == 0)
1.456 + C_head[C_len[q]] = C_next[q];
1.457 + else
1.458 + C_next[C_prev[q]] = C_next[q];
1.459 + if (C_next[q] == 0)
1.460 + ;
1.461 + else
1.462 + C_prev[C_next[q]] = C_prev[q];
1.463 + /* store the pivot value in a separate array */
1.464 + mpq_set(V_piv[p], piv->val);
1.465 + /* remove the pivot from p-th row */
1.466 + if (piv->r_prev == NULL)
1.467 + V_row[p] = piv->r_next;
1.468 + else
1.469 + piv->r_prev->r_next = piv->r_next;
1.470 + if (piv->r_next == NULL)
1.471 + ;
1.472 + else
1.473 + piv->r_next->r_prev = piv->r_prev;
1.474 + R_len[p]--;
1.475 + /* remove the pivot from q-th column */
1.476 + if (piv->c_prev == NULL)
1.477 + V_col[q] = piv->c_next;
1.478 + else
1.479 + piv->c_prev->c_next = piv->c_next;
1.480 + if (piv->c_next == NULL)
1.481 + ;
1.482 + else
1.483 + piv->c_next->c_prev = piv->c_prev;
1.484 + C_len[q]--;
1.485 + /* free the space occupied by the pivot */
1.486 + mpq_clear(piv->val);
1.487 + dmp_free_atom(pool, piv, sizeof(LUXELM));
1.488 + /* walk through p-th (pivot) row, which already does not contain
1.489 + the pivot v[p,q], and do the following... */
1.490 + for (vpj = V_row[p]; vpj != NULL; vpj = vpj->r_next)
1.491 + { /* get column index of v[p,j] */
1.492 + j = vpj->j;
1.493 + /* store v[p,j] in the working array */
1.494 + flag[j] = 1;
1.495 + mpq_set(work[j], vpj->val);
1.496 + /* remove j-th column from the active set; it will return there
1.497 + later with a new length */
1.498 + if (C_prev[j] == 0)
1.499 + C_head[C_len[j]] = C_next[j];
1.500 + else
1.501 + C_next[C_prev[j]] = C_next[j];
1.502 + if (C_next[j] == 0)
1.503 + ;
1.504 + else
1.505 + C_prev[C_next[j]] = C_prev[j];
1.506 + /* v[p,j] leaves the active submatrix, so remove it from j-th
1.507 + column; however, v[p,j] is kept in p-th row */
1.508 + if (vpj->c_prev == NULL)
1.509 + V_col[j] = vpj->c_next;
1.510 + else
1.511 + vpj->c_prev->c_next = vpj->c_next;
1.512 + if (vpj->c_next == NULL)
1.513 + ;
1.514 + else
1.515 + vpj->c_next->c_prev = vpj->c_prev;
1.516 + C_len[j]--;
1.517 + }
1.518 + /* now walk through q-th (pivot) column, which already does not
1.519 + contain the pivot v[p,q], and perform gaussian elimination */
1.520 + while (V_col[q] != NULL)
1.521 + { /* element v[i,q] has to be eliminated */
1.522 + viq = V_col[q];
1.523 + /* get row index of v[i,q] */
1.524 + i = viq->i;
1.525 + /* remove i-th row from the active set; later it will return
1.526 + there with a new length */
1.527 + if (R_prev[i] == 0)
1.528 + R_head[R_len[i]] = R_next[i];
1.529 + else
1.530 + R_next[R_prev[i]] = R_next[i];
1.531 + if (R_next[i] == 0)
1.532 + ;
1.533 + else
1.534 + R_prev[R_next[i]] = R_prev[i];
1.535 + /* compute gaussian multiplier f[i,p] = v[i,q] / v[p,q] and
1.536 + store it in the matrix F */
1.537 + fip = dmp_get_atom(pool, sizeof(LUXELM));
1.538 + fip->i = i, fip->j = p;
1.539 + mpq_init(fip->val);
1.540 + mpq_div(fip->val, viq->val, V_piv[p]);
1.541 + fip->r_prev = NULL;
1.542 + fip->r_next = F_row[i];
1.543 + fip->c_prev = NULL;
1.544 + fip->c_next = F_col[p];
1.545 + if (fip->r_next != NULL) fip->r_next->r_prev = fip;
1.546 + if (fip->c_next != NULL) fip->c_next->c_prev = fip;
1.547 + F_row[i] = F_col[p] = fip;
1.548 + /* v[i,q] has to be eliminated, so remove it from i-th row */
1.549 + if (viq->r_prev == NULL)
1.550 + V_row[i] = viq->r_next;
1.551 + else
1.552 + viq->r_prev->r_next = viq->r_next;
1.553 + if (viq->r_next == NULL)
1.554 + ;
1.555 + else
1.556 + viq->r_next->r_prev = viq->r_prev;
1.557 + R_len[i]--;
1.558 + /* and also from q-th column */
1.559 + V_col[q] = viq->c_next;
1.560 + C_len[q]--;
1.561 + /* free the space occupied by v[i,q] */
1.562 + mpq_clear(viq->val);
1.563 + dmp_free_atom(pool, viq, sizeof(LUXELM));
1.564 + /* perform gaussian transformation:
1.565 + (i-th row) := (i-th row) - f[i,p] * (p-th row)
1.566 + note that now p-th row, which is in the working array,
1.567 + does not contain the pivot v[p,q], and i-th row does not
1.568 + contain the element v[i,q] to be eliminated */
1.569 + /* walk through i-th row and transform existing non-zero
1.570 + elements */
1.571 + for (vij = V_row[i]; vij != NULL; vij = next)
1.572 + { next = vij->r_next;
1.573 + /* get column index of v[i,j] */
1.574 + j = vij->j;
1.575 + /* v[i,j] := v[i,j] - f[i,p] * v[p,j] */
1.576 + if (flag[j])
1.577 + { /* v[p,j] != 0 */
1.578 + flag[j] = 0;
1.579 + mpq_mul(temp, fip->val, work[j]);
1.580 + mpq_sub(vij->val, vij->val, temp);
1.581 + if (mpq_sgn(vij->val) == 0)
1.582 + { /* new v[i,j] is zero, so remove it from the active
1.583 + submatrix */
1.584 + /* remove v[i,j] from i-th row */
1.585 + if (vij->r_prev == NULL)
1.586 + V_row[i] = vij->r_next;
1.587 + else
1.588 + vij->r_prev->r_next = vij->r_next;
1.589 + if (vij->r_next == NULL)
1.590 + ;
1.591 + else
1.592 + vij->r_next->r_prev = vij->r_prev;
1.593 + R_len[i]--;
1.594 + /* remove v[i,j] from j-th column */
1.595 + if (vij->c_prev == NULL)
1.596 + V_col[j] = vij->c_next;
1.597 + else
1.598 + vij->c_prev->c_next = vij->c_next;
1.599 + if (vij->c_next == NULL)
1.600 + ;
1.601 + else
1.602 + vij->c_next->c_prev = vij->c_prev;
1.603 + C_len[j]--;
1.604 + /* free the space occupied by v[i,j] */
1.605 + mpq_clear(vij->val);
1.606 + dmp_free_atom(pool, vij, sizeof(LUXELM));
1.607 + }
1.608 + }
1.609 + }
1.610 + /* now flag is the pattern of the set v[p,*] \ v[i,*] */
1.611 + /* walk through p-th (pivot) row and create new elements in
1.612 + i-th row, which appear due to fill-in */
1.613 + for (vpj = V_row[p]; vpj != NULL; vpj = vpj->r_next)
1.614 + { j = vpj->j;
1.615 + if (flag[j])
1.616 + { /* create new non-zero v[i,j] = 0 - f[i,p] * v[p,j] and
1.617 + add it to i-th row and j-th column */
1.618 + vij = dmp_get_atom(pool, sizeof(LUXELM));
1.619 + vij->i = i, vij->j = j;
1.620 + mpq_init(vij->val);
1.621 + mpq_mul(vij->val, fip->val, work[j]);
1.622 + mpq_neg(vij->val, vij->val);
1.623 + vij->r_prev = NULL;
1.624 + vij->r_next = V_row[i];
1.625 + vij->c_prev = NULL;
1.626 + vij->c_next = V_col[j];
1.627 + if (vij->r_next != NULL) vij->r_next->r_prev = vij;
1.628 + if (vij->c_next != NULL) vij->c_next->c_prev = vij;
1.629 + V_row[i] = V_col[j] = vij;
1.630 + R_len[i]++, C_len[j]++;
1.631 + }
1.632 + else
1.633 + { /* there is no fill-in, because v[i,j] already exists in
1.634 + i-th row; restore the flag, which was reset before */
1.635 + flag[j] = 1;
1.636 + }
1.637 + }
1.638 + /* now i-th row has been completely transformed and can return
1.639 + to the active set with a new length */
1.640 + R_prev[i] = 0;
1.641 + R_next[i] = R_head[R_len[i]];
1.642 + if (R_next[i] != 0) R_prev[R_next[i]] = i;
1.643 + R_head[R_len[i]] = i;
1.644 + }
1.645 + /* at this point q-th (pivot) column must be empty */
1.646 + xassert(C_len[q] == 0);
1.647 + /* walk through p-th (pivot) row again and do the following... */
1.648 + for (vpj = V_row[p]; vpj != NULL; vpj = vpj->r_next)
1.649 + { /* get column index of v[p,j] */
1.650 + j = vpj->j;
1.651 + /* erase v[p,j] from the working array */
1.652 + flag[j] = 0;
1.653 + mpq_set_si(work[j], 0, 1);
1.654 + /* now j-th column has been completely transformed, so it can
1.655 + return to the active list with a new length */
1.656 + C_prev[j] = 0;
1.657 + C_next[j] = C_head[C_len[j]];
1.658 + if (C_next[j] != 0) C_prev[C_next[j]] = j;
1.659 + C_head[C_len[j]] = j;
1.660 + }
1.661 + mpq_clear(temp);
1.662 + /* return to the factorizing routine */
1.663 + return;
1.664 +}
1.665 +
1.666 +/*----------------------------------------------------------------------
1.667 +// lux_decomp - compute LU-factorization.
1.668 +//
1.669 +// SYNOPSIS
1.670 +//
1.671 +// #include "glplux.h"
1.672 +// int lux_decomp(LUX *lux, int (*col)(void *info, int j, int ind[],
1.673 +// mpq_t val[]), void *info);
1.674 +//
1.675 +// DESCRIPTION
1.676 +//
1.677 +// The routine lux_decomp computes LU-factorization of a given square
1.678 +// matrix A.
1.679 +//
1.680 +// The parameter lux specifies LU-factorization data structure built by
1.681 +// means of the routine lux_create.
1.682 +//
1.683 +// The formal routine col specifies the original matrix A. In order to
1.684 +// obtain j-th column of the matrix A the routine lux_decomp calls the
1.685 +// routine col with the parameter j (1 <= j <= n, where n is the order
1.686 +// of A). In response the routine col should store row indices and
1.687 +// numerical values of non-zero elements of j-th column of A to the
1.688 +// locations ind[1], ..., ind[len] and val[1], ..., val[len], resp.,
1.689 +// where len is the number of non-zeros in j-th column, which should be
1.690 +// returned on exit. Neiter zero nor duplicate elements are allowed.
1.691 +//
1.692 +// The parameter info is a transit pointer passed to the formal routine
1.693 +// col; it can be used for various purposes.
1.694 +//
1.695 +// RETURNS
1.696 +//
1.697 +// The routine lux_decomp returns the singularity flag. Zero flag means
1.698 +// that the original matrix A is non-singular while non-zero flag means
1.699 +// that A is (exactly!) singular.
1.700 +//
1.701 +// Note that LU-factorization is valid in both cases, however, in case
1.702 +// of singularity some rows of the matrix V (including pivot elements)
1.703 +// will be empty.
1.704 +//
1.705 +// REPAIRING SINGULAR MATRIX
1.706 +//
1.707 +// If the routine lux_decomp returns non-zero flag, it provides all
1.708 +// necessary information that can be used for "repairing" the matrix A,
1.709 +// where "repairing" means replacing linearly dependent columns of the
1.710 +// matrix A by appropriate columns of the unity matrix. This feature is
1.711 +// needed when the routine lux_decomp is used for reinverting the basis
1.712 +// matrix within the simplex method procedure.
1.713 +//
1.714 +// On exit linearly dependent columns of the matrix U have the numbers
1.715 +// rank+1, rank+2, ..., n, where rank is the exact rank of the matrix A
1.716 +// stored by the routine to the member lux->rank. The correspondence
1.717 +// between columns of A and U is the same as between columns of V and U.
1.718 +// Thus, linearly dependent columns of the matrix A have the numbers
1.719 +// Q_col[rank+1], Q_col[rank+2], ..., Q_col[n], where Q_col is an array
1.720 +// representing the permutation matrix Q in column-like format. It is
1.721 +// understood that each j-th linearly dependent column of the matrix U
1.722 +// should be replaced by the unity vector, where all elements are zero
1.723 +// except the unity diagonal element u[j,j]. On the other hand j-th row
1.724 +// of the matrix U corresponds to the row of the matrix V (and therefore
1.725 +// of the matrix A) with the number P_row[j], where P_row is an array
1.726 +// representing the permutation matrix P in row-like format. Thus, each
1.727 +// j-th linearly dependent column of the matrix U should be replaced by
1.728 +// a column of the unity matrix with the number P_row[j].
1.729 +//
1.730 +// The code that repairs the matrix A may look like follows:
1.731 +//
1.732 +// for (j = rank+1; j <= n; j++)
1.733 +// { replace column Q_col[j] of the matrix A by column P_row[j] of
1.734 +// the unity matrix;
1.735 +// }
1.736 +//
1.737 +// where rank, P_row, and Q_col are members of the structure LUX. */
1.738 +
1.739 +int lux_decomp(LUX *lux, int (*col)(void *info, int j, int ind[],
1.740 + mpq_t val[]), void *info)
1.741 +{ int n = lux->n;
1.742 + LUXELM **V_row = lux->V_row;
1.743 + LUXELM **V_col = lux->V_col;
1.744 + int *P_row = lux->P_row;
1.745 + int *P_col = lux->P_col;
1.746 + int *Q_row = lux->Q_row;
1.747 + int *Q_col = lux->Q_col;
1.748 + LUXELM *piv, *vij;
1.749 + LUXWKA *wka;
1.750 + int i, j, k, p, q, t, *flag;
1.751 + mpq_t *work;
1.752 + /* allocate working area */
1.753 + wka = xmalloc(sizeof(LUXWKA));
1.754 + wka->R_len = xcalloc(1+n, sizeof(int));
1.755 + wka->R_head = xcalloc(1+n, sizeof(int));
1.756 + wka->R_prev = xcalloc(1+n, sizeof(int));
1.757 + wka->R_next = xcalloc(1+n, sizeof(int));
1.758 + wka->C_len = xcalloc(1+n, sizeof(int));
1.759 + wka->C_head = xcalloc(1+n, sizeof(int));
1.760 + wka->C_prev = xcalloc(1+n, sizeof(int));
1.761 + wka->C_next = xcalloc(1+n, sizeof(int));
1.762 + /* initialize LU-factorization data structures */
1.763 + initialize(lux, col, info, wka);
1.764 + /* allocate working arrays */
1.765 + flag = xcalloc(1+n, sizeof(int));
1.766 + work = xcalloc(1+n, sizeof(mpq_t));
1.767 + for (k = 1; k <= n; k++)
1.768 + { flag[k] = 0;
1.769 + mpq_init(work[k]);
1.770 + }
1.771 + /* main elimination loop */
1.772 + for (k = 1; k <= n; k++)
1.773 + { /* choose a pivot element v[p,q] */
1.774 + piv = find_pivot(lux, wka);
1.775 + if (piv == NULL)
1.776 + { /* no pivot can be chosen, because the active submatrix is
1.777 + empty */
1.778 + break;
1.779 + }
1.780 + /* determine row and column indices of the pivot element */
1.781 + p = piv->i, q = piv->j;
1.782 + /* let v[p,q] correspond to u[i',j']; permute k-th and i'-th
1.783 + rows and k-th and j'-th columns of the matrix U = P*V*Q to
1.784 + move the element u[i',j'] to the position u[k,k] */
1.785 + i = P_col[p], j = Q_row[q];
1.786 + xassert(k <= i && i <= n && k <= j && j <= n);
1.787 + /* permute k-th and i-th rows of the matrix U */
1.788 + t = P_row[k];
1.789 + P_row[i] = t, P_col[t] = i;
1.790 + P_row[k] = p, P_col[p] = k;
1.791 + /* permute k-th and j-th columns of the matrix U */
1.792 + t = Q_col[k];
1.793 + Q_col[j] = t, Q_row[t] = j;
1.794 + Q_col[k] = q, Q_row[q] = k;
1.795 + /* eliminate subdiagonal elements of k-th column of the matrix
1.796 + U = P*V*Q using the pivot element u[k,k] = v[p,q] */
1.797 + eliminate(lux, wka, piv, flag, work);
1.798 + }
1.799 + /* determine the rank of A (and V) */
1.800 + lux->rank = k - 1;
1.801 + /* free working arrays */
1.802 + xfree(flag);
1.803 + for (k = 1; k <= n; k++) mpq_clear(work[k]);
1.804 + xfree(work);
1.805 + /* build column lists of the matrix V using its row lists */
1.806 + for (j = 1; j <= n; j++)
1.807 + xassert(V_col[j] == NULL);
1.808 + for (i = 1; i <= n; i++)
1.809 + { for (vij = V_row[i]; vij != NULL; vij = vij->r_next)
1.810 + { j = vij->j;
1.811 + vij->c_prev = NULL;
1.812 + vij->c_next = V_col[j];
1.813 + if (vij->c_next != NULL) vij->c_next->c_prev = vij;
1.814 + V_col[j] = vij;
1.815 + }
1.816 + }
1.817 + /* free working area */
1.818 + xfree(wka->R_len);
1.819 + xfree(wka->R_head);
1.820 + xfree(wka->R_prev);
1.821 + xfree(wka->R_next);
1.822 + xfree(wka->C_len);
1.823 + xfree(wka->C_head);
1.824 + xfree(wka->C_prev);
1.825 + xfree(wka->C_next);
1.826 + xfree(wka);
1.827 + /* return to the calling program */
1.828 + return (lux->rank < n);
1.829 +}
1.830 +
1.831 +/*----------------------------------------------------------------------
1.832 +// lux_f_solve - solve system F*x = b or F'*x = b.
1.833 +//
1.834 +// SYNOPSIS
1.835 +//
1.836 +// #include "glplux.h"
1.837 +// void lux_f_solve(LUX *lux, int tr, mpq_t x[]);
1.838 +//
1.839 +// DESCRIPTION
1.840 +//
1.841 +// The routine lux_f_solve solves either the system F*x = b (if the
1.842 +// flag tr is zero) or the system F'*x = b (if the flag tr is non-zero),
1.843 +// where the matrix F is a component of LU-factorization specified by
1.844 +// the parameter lux, F' is a matrix transposed to F.
1.845 +//
1.846 +// On entry the array x should contain elements of the right-hand side
1.847 +// vector b in locations x[1], ..., x[n], where n is the order of the
1.848 +// matrix F. On exit this array will contain elements of the solution
1.849 +// vector x in the same locations. */
1.850 +
1.851 +void lux_f_solve(LUX *lux, int tr, mpq_t x[])
1.852 +{ int n = lux->n;
1.853 + LUXELM **F_row = lux->F_row;
1.854 + LUXELM **F_col = lux->F_col;
1.855 + int *P_row = lux->P_row;
1.856 + LUXELM *fik, *fkj;
1.857 + int i, j, k;
1.858 + mpq_t temp;
1.859 + mpq_init(temp);
1.860 + if (!tr)
1.861 + { /* solve the system F*x = b */
1.862 + for (j = 1; j <= n; j++)
1.863 + { k = P_row[j];
1.864 + if (mpq_sgn(x[k]) != 0)
1.865 + { for (fik = F_col[k]; fik != NULL; fik = fik->c_next)
1.866 + { mpq_mul(temp, fik->val, x[k]);
1.867 + mpq_sub(x[fik->i], x[fik->i], temp);
1.868 + }
1.869 + }
1.870 + }
1.871 + }
1.872 + else
1.873 + { /* solve the system F'*x = b */
1.874 + for (i = n; i >= 1; i--)
1.875 + { k = P_row[i];
1.876 + if (mpq_sgn(x[k]) != 0)
1.877 + { for (fkj = F_row[k]; fkj != NULL; fkj = fkj->r_next)
1.878 + { mpq_mul(temp, fkj->val, x[k]);
1.879 + mpq_sub(x[fkj->j], x[fkj->j], temp);
1.880 + }
1.881 + }
1.882 + }
1.883 + }
1.884 + mpq_clear(temp);
1.885 + return;
1.886 +}
1.887 +
1.888 +/*----------------------------------------------------------------------
1.889 +// lux_v_solve - solve system V*x = b or V'*x = b.
1.890 +//
1.891 +// SYNOPSIS
1.892 +//
1.893 +// #include "glplux.h"
1.894 +// void lux_v_solve(LUX *lux, int tr, double x[]);
1.895 +//
1.896 +// DESCRIPTION
1.897 +//
1.898 +// The routine lux_v_solve solves either the system V*x = b (if the
1.899 +// flag tr is zero) or the system V'*x = b (if the flag tr is non-zero),
1.900 +// where the matrix V is a component of LU-factorization specified by
1.901 +// the parameter lux, V' is a matrix transposed to V.
1.902 +//
1.903 +// On entry the array x should contain elements of the right-hand side
1.904 +// vector b in locations x[1], ..., x[n], where n is the order of the
1.905 +// matrix V. On exit this array will contain elements of the solution
1.906 +// vector x in the same locations. */
1.907 +
1.908 +void lux_v_solve(LUX *lux, int tr, mpq_t x[])
1.909 +{ int n = lux->n;
1.910 + mpq_t *V_piv = lux->V_piv;
1.911 + LUXELM **V_row = lux->V_row;
1.912 + LUXELM **V_col = lux->V_col;
1.913 + int *P_row = lux->P_row;
1.914 + int *Q_col = lux->Q_col;
1.915 + LUXELM *vij;
1.916 + int i, j, k;
1.917 + mpq_t *b, temp;
1.918 + b = xcalloc(1+n, sizeof(mpq_t));
1.919 + for (k = 1; k <= n; k++)
1.920 + mpq_init(b[k]), mpq_set(b[k], x[k]), mpq_set_si(x[k], 0, 1);
1.921 + mpq_init(temp);
1.922 + if (!tr)
1.923 + { /* solve the system V*x = b */
1.924 + for (k = n; k >= 1; k--)
1.925 + { i = P_row[k], j = Q_col[k];
1.926 + if (mpq_sgn(b[i]) != 0)
1.927 + { mpq_set(x[j], b[i]);
1.928 + mpq_div(x[j], x[j], V_piv[i]);
1.929 + for (vij = V_col[j]; vij != NULL; vij = vij->c_next)
1.930 + { mpq_mul(temp, vij->val, x[j]);
1.931 + mpq_sub(b[vij->i], b[vij->i], temp);
1.932 + }
1.933 + }
1.934 + }
1.935 + }
1.936 + else
1.937 + { /* solve the system V'*x = b */
1.938 + for (k = 1; k <= n; k++)
1.939 + { i = P_row[k], j = Q_col[k];
1.940 + if (mpq_sgn(b[j]) != 0)
1.941 + { mpq_set(x[i], b[j]);
1.942 + mpq_div(x[i], x[i], V_piv[i]);
1.943 + for (vij = V_row[i]; vij != NULL; vij = vij->r_next)
1.944 + { mpq_mul(temp, vij->val, x[i]);
1.945 + mpq_sub(b[vij->j], b[vij->j], temp);
1.946 + }
1.947 + }
1.948 + }
1.949 + }
1.950 + for (k = 1; k <= n; k++) mpq_clear(b[k]);
1.951 + mpq_clear(temp);
1.952 + xfree(b);
1.953 + return;
1.954 +}
1.955 +
1.956 +/*----------------------------------------------------------------------
1.957 +// lux_solve - solve system A*x = b or A'*x = b.
1.958 +//
1.959 +// SYNOPSIS
1.960 +//
1.961 +// #include "glplux.h"
1.962 +// void lux_solve(LUX *lux, int tr, mpq_t x[]);
1.963 +//
1.964 +// DESCRIPTION
1.965 +//
1.966 +// The routine lux_solve solves either the system A*x = b (if the flag
1.967 +// tr is zero) or the system A'*x = b (if the flag tr is non-zero),
1.968 +// where the parameter lux specifies LU-factorization of the matrix A,
1.969 +// A' is a matrix transposed to A.
1.970 +//
1.971 +// On entry the array x should contain elements of the right-hand side
1.972 +// vector b in locations x[1], ..., x[n], where n is the order of the
1.973 +// matrix A. On exit this array will contain elements of the solution
1.974 +// vector x in the same locations. */
1.975 +
1.976 +void lux_solve(LUX *lux, int tr, mpq_t x[])
1.977 +{ if (lux->rank < lux->n)
1.978 + xfault("lux_solve: LU-factorization has incomplete rank\n");
1.979 + if (!tr)
1.980 + { /* A = F*V, therefore inv(A) = inv(V)*inv(F) */
1.981 + lux_f_solve(lux, 0, x);
1.982 + lux_v_solve(lux, 0, x);
1.983 + }
1.984 + else
1.985 + { /* A' = V'*F', therefore inv(A') = inv(F')*inv(V') */
1.986 + lux_v_solve(lux, 1, x);
1.987 + lux_f_solve(lux, 1, x);
1.988 + }
1.989 + return;
1.990 +}
1.991 +
1.992 +/*----------------------------------------------------------------------
1.993 +// lux_delete - delete LU-factorization.
1.994 +//
1.995 +// SYNOPSIS
1.996 +//
1.997 +// #include "glplux.h"
1.998 +// void lux_delete(LUX *lux);
1.999 +//
1.1000 +// DESCRIPTION
1.1001 +//
1.1002 +// The routine lux_delete deletes LU-factorization data structure,
1.1003 +// which the parameter lux points to, freeing all the memory allocated
1.1004 +// to this object. */
1.1005 +
1.1006 +void lux_delete(LUX *lux)
1.1007 +{ int n = lux->n;
1.1008 + LUXELM *fij, *vij;
1.1009 + int i;
1.1010 + for (i = 1; i <= n; i++)
1.1011 + { for (fij = lux->F_row[i]; fij != NULL; fij = fij->r_next)
1.1012 + mpq_clear(fij->val);
1.1013 + mpq_clear(lux->V_piv[i]);
1.1014 + for (vij = lux->V_row[i]; vij != NULL; vij = vij->r_next)
1.1015 + mpq_clear(vij->val);
1.1016 + }
1.1017 + dmp_delete_pool(lux->pool);
1.1018 + xfree(lux->F_row);
1.1019 + xfree(lux->F_col);
1.1020 + xfree(lux->V_piv);
1.1021 + xfree(lux->V_row);
1.1022 + xfree(lux->V_col);
1.1023 + xfree(lux->P_row);
1.1024 + xfree(lux->P_col);
1.1025 + xfree(lux->Q_row);
1.1026 + xfree(lux->Q_col);
1.1027 + xfree(lux);
1.1028 + return;
1.1029 +}
1.1030 +
1.1031 +/* eof */