alpar@9: /* glplux.c */
alpar@9: 
alpar@9: /***********************************************************************
alpar@9: *  This code is part of GLPK (GNU Linear Programming Kit).
alpar@9: *
alpar@9: *  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
alpar@9: *  2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics,
alpar@9: *  Moscow Aviation Institute, Moscow, Russia. All rights reserved.
alpar@9: *  E-mail: <mao@gnu.org>.
alpar@9: *
alpar@9: *  GLPK is free software: you can redistribute it and/or modify it
alpar@9: *  under the terms of the GNU General Public License as published by
alpar@9: *  the Free Software Foundation, either version 3 of the License, or
alpar@9: *  (at your option) any later version.
alpar@9: *
alpar@9: *  GLPK is distributed in the hope that it will be useful, but WITHOUT
alpar@9: *  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
alpar@9: *  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
alpar@9: *  License for more details.
alpar@9: *
alpar@9: *  You should have received a copy of the GNU General Public License
alpar@9: *  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
alpar@9: ***********************************************************************/
alpar@9: 
alpar@9: #include "glplux.h"
alpar@9: #define xfault xerror
alpar@9: #define dmp_create_poolx(size) dmp_create_pool()
alpar@9: 
alpar@9: /*----------------------------------------------------------------------
alpar@9: // lux_create - create LU-factorization.
alpar@9: //
alpar@9: // SYNOPSIS
alpar@9: //
alpar@9: // #include "glplux.h"
alpar@9: // LUX *lux_create(int n);
alpar@9: //
alpar@9: // DESCRIPTION
alpar@9: //
alpar@9: // The routine lux_create creates LU-factorization data structure for
alpar@9: // a matrix of the order n. Initially the factorization corresponds to
alpar@9: // the unity matrix (F = V = P = Q = I, so A = I).
alpar@9: //
alpar@9: // RETURNS
alpar@9: //
alpar@9: // The routine returns a pointer to the created LU-factorization data
alpar@9: // structure, which represents the unity matrix of the order n. */
alpar@9: 
alpar@9: LUX *lux_create(int n)
alpar@9: {     LUX *lux;
alpar@9:       int k;
alpar@9:       if (n < 1)
alpar@9:          xfault("lux_create: n = %d; invalid parameter\n", n);
alpar@9:       lux = xmalloc(sizeof(LUX));
alpar@9:       lux->n = n;
alpar@9:       lux->pool = dmp_create_poolx(sizeof(LUXELM));
alpar@9:       lux->F_row = xcalloc(1+n, sizeof(LUXELM *));
alpar@9:       lux->F_col = xcalloc(1+n, sizeof(LUXELM *));
alpar@9:       lux->V_piv = xcalloc(1+n, sizeof(mpq_t));
alpar@9:       lux->V_row = xcalloc(1+n, sizeof(LUXELM *));
alpar@9:       lux->V_col = xcalloc(1+n, sizeof(LUXELM *));
alpar@9:       lux->P_row = xcalloc(1+n, sizeof(int));
alpar@9:       lux->P_col = xcalloc(1+n, sizeof(int));
alpar@9:       lux->Q_row = xcalloc(1+n, sizeof(int));
alpar@9:       lux->Q_col = xcalloc(1+n, sizeof(int));
alpar@9:       for (k = 1; k <= n; k++)
alpar@9:       {  lux->F_row[k] = lux->F_col[k] = NULL;
alpar@9:          mpq_init(lux->V_piv[k]);
alpar@9:          mpq_set_si(lux->V_piv[k], 1, 1);
alpar@9:          lux->V_row[k] = lux->V_col[k] = NULL;
alpar@9:          lux->P_row[k] = lux->P_col[k] = k;
alpar@9:          lux->Q_row[k] = lux->Q_col[k] = k;
alpar@9:       }
alpar@9:       lux->rank = n;
alpar@9:       return lux;
alpar@9: }
alpar@9: 
alpar@9: /*----------------------------------------------------------------------
alpar@9: // initialize - initialize LU-factorization data structures.
alpar@9: //
alpar@9: // This routine initializes data structures for subsequent computing
alpar@9: // the LU-factorization of a given matrix A, which is specified by the
alpar@9: // formal routine col. On exit V = A and F = P = Q = I, where I is the
alpar@9: // unity matrix. */
alpar@9: 
alpar@9: static void initialize(LUX *lux, int (*col)(void *info, int j,
alpar@9:       int ind[], mpq_t val[]), void *info, LUXWKA *wka)
alpar@9: {     int n = lux->n;
alpar@9:       DMP *pool = lux->pool;
alpar@9:       LUXELM **F_row = lux->F_row;
alpar@9:       LUXELM **F_col = lux->F_col;
alpar@9:       mpq_t *V_piv = lux->V_piv;
alpar@9:       LUXELM **V_row = lux->V_row;
alpar@9:       LUXELM **V_col = lux->V_col;
alpar@9:       int *P_row = lux->P_row;
alpar@9:       int *P_col = lux->P_col;
alpar@9:       int *Q_row = lux->Q_row;
alpar@9:       int *Q_col = lux->Q_col;
alpar@9:       int *R_len = wka->R_len;
alpar@9:       int *R_head = wka->R_head;
alpar@9:       int *R_prev = wka->R_prev;
alpar@9:       int *R_next = wka->R_next;
alpar@9:       int *C_len = wka->C_len;
alpar@9:       int *C_head = wka->C_head;
alpar@9:       int *C_prev = wka->C_prev;
alpar@9:       int *C_next = wka->C_next;
alpar@9:       LUXELM *fij, *vij;
alpar@9:       int i, j, k, len, *ind;
alpar@9:       mpq_t *val;
alpar@9:       /* F := I */
alpar@9:       for (i = 1; i <= n; i++)
alpar@9:       {  while (F_row[i] != NULL)
alpar@9:          {  fij = F_row[i], F_row[i] = fij->r_next;
alpar@9:             mpq_clear(fij->val);
alpar@9:             dmp_free_atom(pool, fij, sizeof(LUXELM));
alpar@9:          }
alpar@9:       }
alpar@9:       for (j = 1; j <= n; j++) F_col[j] = NULL;
alpar@9:       /* V := 0 */
alpar@9:       for (k = 1; k <= n; k++) mpq_set_si(V_piv[k], 0, 1);
alpar@9:       for (i = 1; i <= n; i++)
alpar@9:       {  while (V_row[i] != NULL)
alpar@9:          {  vij = V_row[i], V_row[i] = vij->r_next;
alpar@9:             mpq_clear(vij->val);
alpar@9:             dmp_free_atom(pool, vij, sizeof(LUXELM));
alpar@9:          }
alpar@9:       }
alpar@9:       for (j = 1; j <= n; j++) V_col[j] = NULL;
alpar@9:       /* V := A */
alpar@9:       ind = xcalloc(1+n, sizeof(int));
alpar@9:       val = xcalloc(1+n, sizeof(mpq_t));
alpar@9:       for (k = 1; k <= n; k++) mpq_init(val[k]);
alpar@9:       for (j = 1; j <= n; j++)
alpar@9:       {  /* obtain j-th column of matrix A */
alpar@9:          len = col(info, j, ind, val);
alpar@9:          if (!(0 <= len && len <= n))
alpar@9:             xfault("lux_decomp: j = %d: len = %d; invalid column length"
alpar@9:                "\n", j, len);
alpar@9:          /* copy elements of j-th column to matrix V */
alpar@9:          for (k = 1; k <= len; k++)
alpar@9:          {  /* get row index of a[i,j] */
alpar@9:             i = ind[k];
alpar@9:             if (!(1 <= i && i <= n))
alpar@9:                xfault("lux_decomp: j = %d: i = %d; row index out of ran"
alpar@9:                   "ge\n", j, i);
alpar@9:             /* check for duplicate indices */
alpar@9:             if (V_row[i] != NULL && V_row[i]->j == j)
alpar@9:                xfault("lux_decomp: j = %d: i = %d; duplicate row indice"
alpar@9:                   "s not allowed\n", j, i);
alpar@9:             /* check for zero value */
alpar@9:             if (mpq_sgn(val[k]) == 0)
alpar@9:                xfault("lux_decomp: j = %d: i = %d; zero elements not al"
alpar@9:                   "lowed\n", j, i);
alpar@9:             /* add new element v[i,j] = a[i,j] to V */
alpar@9:             vij = dmp_get_atom(pool, sizeof(LUXELM));
alpar@9:             vij->i = i, vij->j = j;
alpar@9:             mpq_init(vij->val);
alpar@9:             mpq_set(vij->val, val[k]);
alpar@9:             vij->r_prev = NULL;
alpar@9:             vij->r_next = V_row[i];
alpar@9:             vij->c_prev = NULL;
alpar@9:             vij->c_next = V_col[j];
alpar@9:             if (vij->r_next != NULL) vij->r_next->r_prev = vij;
alpar@9:             if (vij->c_next != NULL) vij->c_next->c_prev = vij;
alpar@9:             V_row[i] = V_col[j] = vij;
alpar@9:          }
alpar@9:       }
alpar@9:       xfree(ind);
alpar@9:       for (k = 1; k <= n; k++) mpq_clear(val[k]);
alpar@9:       xfree(val);
alpar@9:       /* P := Q := I */
alpar@9:       for (k = 1; k <= n; k++)
alpar@9:          P_row[k] = P_col[k] = Q_row[k] = Q_col[k] = k;
alpar@9:       /* the rank of A and V is not determined yet */
alpar@9:       lux->rank = -1;
alpar@9:       /* initially the entire matrix V is active */
alpar@9:       /* determine its row lengths */
alpar@9:       for (i = 1; i <= n; i++)
alpar@9:       {  len = 0;
alpar@9:          for (vij = V_row[i]; vij != NULL; vij = vij->r_next) len++;
alpar@9:          R_len[i] = len;
alpar@9:       }
alpar@9:       /* build linked lists of active rows */
alpar@9:       for (len = 0; len <= n; len++) R_head[len] = 0;
alpar@9:       for (i = 1; i <= n; i++)
alpar@9:       {  len = R_len[i];
alpar@9:          R_prev[i] = 0;
alpar@9:          R_next[i] = R_head[len];
alpar@9:          if (R_next[i] != 0) R_prev[R_next[i]] = i;
alpar@9:          R_head[len] = i;
alpar@9:       }
alpar@9:       /* determine its column lengths */
alpar@9:       for (j = 1; j <= n; j++)
alpar@9:       {  len = 0;
alpar@9:          for (vij = V_col[j]; vij != NULL; vij = vij->c_next) len++;
alpar@9:          C_len[j] = len;
alpar@9:       }
alpar@9:       /* build linked lists of active columns */
alpar@9:       for (len = 0; len <= n; len++) C_head[len] = 0;
alpar@9:       for (j = 1; j <= n; j++)
alpar@9:       {  len = C_len[j];
alpar@9:          C_prev[j] = 0;
alpar@9:          C_next[j] = C_head[len];
alpar@9:          if (C_next[j] != 0) C_prev[C_next[j]] = j;
alpar@9:          C_head[len] = j;
alpar@9:       }
alpar@9:       return;
alpar@9: }
alpar@9: 
alpar@9: /*----------------------------------------------------------------------
alpar@9: // find_pivot - choose a pivot element.
alpar@9: //
alpar@9: // This routine chooses a pivot element v[p,q] in the active submatrix
alpar@9: // of matrix U = P*V*Q.
alpar@9: //
alpar@9: // It is assumed that on entry the matrix U has the following partially
alpar@9: // triangularized form:
alpar@9: //
alpar@9: //       1       k         n
alpar@9: //    1  x x x x x x x x x x
alpar@9: //       . x x x x x x x x x
alpar@9: //       . . x x x x x x x x
alpar@9: //       . . . x x x x x x x
alpar@9: //    k  . . . . * * * * * *
alpar@9: //       . . . . * * * * * *
alpar@9: //       . . . . * * * * * *
alpar@9: //       . . . . * * * * * *
alpar@9: //       . . . . * * * * * *
alpar@9: //    n  . . . . * * * * * *
alpar@9: //
alpar@9: // where rows and columns k, k+1, ..., n belong to the active submatrix
alpar@9: // (elements of the active submatrix are marked by '*').
alpar@9: //
alpar@9: // Since the matrix U = P*V*Q is not stored, the routine works with the
alpar@9: // matrix V. It is assumed that the row-wise representation corresponds
alpar@9: // to the matrix V, but the column-wise representation corresponds to
alpar@9: // the active submatrix of the matrix V, i.e. elements of the matrix V,
alpar@9: // which does not belong to the active submatrix, are missing from the
alpar@9: // column linked lists. It is also assumed that each active row of the
alpar@9: // matrix V is in the set R[len], where len is number of non-zeros in
alpar@9: // the row, and each active column of the matrix V is in the set C[len],
alpar@9: // where len is number of non-zeros in the column (in the latter case
alpar@9: // only elements of the active submatrix are counted; such elements are
alpar@9: // marked by '*' on the figure above).
alpar@9: //
alpar@9: // Due to exact arithmetic any non-zero element of the active submatrix
alpar@9: // can be chosen as a pivot. However, to keep sparsity of the matrix V
alpar@9: // the routine uses Markowitz strategy, trying to choose such element
alpar@9: // v[p,q], which has smallest Markowitz cost (nr[p]-1) * (nc[q]-1),
alpar@9: // where nr[p] and nc[q] are the number of non-zero elements, resp., in
alpar@9: // p-th row and in q-th column of the active submatrix.
alpar@9: //
alpar@9: // In order to reduce the search, i.e. not to walk through all elements
alpar@9: // of the active submatrix, the routine exploits a technique proposed by
alpar@9: // I.Duff. This technique is based on using the sets R[len] and C[len]
alpar@9: // of active rows and columns.
alpar@9: //
alpar@9: // On exit the routine returns a pointer to a pivot v[p,q] chosen, or
alpar@9: // NULL, if the active submatrix is empty. */
alpar@9: 
alpar@9: static LUXELM *find_pivot(LUX *lux, LUXWKA *wka)
alpar@9: {     int n = lux->n;
alpar@9:       LUXELM **V_row = lux->V_row;
alpar@9:       LUXELM **V_col = lux->V_col;
alpar@9:       int *R_len = wka->R_len;
alpar@9:       int *R_head = wka->R_head;
alpar@9:       int *R_next = wka->R_next;
alpar@9:       int *C_len = wka->C_len;
alpar@9:       int *C_head = wka->C_head;
alpar@9:       int *C_next = wka->C_next;
alpar@9:       LUXELM *piv, *some, *vij;
alpar@9:       int i, j, len, min_len, ncand, piv_lim = 5;
alpar@9:       double best, cost;
alpar@9:       /* nothing is chosen so far */
alpar@9:       piv = NULL, best = DBL_MAX, ncand = 0;
alpar@9:       /* if in the active submatrix there is a column that has the only
alpar@9:          non-zero (column singleton), choose it as a pivot */
alpar@9:       j = C_head[1];
alpar@9:       if (j != 0)
alpar@9:       {  xassert(C_len[j] == 1);
alpar@9:          piv = V_col[j];
alpar@9:          xassert(piv != NULL && piv->c_next == NULL);
alpar@9:          goto done;
alpar@9:       }
alpar@9:       /* if in the active submatrix there is a row that has the only
alpar@9:          non-zero (row singleton), choose it as a pivot */
alpar@9:       i = R_head[1];
alpar@9:       if (i != 0)
alpar@9:       {  xassert(R_len[i] == 1);
alpar@9:          piv = V_row[i];
alpar@9:          xassert(piv != NULL && piv->r_next == NULL);
alpar@9:          goto done;
alpar@9:       }
alpar@9:       /* there are no singletons in the active submatrix; walk through
alpar@9:          other non-empty rows and columns */
alpar@9:       for (len = 2; len <= n; len++)
alpar@9:       {  /* consider active columns having len non-zeros */
alpar@9:          for (j = C_head[len]; j != 0; j = C_next[j])
alpar@9:          {  /* j-th column has len non-zeros */
alpar@9:             /* find an element in the row of minimal length */
alpar@9:             some = NULL, min_len = INT_MAX;
alpar@9:             for (vij = V_col[j]; vij != NULL; vij = vij->c_next)
alpar@9:             {  if (min_len > R_len[vij->i])
alpar@9:                   some = vij, min_len = R_len[vij->i];
alpar@9:                /* if Markowitz cost of this element is not greater than
alpar@9:                   (len-1)**2, it can be chosen right now; this heuristic
alpar@9:                   reduces the search and works well in many cases */
alpar@9:                if (min_len <= len)
alpar@9:                {  piv = some;
alpar@9:                   goto done;
alpar@9:                }
alpar@9:             }
alpar@9:             /* j-th column has been scanned */
alpar@9:             /* the minimal element found is a next pivot candidate */
alpar@9:             xassert(some != NULL);
alpar@9:             ncand++;
alpar@9:             /* compute its Markowitz cost */
alpar@9:             cost = (double)(min_len - 1) * (double)(len - 1);
alpar@9:             /* choose between the current candidate and this element */
alpar@9:             if (cost < best) piv = some, best = cost;
alpar@9:             /* if piv_lim candidates have been considered, there is a
alpar@9:                doubt that a much better candidate exists; therefore it
alpar@9:                is the time to terminate the search */
alpar@9:             if (ncand == piv_lim) goto done;
alpar@9:          }
alpar@9:          /* now consider active rows having len non-zeros */
alpar@9:          for (i = R_head[len]; i != 0; i = R_next[i])
alpar@9:          {  /* i-th row has len non-zeros */
alpar@9:             /* find an element in the column of minimal length */
alpar@9:             some = NULL, min_len = INT_MAX;
alpar@9:             for (vij = V_row[i]; vij != NULL; vij = vij->r_next)
alpar@9:             {  if (min_len > C_len[vij->j])
alpar@9:                   some = vij, min_len = C_len[vij->j];
alpar@9:                /* if Markowitz cost of this element is not greater than
alpar@9:                   (len-1)**2, it can be chosen right now; this heuristic
alpar@9:                   reduces the search and works well in many cases */
alpar@9:                if (min_len <= len)
alpar@9:                {  piv = some;
alpar@9:                   goto done;
alpar@9:                }
alpar@9:             }
alpar@9:             /* i-th row has been scanned */
alpar@9:             /* the minimal element found is a next pivot candidate */
alpar@9:             xassert(some != NULL);
alpar@9:             ncand++;
alpar@9:             /* compute its Markowitz cost */
alpar@9:             cost = (double)(len - 1) * (double)(min_len - 1);
alpar@9:             /* choose between the current candidate and this element */
alpar@9:             if (cost < best) piv = some, best = cost;
alpar@9:             /* if piv_lim candidates have been considered, there is a
alpar@9:                doubt that a much better candidate exists; therefore it
alpar@9:                is the time to terminate the search */
alpar@9:             if (ncand == piv_lim) goto done;
alpar@9:          }
alpar@9:       }
alpar@9: done: /* bring the pivot v[p,q] to the factorizing routine */
alpar@9:       return piv;
alpar@9: }
alpar@9: 
alpar@9: /*----------------------------------------------------------------------
alpar@9: // eliminate - perform gaussian elimination.
alpar@9: //
alpar@9: // This routine performs elementary gaussian transformations in order
alpar@9: // to eliminate subdiagonal elements in the k-th column of the matrix
alpar@9: // U = P*V*Q using the pivot element u[k,k], where k is the number of
alpar@9: // the current elimination step.
alpar@9: //
alpar@9: // The parameter piv specifies the pivot element v[p,q] = u[k,k].
alpar@9: //
alpar@9: // Each time when the routine applies the elementary transformation to
alpar@9: // a non-pivot row of the matrix V, it stores the corresponding element
alpar@9: // to the matrix F in order to keep the main equality A = F*V.
alpar@9: //
alpar@9: // The routine assumes that on entry the matrices L = P*F*inv(P) and
alpar@9: // U = P*V*Q are the following:
alpar@9: //
alpar@9: //       1       k                  1       k         n
alpar@9: //    1  1 . . . . . . . . .     1  x x x x x x x x x x
alpar@9: //       x 1 . . . . . . . .        . x x x x x x x x x
alpar@9: //       x x 1 . . . . . . .        . . x x x x x x x x
alpar@9: //       x x x 1 . . . . . .        . . . x x x x x x x
alpar@9: //    k  x x x x 1 . . . . .     k  . . . . * * * * * *
alpar@9: //       x x x x _ 1 . . . .        . . . . # * * * * *
alpar@9: //       x x x x _ . 1 . . .        . . . . # * * * * *
alpar@9: //       x x x x _ . . 1 . .        . . . . # * * * * *
alpar@9: //       x x x x _ . . . 1 .        . . . . # * * * * *
alpar@9: //    n  x x x x _ . . . . 1     n  . . . . # * * * * *
alpar@9: //
alpar@9: //            matrix L                   matrix U
alpar@9: //
alpar@9: // where rows and columns of the matrix U with numbers k, k+1, ..., n
alpar@9: // form the active submatrix (eliminated elements are marked by '#' and
alpar@9: // other elements of the active submatrix are marked by '*'). Note that
alpar@9: // each eliminated non-zero element u[i,k] of the matrix U gives the
alpar@9: // corresponding element l[i,k] of the matrix L (marked by '_').
alpar@9: //
alpar@9: // Actually all operations are performed on the matrix V. Should note
alpar@9: // that the row-wise representation corresponds to the matrix V, but the
alpar@9: // column-wise representation corresponds to the active submatrix of the
alpar@9: // matrix V, i.e. elements of the matrix V, which doesn't belong to the
alpar@9: // active submatrix, are missing from the column linked lists.
alpar@9: //
alpar@9: // Let u[k,k] = v[p,q] be the pivot. In order to eliminate subdiagonal
alpar@9: // elements u[i',k] = v[i,q], i' = k+1, k+2, ..., n, the routine applies
alpar@9: // the following elementary gaussian transformations:
alpar@9: //
alpar@9: //    (i-th row of V) := (i-th row of V) - f[i,p] * (p-th row of V),
alpar@9: //
alpar@9: // where f[i,p] = v[i,q] / v[p,q] is a gaussian multiplier.
alpar@9: //
alpar@9: // Additionally, in order to keep the main equality A = F*V, each time
alpar@9: // when the routine applies the transformation to i-th row of the matrix
alpar@9: // V, it also adds f[i,p] as a new element to the matrix F.
alpar@9: //
alpar@9: // IMPORTANT: On entry the working arrays flag and work should contain
alpar@9: // zeros. This status is provided by the routine on exit. */
alpar@9: 
alpar@9: static void eliminate(LUX *lux, LUXWKA *wka, LUXELM *piv, int flag[],
alpar@9:       mpq_t work[])
alpar@9: {     DMP *pool = lux->pool;
alpar@9:       LUXELM **F_row = lux->F_row;
alpar@9:       LUXELM **F_col = lux->F_col;
alpar@9:       mpq_t *V_piv = lux->V_piv;
alpar@9:       LUXELM **V_row = lux->V_row;
alpar@9:       LUXELM **V_col = lux->V_col;
alpar@9:       int *R_len = wka->R_len;
alpar@9:       int *R_head = wka->R_head;
alpar@9:       int *R_prev = wka->R_prev;
alpar@9:       int *R_next = wka->R_next;
alpar@9:       int *C_len = wka->C_len;
alpar@9:       int *C_head = wka->C_head;
alpar@9:       int *C_prev = wka->C_prev;
alpar@9:       int *C_next = wka->C_next;
alpar@9:       LUXELM *fip, *vij, *vpj, *viq, *next;
alpar@9:       mpq_t temp;
alpar@9:       int i, j, p, q;
alpar@9:       mpq_init(temp);
alpar@9:       /* determine row and column indices of the pivot v[p,q] */
alpar@9:       xassert(piv != NULL);
alpar@9:       p = piv->i, q = piv->j;
alpar@9:       /* remove p-th (pivot) row from the active set; it will never
alpar@9:          return there */
alpar@9:       if (R_prev[p] == 0)
alpar@9:          R_head[R_len[p]] = R_next[p];
alpar@9:       else
alpar@9:          R_next[R_prev[p]] = R_next[p];
alpar@9:       if (R_next[p] == 0)
alpar@9:          ;
alpar@9:       else
alpar@9:          R_prev[R_next[p]] = R_prev[p];
alpar@9:       /* remove q-th (pivot) column from the active set; it will never
alpar@9:          return there */
alpar@9:       if (C_prev[q] == 0)
alpar@9:          C_head[C_len[q]] = C_next[q];
alpar@9:       else
alpar@9:          C_next[C_prev[q]] = C_next[q];
alpar@9:       if (C_next[q] == 0)
alpar@9:          ;
alpar@9:       else
alpar@9:          C_prev[C_next[q]] = C_prev[q];
alpar@9:       /* store the pivot value in a separate array */
alpar@9:       mpq_set(V_piv[p], piv->val);
alpar@9:       /* remove the pivot from p-th row */
alpar@9:       if (piv->r_prev == NULL)
alpar@9:          V_row[p] = piv->r_next;
alpar@9:       else
alpar@9:          piv->r_prev->r_next = piv->r_next;
alpar@9:       if (piv->r_next == NULL)
alpar@9:          ;
alpar@9:       else
alpar@9:          piv->r_next->r_prev = piv->r_prev;
alpar@9:       R_len[p]--;
alpar@9:       /* remove the pivot from q-th column */
alpar@9:       if (piv->c_prev == NULL)
alpar@9:          V_col[q] = piv->c_next;
alpar@9:       else
alpar@9:          piv->c_prev->c_next = piv->c_next;
alpar@9:       if (piv->c_next == NULL)
alpar@9:          ;
alpar@9:       else
alpar@9:          piv->c_next->c_prev = piv->c_prev;
alpar@9:       C_len[q]--;
alpar@9:       /* free the space occupied by the pivot */
alpar@9:       mpq_clear(piv->val);
alpar@9:       dmp_free_atom(pool, piv, sizeof(LUXELM));
alpar@9:       /* walk through p-th (pivot) row, which already does not contain
alpar@9:          the pivot v[p,q], and do the following... */
alpar@9:       for (vpj = V_row[p]; vpj != NULL; vpj = vpj->r_next)
alpar@9:       {  /* get column index of v[p,j] */
alpar@9:          j = vpj->j;
alpar@9:          /* store v[p,j] in the working array */
alpar@9:          flag[j] = 1;
alpar@9:          mpq_set(work[j], vpj->val);
alpar@9:          /* remove j-th column from the active set; it will return there
alpar@9:             later with a new length */
alpar@9:          if (C_prev[j] == 0)
alpar@9:             C_head[C_len[j]] = C_next[j];
alpar@9:          else
alpar@9:             C_next[C_prev[j]] = C_next[j];
alpar@9:          if (C_next[j] == 0)
alpar@9:             ;
alpar@9:          else
alpar@9:             C_prev[C_next[j]] = C_prev[j];
alpar@9:          /* v[p,j] leaves the active submatrix, so remove it from j-th
alpar@9:             column; however, v[p,j] is kept in p-th row */
alpar@9:          if (vpj->c_prev == NULL)
alpar@9:             V_col[j] = vpj->c_next;
alpar@9:          else
alpar@9:             vpj->c_prev->c_next = vpj->c_next;
alpar@9:          if (vpj->c_next == NULL)
alpar@9:             ;
alpar@9:          else
alpar@9:             vpj->c_next->c_prev = vpj->c_prev;
alpar@9:          C_len[j]--;
alpar@9:       }
alpar@9:       /* now walk through q-th (pivot) column, which already does not
alpar@9:          contain the pivot v[p,q], and perform gaussian elimination */
alpar@9:       while (V_col[q] != NULL)
alpar@9:       {  /* element v[i,q] has to be eliminated */
alpar@9:          viq = V_col[q];
alpar@9:          /* get row index of v[i,q] */
alpar@9:          i = viq->i;
alpar@9:          /* remove i-th row from the active set; later it will return
alpar@9:             there with a new length */
alpar@9:          if (R_prev[i] == 0)
alpar@9:             R_head[R_len[i]] = R_next[i];
alpar@9:          else
alpar@9:             R_next[R_prev[i]] = R_next[i];
alpar@9:          if (R_next[i] == 0)
alpar@9:             ;
alpar@9:          else
alpar@9:             R_prev[R_next[i]] = R_prev[i];
alpar@9:          /* compute gaussian multiplier f[i,p] = v[i,q] / v[p,q] and
alpar@9:             store it in the matrix F */
alpar@9:          fip = dmp_get_atom(pool, sizeof(LUXELM));
alpar@9:          fip->i = i, fip->j = p;
alpar@9:          mpq_init(fip->val);
alpar@9:          mpq_div(fip->val, viq->val, V_piv[p]);
alpar@9:          fip->r_prev = NULL;
alpar@9:          fip->r_next = F_row[i];
alpar@9:          fip->c_prev = NULL;
alpar@9:          fip->c_next = F_col[p];
alpar@9:          if (fip->r_next != NULL) fip->r_next->r_prev = fip;
alpar@9:          if (fip->c_next != NULL) fip->c_next->c_prev = fip;
alpar@9:          F_row[i] = F_col[p] = fip;
alpar@9:          /* v[i,q] has to be eliminated, so remove it from i-th row */
alpar@9:          if (viq->r_prev == NULL)
alpar@9:             V_row[i] = viq->r_next;
alpar@9:          else
alpar@9:             viq->r_prev->r_next = viq->r_next;
alpar@9:          if (viq->r_next == NULL)
alpar@9:             ;
alpar@9:          else
alpar@9:             viq->r_next->r_prev = viq->r_prev;
alpar@9:          R_len[i]--;
alpar@9:          /* and also from q-th column */
alpar@9:          V_col[q] = viq->c_next;
alpar@9:          C_len[q]--;
alpar@9:          /* free the space occupied by v[i,q] */
alpar@9:          mpq_clear(viq->val);
alpar@9:          dmp_free_atom(pool, viq, sizeof(LUXELM));
alpar@9:          /* perform gaussian transformation:
alpar@9:             (i-th row) := (i-th row) - f[i,p] * (p-th row)
alpar@9:             note that now p-th row, which is in the working array,
alpar@9:             does not contain the pivot v[p,q], and i-th row does not
alpar@9:             contain the element v[i,q] to be eliminated */
alpar@9:          /* walk through i-th row and transform existing non-zero
alpar@9:             elements */
alpar@9:          for (vij = V_row[i]; vij != NULL; vij = next)
alpar@9:          {  next = vij->r_next;
alpar@9:             /* get column index of v[i,j] */
alpar@9:             j = vij->j;
alpar@9:             /* v[i,j] := v[i,j] - f[i,p] * v[p,j] */
alpar@9:             if (flag[j])
alpar@9:             {  /* v[p,j] != 0 */
alpar@9:                flag[j] = 0;
alpar@9:                mpq_mul(temp, fip->val, work[j]);
alpar@9:                mpq_sub(vij->val, vij->val, temp);
alpar@9:                if (mpq_sgn(vij->val) == 0)
alpar@9:                {  /* new v[i,j] is zero, so remove it from the active
alpar@9:                      submatrix */
alpar@9:                   /* remove v[i,j] from i-th row */
alpar@9:                   if (vij->r_prev == NULL)
alpar@9:                      V_row[i] = vij->r_next;
alpar@9:                   else
alpar@9:                      vij->r_prev->r_next = vij->r_next;
alpar@9:                   if (vij->r_next == NULL)
alpar@9:                      ;
alpar@9:                   else
alpar@9:                      vij->r_next->r_prev = vij->r_prev;
alpar@9:                   R_len[i]--;
alpar@9:                   /* remove v[i,j] from j-th column */
alpar@9:                   if (vij->c_prev == NULL)
alpar@9:                      V_col[j] = vij->c_next;
alpar@9:                   else
alpar@9:                      vij->c_prev->c_next = vij->c_next;
alpar@9:                   if (vij->c_next == NULL)
alpar@9:                      ;
alpar@9:                   else
alpar@9:                      vij->c_next->c_prev = vij->c_prev;
alpar@9:                   C_len[j]--;
alpar@9:                   /* free the space occupied by v[i,j] */
alpar@9:                   mpq_clear(vij->val);
alpar@9:                   dmp_free_atom(pool, vij, sizeof(LUXELM));
alpar@9:                }
alpar@9:             }
alpar@9:          }
alpar@9:          /* now flag is the pattern of the set v[p,*] \ v[i,*] */
alpar@9:          /* walk through p-th (pivot) row and create new elements in
alpar@9:             i-th row, which appear due to fill-in */
alpar@9:          for (vpj = V_row[p]; vpj != NULL; vpj = vpj->r_next)
alpar@9:          {  j = vpj->j;
alpar@9:             if (flag[j])
alpar@9:             {  /* create new non-zero v[i,j] = 0 - f[i,p] * v[p,j] and
alpar@9:                   add it to i-th row and j-th column */
alpar@9:                vij = dmp_get_atom(pool, sizeof(LUXELM));
alpar@9:                vij->i = i, vij->j = j;
alpar@9:                mpq_init(vij->val);
alpar@9:                mpq_mul(vij->val, fip->val, work[j]);
alpar@9:                mpq_neg(vij->val, vij->val);
alpar@9:                vij->r_prev = NULL;
alpar@9:                vij->r_next = V_row[i];
alpar@9:                vij->c_prev = NULL;
alpar@9:                vij->c_next = V_col[j];
alpar@9:                if (vij->r_next != NULL) vij->r_next->r_prev = vij;
alpar@9:                if (vij->c_next != NULL) vij->c_next->c_prev = vij;
alpar@9:                V_row[i] = V_col[j] = vij;
alpar@9:                R_len[i]++, C_len[j]++;
alpar@9:             }
alpar@9:             else
alpar@9:             {  /* there is no fill-in, because v[i,j] already exists in
alpar@9:                   i-th row; restore the flag, which was reset before */
alpar@9:                flag[j] = 1;
alpar@9:             }
alpar@9:          }
alpar@9:          /* now i-th row has been completely transformed and can return
alpar@9:             to the active set with a new length */
alpar@9:          R_prev[i] = 0;
alpar@9:          R_next[i] = R_head[R_len[i]];
alpar@9:          if (R_next[i] != 0) R_prev[R_next[i]] = i;
alpar@9:          R_head[R_len[i]] = i;
alpar@9:       }
alpar@9:       /* at this point q-th (pivot) column must be empty */
alpar@9:       xassert(C_len[q] == 0);
alpar@9:       /* walk through p-th (pivot) row again and do the following... */
alpar@9:       for (vpj = V_row[p]; vpj != NULL; vpj = vpj->r_next)
alpar@9:       {  /* get column index of v[p,j] */
alpar@9:          j = vpj->j;
alpar@9:          /* erase v[p,j] from the working array */
alpar@9:          flag[j] = 0;
alpar@9:          mpq_set_si(work[j], 0, 1);
alpar@9:          /* now j-th column has been completely transformed, so it can
alpar@9:             return to the active list with a new length */
alpar@9:          C_prev[j] = 0;
alpar@9:          C_next[j] = C_head[C_len[j]];
alpar@9:          if (C_next[j] != 0) C_prev[C_next[j]] = j;
alpar@9:          C_head[C_len[j]] = j;
alpar@9:       }
alpar@9:       mpq_clear(temp);
alpar@9:       /* return to the factorizing routine */
alpar@9:       return;
alpar@9: }
alpar@9: 
alpar@9: /*----------------------------------------------------------------------
alpar@9: // lux_decomp - compute LU-factorization.
alpar@9: //
alpar@9: // SYNOPSIS
alpar@9: //
alpar@9: // #include "glplux.h"
alpar@9: // int lux_decomp(LUX *lux, int (*col)(void *info, int j, int ind[],
alpar@9: //    mpq_t val[]), void *info);
alpar@9: //
alpar@9: // DESCRIPTION
alpar@9: //
alpar@9: // The routine lux_decomp computes LU-factorization of a given square
alpar@9: // matrix A.
alpar@9: //
alpar@9: // The parameter lux specifies LU-factorization data structure built by
alpar@9: // means of the routine lux_create.
alpar@9: //
alpar@9: // The formal routine col specifies the original matrix A. In order to
alpar@9: // obtain j-th column of the matrix A the routine lux_decomp calls the
alpar@9: // routine col with the parameter j (1 <= j <= n, where n is the order
alpar@9: // of A). In response the routine col should store row indices and
alpar@9: // numerical values of non-zero elements of j-th column of A to the
alpar@9: // locations ind[1], ..., ind[len] and val[1], ..., val[len], resp.,
alpar@9: // where len is the number of non-zeros in j-th column, which should be
alpar@9: // returned on exit. Neiter zero nor duplicate elements are allowed.
alpar@9: //
alpar@9: // The parameter info is a transit pointer passed to the formal routine
alpar@9: // col; it can be used for various purposes.
alpar@9: //
alpar@9: // RETURNS
alpar@9: //
alpar@9: // The routine lux_decomp returns the singularity flag. Zero flag means
alpar@9: // that the original matrix A is non-singular while non-zero flag means
alpar@9: // that A is (exactly!) singular.
alpar@9: //
alpar@9: // Note that LU-factorization is valid in both cases, however, in case
alpar@9: // of singularity some rows of the matrix V (including pivot elements)
alpar@9: // will be empty.
alpar@9: //
alpar@9: // REPAIRING SINGULAR MATRIX
alpar@9: //
alpar@9: // If the routine lux_decomp returns non-zero flag, it provides all
alpar@9: // necessary information that can be used for "repairing" the matrix A,
alpar@9: // where "repairing" means replacing linearly dependent columns of the
alpar@9: // matrix A by appropriate columns of the unity matrix. This feature is
alpar@9: // needed when the routine lux_decomp is used for reinverting the basis
alpar@9: // matrix within the simplex method procedure.
alpar@9: //
alpar@9: // On exit linearly dependent columns of the matrix U have the numbers
alpar@9: // rank+1, rank+2, ..., n, where rank is the exact rank of the matrix A
alpar@9: // stored by the routine to the member lux->rank. The correspondence
alpar@9: // between columns of A and U is the same as between columns of V and U.
alpar@9: // Thus, linearly dependent columns of the matrix A have the numbers
alpar@9: // Q_col[rank+1], Q_col[rank+2], ..., Q_col[n], where Q_col is an array
alpar@9: // representing the permutation matrix Q in column-like format. It is
alpar@9: // understood that each j-th linearly dependent column of the matrix U
alpar@9: // should be replaced by the unity vector, where all elements are zero
alpar@9: // except the unity diagonal element u[j,j]. On the other hand j-th row
alpar@9: // of the matrix U corresponds to the row of the matrix V (and therefore
alpar@9: // of the matrix A) with the number P_row[j], where P_row is an array
alpar@9: // representing the permutation matrix P in row-like format. Thus, each
alpar@9: // j-th linearly dependent column of the matrix U should be replaced by
alpar@9: // a column of the unity matrix with the number P_row[j].
alpar@9: //
alpar@9: // The code that repairs the matrix A may look like follows:
alpar@9: //
alpar@9: //    for (j = rank+1; j <= n; j++)
alpar@9: //    {  replace column Q_col[j] of the matrix A by column P_row[j] of
alpar@9: //       the unity matrix;
alpar@9: //    }
alpar@9: //
alpar@9: // where rank, P_row, and Q_col are members of the structure LUX. */
alpar@9: 
alpar@9: int lux_decomp(LUX *lux, int (*col)(void *info, int j, int ind[],
alpar@9:       mpq_t val[]), void *info)
alpar@9: {     int n = lux->n;
alpar@9:       LUXELM **V_row = lux->V_row;
alpar@9:       LUXELM **V_col = lux->V_col;
alpar@9:       int *P_row = lux->P_row;
alpar@9:       int *P_col = lux->P_col;
alpar@9:       int *Q_row = lux->Q_row;
alpar@9:       int *Q_col = lux->Q_col;
alpar@9:       LUXELM *piv, *vij;
alpar@9:       LUXWKA *wka;
alpar@9:       int i, j, k, p, q, t, *flag;
alpar@9:       mpq_t *work;
alpar@9:       /* allocate working area */
alpar@9:       wka = xmalloc(sizeof(LUXWKA));
alpar@9:       wka->R_len = xcalloc(1+n, sizeof(int));
alpar@9:       wka->R_head = xcalloc(1+n, sizeof(int));
alpar@9:       wka->R_prev = xcalloc(1+n, sizeof(int));
alpar@9:       wka->R_next = xcalloc(1+n, sizeof(int));
alpar@9:       wka->C_len = xcalloc(1+n, sizeof(int));
alpar@9:       wka->C_head = xcalloc(1+n, sizeof(int));
alpar@9:       wka->C_prev = xcalloc(1+n, sizeof(int));
alpar@9:       wka->C_next = xcalloc(1+n, sizeof(int));
alpar@9:       /* initialize LU-factorization data structures */
alpar@9:       initialize(lux, col, info, wka);
alpar@9:       /* allocate working arrays */
alpar@9:       flag = xcalloc(1+n, sizeof(int));
alpar@9:       work = xcalloc(1+n, sizeof(mpq_t));
alpar@9:       for (k = 1; k <= n; k++)
alpar@9:       {  flag[k] = 0;
alpar@9:          mpq_init(work[k]);
alpar@9:       }
alpar@9:       /* main elimination loop */
alpar@9:       for (k = 1; k <= n; k++)
alpar@9:       {  /* choose a pivot element v[p,q] */
alpar@9:          piv = find_pivot(lux, wka);
alpar@9:          if (piv == NULL)
alpar@9:          {  /* no pivot can be chosen, because the active submatrix is
alpar@9:                empty */
alpar@9:             break;
alpar@9:          }
alpar@9:          /* determine row and column indices of the pivot element */
alpar@9:          p = piv->i, q = piv->j;
alpar@9:          /* let v[p,q] correspond to u[i',j']; permute k-th and i'-th
alpar@9:             rows and k-th and j'-th columns of the matrix U = P*V*Q to
alpar@9:             move the element u[i',j'] to the position u[k,k] */
alpar@9:          i = P_col[p], j = Q_row[q];
alpar@9:          xassert(k <= i && i <= n && k <= j && j <= n);
alpar@9:          /* permute k-th and i-th rows of the matrix U */
alpar@9:          t = P_row[k];
alpar@9:          P_row[i] = t, P_col[t] = i;
alpar@9:          P_row[k] = p, P_col[p] = k;
alpar@9:          /* permute k-th and j-th columns of the matrix U */
alpar@9:          t = Q_col[k];
alpar@9:          Q_col[j] = t, Q_row[t] = j;
alpar@9:          Q_col[k] = q, Q_row[q] = k;
alpar@9:          /* eliminate subdiagonal elements of k-th column of the matrix
alpar@9:             U = P*V*Q using the pivot element u[k,k] = v[p,q] */
alpar@9:          eliminate(lux, wka, piv, flag, work);
alpar@9:       }
alpar@9:       /* determine the rank of A (and V) */
alpar@9:       lux->rank = k - 1;
alpar@9:       /* free working arrays */
alpar@9:       xfree(flag);
alpar@9:       for (k = 1; k <= n; k++) mpq_clear(work[k]);
alpar@9:       xfree(work);
alpar@9:       /* build column lists of the matrix V using its row lists */
alpar@9:       for (j = 1; j <= n; j++)
alpar@9:          xassert(V_col[j] == NULL);
alpar@9:       for (i = 1; i <= n; i++)
alpar@9:       {  for (vij = V_row[i]; vij != NULL; vij = vij->r_next)
alpar@9:          {  j = vij->j;
alpar@9:             vij->c_prev = NULL;
alpar@9:             vij->c_next = V_col[j];
alpar@9:             if (vij->c_next != NULL) vij->c_next->c_prev = vij;
alpar@9:             V_col[j] = vij;
alpar@9:          }
alpar@9:       }
alpar@9:       /* free working area */
alpar@9:       xfree(wka->R_len);
alpar@9:       xfree(wka->R_head);
alpar@9:       xfree(wka->R_prev);
alpar@9:       xfree(wka->R_next);
alpar@9:       xfree(wka->C_len);
alpar@9:       xfree(wka->C_head);
alpar@9:       xfree(wka->C_prev);
alpar@9:       xfree(wka->C_next);
alpar@9:       xfree(wka);
alpar@9:       /* return to the calling program */
alpar@9:       return (lux->rank < n);
alpar@9: }
alpar@9: 
alpar@9: /*----------------------------------------------------------------------
alpar@9: // lux_f_solve - solve system F*x = b or F'*x = b.
alpar@9: //
alpar@9: // SYNOPSIS
alpar@9: //
alpar@9: // #include "glplux.h"
alpar@9: // void lux_f_solve(LUX *lux, int tr, mpq_t x[]);
alpar@9: //
alpar@9: // DESCRIPTION
alpar@9: //
alpar@9: // The routine lux_f_solve solves either the system F*x = b (if the
alpar@9: // flag tr is zero) or the system F'*x = b (if the flag tr is non-zero),
alpar@9: // where the matrix F is a component of LU-factorization specified by
alpar@9: // the parameter lux, F' is a matrix transposed to F.
alpar@9: //
alpar@9: // On entry the array x should contain elements of the right-hand side
alpar@9: // vector b in locations x[1], ..., x[n], where n is the order of the
alpar@9: // matrix F. On exit this array will contain elements of the solution
alpar@9: // vector x in the same locations. */
alpar@9: 
alpar@9: void lux_f_solve(LUX *lux, int tr, mpq_t x[])
alpar@9: {     int n = lux->n;
alpar@9:       LUXELM **F_row = lux->F_row;
alpar@9:       LUXELM **F_col = lux->F_col;
alpar@9:       int *P_row = lux->P_row;
alpar@9:       LUXELM *fik, *fkj;
alpar@9:       int i, j, k;
alpar@9:       mpq_t temp;
alpar@9:       mpq_init(temp);
alpar@9:       if (!tr)
alpar@9:       {  /* solve the system F*x = b */
alpar@9:          for (j = 1; j <= n; j++)
alpar@9:          {  k = P_row[j];
alpar@9:             if (mpq_sgn(x[k]) != 0)
alpar@9:             {  for (fik = F_col[k]; fik != NULL; fik = fik->c_next)
alpar@9:                {  mpq_mul(temp, fik->val, x[k]);
alpar@9:                   mpq_sub(x[fik->i], x[fik->i], temp);
alpar@9:                }
alpar@9:             }
alpar@9:          }
alpar@9:       }
alpar@9:       else
alpar@9:       {  /* solve the system F'*x = b */
alpar@9:          for (i = n; i >= 1; i--)
alpar@9:          {  k = P_row[i];
alpar@9:             if (mpq_sgn(x[k]) != 0)
alpar@9:             {  for (fkj = F_row[k]; fkj != NULL; fkj = fkj->r_next)
alpar@9:                {  mpq_mul(temp, fkj->val, x[k]);
alpar@9:                   mpq_sub(x[fkj->j], x[fkj->j], temp);
alpar@9:                }
alpar@9:             }
alpar@9:          }
alpar@9:       }
alpar@9:       mpq_clear(temp);
alpar@9:       return;
alpar@9: }
alpar@9: 
alpar@9: /*----------------------------------------------------------------------
alpar@9: // lux_v_solve - solve system V*x = b or V'*x = b.
alpar@9: //
alpar@9: // SYNOPSIS
alpar@9: //
alpar@9: // #include "glplux.h"
alpar@9: // void lux_v_solve(LUX *lux, int tr, double x[]);
alpar@9: //
alpar@9: // DESCRIPTION
alpar@9: //
alpar@9: // The routine lux_v_solve solves either the system V*x = b (if the
alpar@9: // flag tr is zero) or the system V'*x = b (if the flag tr is non-zero),
alpar@9: // where the matrix V is a component of LU-factorization specified by
alpar@9: // the parameter lux, V' is a matrix transposed to V.
alpar@9: //
alpar@9: // On entry the array x should contain elements of the right-hand side
alpar@9: // vector b in locations x[1], ..., x[n], where n is the order of the
alpar@9: // matrix V. On exit this array will contain elements of the solution
alpar@9: // vector x in the same locations. */
alpar@9: 
alpar@9: void lux_v_solve(LUX *lux, int tr, mpq_t x[])
alpar@9: {     int n = lux->n;
alpar@9:       mpq_t *V_piv = lux->V_piv;
alpar@9:       LUXELM **V_row = lux->V_row;
alpar@9:       LUXELM **V_col = lux->V_col;
alpar@9:       int *P_row = lux->P_row;
alpar@9:       int *Q_col = lux->Q_col;
alpar@9:       LUXELM *vij;
alpar@9:       int i, j, k;
alpar@9:       mpq_t *b, temp;
alpar@9:       b = xcalloc(1+n, sizeof(mpq_t));
alpar@9:       for (k = 1; k <= n; k++)
alpar@9:          mpq_init(b[k]), mpq_set(b[k], x[k]), mpq_set_si(x[k], 0, 1);
alpar@9:       mpq_init(temp);
alpar@9:       if (!tr)
alpar@9:       {  /* solve the system V*x = b */
alpar@9:          for (k = n; k >= 1; k--)
alpar@9:          {  i = P_row[k], j = Q_col[k];
alpar@9:             if (mpq_sgn(b[i]) != 0)
alpar@9:             {  mpq_set(x[j], b[i]);
alpar@9:                mpq_div(x[j], x[j], V_piv[i]);
alpar@9:                for (vij = V_col[j]; vij != NULL; vij = vij->c_next)
alpar@9:                {  mpq_mul(temp, vij->val, x[j]);
alpar@9:                   mpq_sub(b[vij->i], b[vij->i], temp);
alpar@9:                }
alpar@9:             }
alpar@9:          }
alpar@9:       }
alpar@9:       else
alpar@9:       {  /* solve the system V'*x = b */
alpar@9:          for (k = 1; k <= n; k++)
alpar@9:          {  i = P_row[k], j = Q_col[k];
alpar@9:             if (mpq_sgn(b[j]) != 0)
alpar@9:             {  mpq_set(x[i], b[j]);
alpar@9:                mpq_div(x[i], x[i], V_piv[i]);
alpar@9:                for (vij = V_row[i]; vij != NULL; vij = vij->r_next)
alpar@9:                {  mpq_mul(temp, vij->val, x[i]);
alpar@9:                   mpq_sub(b[vij->j], b[vij->j], temp);
alpar@9:                }
alpar@9:             }
alpar@9:          }
alpar@9:       }
alpar@9:       for (k = 1; k <= n; k++) mpq_clear(b[k]);
alpar@9:       mpq_clear(temp);
alpar@9:       xfree(b);
alpar@9:       return;
alpar@9: }
alpar@9: 
alpar@9: /*----------------------------------------------------------------------
alpar@9: // lux_solve - solve system A*x = b or A'*x = b.
alpar@9: //
alpar@9: // SYNOPSIS
alpar@9: //
alpar@9: // #include "glplux.h"
alpar@9: // void lux_solve(LUX *lux, int tr, mpq_t x[]);
alpar@9: //
alpar@9: // DESCRIPTION
alpar@9: //
alpar@9: // The routine lux_solve solves either the system A*x = b (if the flag
alpar@9: // tr is zero) or the system A'*x = b (if the flag tr is non-zero),
alpar@9: // where the parameter lux specifies LU-factorization of the matrix A,
alpar@9: // A' is a matrix transposed to A.
alpar@9: //
alpar@9: // On entry the array x should contain elements of the right-hand side
alpar@9: // vector b in locations x[1], ..., x[n], where n is the order of the
alpar@9: // matrix A. On exit this array will contain elements of the solution
alpar@9: // vector x in the same locations. */
alpar@9: 
alpar@9: void lux_solve(LUX *lux, int tr, mpq_t x[])
alpar@9: {     if (lux->rank < lux->n)
alpar@9:          xfault("lux_solve: LU-factorization has incomplete rank\n");
alpar@9:       if (!tr)
alpar@9:       {  /* A = F*V, therefore inv(A) = inv(V)*inv(F) */
alpar@9:          lux_f_solve(lux, 0, x);
alpar@9:          lux_v_solve(lux, 0, x);
alpar@9:       }
alpar@9:       else
alpar@9:       {  /* A' = V'*F', therefore inv(A') = inv(F')*inv(V') */
alpar@9:          lux_v_solve(lux, 1, x);
alpar@9:          lux_f_solve(lux, 1, x);
alpar@9:       }
alpar@9:       return;
alpar@9: }
alpar@9: 
alpar@9: /*----------------------------------------------------------------------
alpar@9: // lux_delete - delete LU-factorization.
alpar@9: //
alpar@9: // SYNOPSIS
alpar@9: //
alpar@9: // #include "glplux.h"
alpar@9: // void lux_delete(LUX *lux);
alpar@9: //
alpar@9: // DESCRIPTION
alpar@9: //
alpar@9: // The routine lux_delete deletes LU-factorization data structure,
alpar@9: // which the parameter lux points to, freeing all the memory allocated
alpar@9: // to this object. */
alpar@9: 
alpar@9: void lux_delete(LUX *lux)
alpar@9: {     int n = lux->n;
alpar@9:       LUXELM *fij, *vij;
alpar@9:       int i;
alpar@9:       for (i = 1; i <= n; i++)
alpar@9:       {  for (fij = lux->F_row[i]; fij != NULL; fij = fij->r_next)
alpar@9:             mpq_clear(fij->val);
alpar@9:          mpq_clear(lux->V_piv[i]);
alpar@9:          for (vij = lux->V_row[i]; vij != NULL; vij = vij->r_next)
alpar@9:             mpq_clear(vij->val);
alpar@9:       }
alpar@9:       dmp_delete_pool(lux->pool);
alpar@9:       xfree(lux->F_row);
alpar@9:       xfree(lux->F_col);
alpar@9:       xfree(lux->V_piv);
alpar@9:       xfree(lux->V_row);
alpar@9:       xfree(lux->V_col);
alpar@9:       xfree(lux->P_row);
alpar@9:       xfree(lux->P_col);
alpar@9:       xfree(lux->Q_row);
alpar@9:       xfree(lux->Q_col);
alpar@9:       xfree(lux);
alpar@9:       return;
alpar@9: }
alpar@9: 
alpar@9: /* eof */