1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/deps/glpk/src/glplux.c	Sun Nov 06 20:59:10 2011 +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, 2011 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 */