path: root/test/monniaux/glpk-4.65/src/draft/lux.c

/* lux.c (LU-factorization, rational arithmetic) */

/***********************************************************************
*  This code is part of GLPK (GNU Linear Programming Kit).
*
*  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
*  2009, 2010, 2011, 2013 Andrew Makhorin, Department for Applied
*  Informatics, Moscow Aviation Institute, Moscow, Russia. All rights
*  reserved. E-mail: <mao@gnu.org>.
*
*  GLPK is free software: you can redistribute it and/or modify it
*  under the terms of the GNU General Public License as published by
*  the Free Software Foundation, either version 3 of the License, or
*  (at your option) any later version.
*
*  GLPK is distributed in the hope that it will be useful, but WITHOUT
*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
*  License for more details.
*
*  You should have received a copy of the GNU General Public License
*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
***********************************************************************/

#include "env.h"
#include "lux.h"

#define xfault xerror
#define dmp_create_poolx(size) dmp_create_pool()

/***********************************************************************
*  lux_create - create LU-factorization
*
*  SYNOPSIS
*
*  #include "lux.h"
*  LUX *lux_create(int n);
*
*  DESCRIPTION
*
*  The routine lux_create creates LU-factorization data structure for
*  a matrix of the order n. Initially the factorization corresponds to
*  the unity matrix (F = V = P = Q = I, so A = I).
*
*  RETURNS
*
*  The routine returns a pointer to the created LU-factorization data
*  structure, which represents the unity matrix of the order n. */

LUX *lux_create(int n)
{     LUX *lux;
      int k;
      if (n < 1)
         xfault("lux_create: n = %d; invalid parameter\n", n);
      lux = xmalloc(sizeof(LUX));
      lux->n = n;
      lux->pool = dmp_create_poolx(sizeof(LUXELM));
      lux->F_row = xcalloc(1+n, sizeof(LUXELM *));
      lux->F_col = xcalloc(1+n, sizeof(LUXELM *));
      lux->V_piv = xcalloc(1+n, sizeof(mpq_t));
      lux->V_row = xcalloc(1+n, sizeof(LUXELM *));
      lux->V_col = xcalloc(1+n, sizeof(LUXELM *));
      lux->P_row = xcalloc(1+n, sizeof(int));
      lux->P_col = xcalloc(1+n, sizeof(int));
      lux->Q_row = xcalloc(1+n, sizeof(int));
      lux->Q_col = xcalloc(1+n, sizeof(int));
      for (k = 1; k <= n; k++)
      {  lux->F_row[k] = lux->F_col[k] = NULL;
         mpq_init(lux->V_piv[k]);
         mpq_set_si(lux->V_piv[k], 1, 1);
         lux->V_row[k] = lux->V_col[k] = NULL;
         lux->P_row[k] = lux->P_col[k] = k;
         lux->Q_row[k] = lux->Q_col[k] = k;
      }
      lux->rank = n;
      return lux;
}

/***********************************************************************
*  initialize - initialize LU-factorization data structures
*
*  This routine initializes data structures for subsequent computing
*  the LU-factorization of a given matrix A, which is specified by the
*  formal routine col. On exit V = A and F = P = Q = I, where I is the
*  unity matrix. */

static void initialize(LUX *lux, int (*col)(void *info, int j,
      int ind[], mpq_t val[]), void *info, LUXWKA *wka)
{     int n = lux->n;
      DMP *pool = lux->pool;
      LUXELM **F_row = lux->F_row;
      LUXELM **F_col = lux->F_col;
      mpq_t *V_piv = lux->V_piv;
      LUXELM **V_row = lux->V_row;
      LUXELM **V_col = lux->V_col;
      int *P_row = lux->P_row;
      int *P_col = lux->P_col;
      int *Q_row = lux->Q_row;
      int *Q_col = lux->Q_col;
      int *R_len = wka->R_len;
      int *R_head = wka->R_head;
      int *R_prev = wka->R_prev;
      int *R_next = wka->R_next;
      int *C_len = wka->C_len;
      int *C_head = wka->C_head;
      int *C_prev = wka->C_prev;
      int *C_next = wka->C_next;
      LUXELM *fij, *vij;
      int i, j, k, len, *ind;
      mpq_t *val;
      /* F := I */
      for (i = 1; i <= n; i++)
      {  while (F_row[i] != NULL)
         {  fij = F_row[i], F_row[i] = fij->r_next;
            mpq_clear(fij->val);
            dmp_free_atom(pool, fij, sizeof(LUXELM));
         }
      }
      for (j = 1; j <= n; j++) F_col[j] = NULL;
      /* V := 0 */
      for (k = 1; k <= n; k++) mpq_set_si(V_piv[k], 0, 1);
      for (i = 1; i <= n; i++)
      {  while (V_row[i] != NULL)
         {  vij = V_row[i], V_row[i] = vij->r_next;
            mpq_clear(vij->val);
            dmp_free_atom(pool, vij, sizeof(LUXELM));
         }
      }
      for (j = 1; j <= n; j++) V_col[j] = NULL;
      /* V := A */
      ind = xcalloc(1+n, sizeof(int));
      val = xcalloc(1+n, sizeof(mpq_t));
      for (k = 1; k <= n; k++) mpq_init(val[k]);
      for (j = 1; j <= n; j++)
      {  /* obtain j-th column of matrix A */
         len = col(info, j, ind, val);
         if (!(0 <= len && len <= n))
            xfault("lux_decomp: j = %d: len = %d; invalid column length"
               "\n", j, len);
         /* copy elements of j-th column to matrix V */
         for (k = 1; k <= len; k++)
         {  /* get row index of a[i,j] */
            i = ind[k];
            if (!(1 <= i && i <= n))
               xfault("lux_decomp: j = %d: i = %d; row index out of ran"
                  "ge\n", j, i);
            /* check for duplicate indices */
            if (V_row[i] != NULL && V_row[i]->j == j)
               xfault("lux_decomp: j = %d: i = %d; duplicate row indice"
                  "s not allowed\n", j, i);
            /* check for zero value */
            if (mpq_sgn(val[k]) == 0)
               xfault("lux_decomp: j = %d: i = %d; zero elements not al"
                  "lowed\n", j, i);
            /* add new element v[i,j] = a[i,j] to V */
            vij = dmp_get_atom(pool, sizeof(LUXELM));
            vij->i = i, vij->j = j;
            mpq_init(vij->val);
            mpq_set(vij->val, val[k]);
            vij->r_prev = NULL;
            vij->r_next = V_row[i];
            vij->c_prev = NULL;
            vij->c_next = V_col[j];
            if (vij->r_next != NULL) vij->r_next->r_prev = vij;
            if (vij->c_next != NULL) vij->c_next->c_prev = vij;
            V_row[i] = V_col[j] = vij;
         }
      }
      xfree(ind);
      for (k = 1; k <= n; k++) mpq_clear(val[k]);
      xfree(val);
      /* P := Q := I */
      for (k = 1; k <= n; k++)
         P_row[k] = P_col[k] = Q_row[k] = Q_col[k] = k;
      /* the rank of A and V is not determined yet */
      lux->rank = -1;
      /* initially the entire matrix V is active */
      /* determine its row lengths */
      for (i = 1; i <= n; i++)
      {  len = 0;
         for (vij = V_row[i]; vij != NULL; vij = vij->r_next) len++;
         R_len[i] = len;
      }
      /* build linked lists of active rows */
      for (len = 0; len <= n; len++) R_head[len] = 0;
      for (i = 1; i <= n; i++)
      {  len = R_len[i];
         R_prev[i] = 0;
         R_next[i] = R_head[len];
         if (R_next[i] != 0) R_prev[R_next[i]] = i;
         R_head[len] = i;
      }
      /* determine its column lengths */
      for (j = 1; j <= n; j++)
      {  len = 0;
         for (vij = V_col[j]; vij != NULL; vij = vij->c_next) len++;
         C_len[j] = len;
      }
      /* build linked lists of active columns */
      for (len = 0; len <= n; len++) C_head[len] = 0;
      for (j = 1; j <= n; j++)
      {  len = C_len[j];
         C_prev[j] = 0;
         C_next[j] = C_head[len];
         if (C_next[j] != 0) C_prev[C_next[j]] = j;
         C_head[len] = j;
      }
      return;
}

/***********************************************************************
*  find_pivot - choose a pivot element
*
*  This routine chooses a pivot element v[p,q] in the active submatrix
*  of matrix U = P*V*Q.
*
*  It is assumed that on entry the matrix U has the following partially
*  triangularized form:
*
*        1       k         n
*     1  x x x x x x x x x x
*        . x x x x x x x x x
*        . . x x x x x x x x
*        . . . x x x x x x x
*     k  . . . . * * * * * *
*        . . . . * * * * * *
*        . . . . * * * * * *
*        . . . . * * * * * *
*        . . . . * * * * * *
*     n  . . . . * * * * * *
*
*  where rows and columns k, k+1, ..., n belong to the active submatrix
*  (elements of the active submatrix are marked by '*').
*
*  Since the matrix U = P*V*Q is not stored, the routine works with the
*  matrix V. It is assumed that the row-wise representation corresponds
*  to the matrix V, but the column-wise representation corresponds to
*  the active submatrix of the matrix V, i.e. elements of the matrix V,
*  which does not belong to the active submatrix, are missing from the
*  column linked lists. It is also assumed that each active row of the
*  matrix V is in the set R[len], where len is number of non-zeros in
*  the row, and each active column of the matrix V is in the set C[len],
*  where len is number of non-zeros in the column (in the latter case
*  only elements of the active submatrix are counted; such elements are
*  marked by '*' on the figure above).
*
*  Due to exact arithmetic any non-zero element of the active submatrix
*  can be chosen as a pivot. However, to keep sparsity of the matrix V
*  the routine uses Markowitz strategy, trying to choose such element
*  v[p,q], which has smallest Markowitz cost (nr[p]-1) * (nc[q]-1),
*  where nr[p] and nc[q] are the number of non-zero elements, resp., in
*  p-th row and in q-th column of the active submatrix.
*
*  In order to reduce the search, i.e. not to walk through all elements
*  of the active submatrix, the routine exploits a technique proposed by
*  I.Duff. This technique is based on using the sets R[len] and C[len]
*  of active rows and columns.
*
*  On exit the routine returns a pointer to a pivot v[p,q] chosen, or
*  NULL, if the active submatrix is empty. */

static LUXELM *find_pivot(LUX *lux, LUXWKA *wka)
{     int n = lux->n;
      LUXELM **V_row = lux->V_row;
      LUXELM **V_col = lux->V_col;
      int *R_len = wka->R_len;
      int *R_head = wka->R_head;
      int *R_next = wka->R_next;
      int *C_len = wka->C_len;
      int *C_head = wka->C_head;
      int *C_next = wka->C_next;
      LUXELM *piv, *some, *vij;
      int i, j, len, min_len, ncand, piv_lim = 5;
      double best, cost;
      /* nothing is chosen so far */
      piv = NULL, best = DBL_MAX, ncand = 0;
      /* if in the active submatrix there is a column that has the only
         non-zero (column singleton), choose it as a pivot */
      j = C_head[1];
      if (j != 0)
      {  xassert(C_len[j] == 1);
         piv = V_col[j];
         xassert(piv != NULL && piv->c_next == NULL);
         goto done;
      }
      /* if in the active submatrix there is a row that has the only
         non-zero (row singleton), choose it as a pivot */
      i = R_head[1];
      if (i != 0)
      {  xassert(R_len[i] == 1);
         piv = V_row[i];
         xassert(piv != NULL && piv->r_next == NULL);
         goto done;
      }
      /* there are no singletons in the active submatrix; walk through
         other non-empty rows and columns */
      for (len = 2; len <= n; len++)
      {  /* consider active columns having len non-zeros */
         for (j = C_head[len]; j != 0; j = C_next[j])
         {  /* j-th column has len non-zeros */
            /* find an element in the row of minimal length */
            some = NULL, min_len = INT_MAX;
            for (vij = V_col[j]; vij != NULL; vij = vij->c_next)
            {  if (min_len > R_len[vij->i])
                  some = vij, min_len = R_len[vij->i];
               /* if Markowitz cost of this element is not greater than
                  (len-1)**2, it can be chosen right now; this heuristic
                  reduces the search and works well in many cases */
               if (min_len <= len)
               {  piv = some;
                  goto done;
               }
            }
            /* j-th column has been scanned */
            /* the minimal element found is a next pivot candidate */
            xassert(some != NULL);
            ncand++;
            /* compute its Markowitz cost */
            cost = (double)(min_len - 1) * (double)(len - 1);
            /* choose between the current candidate and this element */
            if (cost < best) piv = some, best = cost;
            /* if piv_lim candidates have been considered, there is a
               doubt that a much better candidate exists; therefore it
               is the time to terminate the search */
            if (ncand == piv_lim) goto done;
         }
         /* now consider active rows having len non-zeros */
         for (i = R_head[len]; i != 0; i = R_next[i])
         {  /* i-th row has len non-zeros */
            /* find an element in the column of minimal length */
            some = NULL, min_len = INT_MAX;
            for (vij = V_row[i]; vij != NULL; vij = vij->r_next)
            {  if (min_len > C_len[vij->j])
                  some = vij, min_len = C_len[vij->j];
               /* if Markowitz cost of this element is not greater than
                  (len-1)**2, it can be chosen right now; this heuristic
                  reduces the search and works well in many cases */
               if (min_len <= len)
               {  piv = some;
                  goto done;
               }
            }
            /* i-th row has been scanned */
            /* the minimal element found is a next pivot candidate */
            xassert(some != NULL);
            ncand++;
            /* compute its Markowitz cost */
            cost = (double)(len - 1) * (double)(min_len - 1);
            /* choose between the current candidate and this element */
            if (cost < best) piv = some, best = cost;
            /* if piv_lim candidates have been considered, there is a
               doubt that a much better candidate exists; therefore it
               is the time to terminate the search */
            if (ncand == piv_lim) goto done;
         }
      }
done: /* bring the pivot v[p,q] to the factorizing routine */
      return piv;
}

/***********************************************************************
*  eliminate - perform gaussian elimination
*
*  This routine performs elementary gaussian transformations in order
*  to eliminate subdiagonal elements in the k-th column of the matrix
*  U = P*V*Q using the pivot element u[k,k], where k is the number of
*  the current elimination step.
*
*  The parameter piv specifies the pivot element v[p,q] = u[k,k].
*
*  Each time when the routine applies the elementary transformation to
*  a non-pivot row of the matrix V, it stores the corresponding element
*  to the matrix F in order to keep the main equality A = F*V.
*
*  The routine assumes that on entry the matrices L = P*F*inv(P) and
*  U = P*V*Q are the following:
*
*        1       k                  1       k         n
*     1  1 . . . . . . . . .     1  x x x x x x x x x x
*        x 1 . . . . . . . .        . x x x x x x x x x
*        x x 1 . . . . . . .        . . x x x x x x x x
*        x x x 1 . . . . . .        . . . x x x x x x x
*     k  x x x x 1 . . . . .     k  . . . . * * * * * *
*        x x x x _ 1 . . . .        . . . . # * * * * *
*        x x x x _ . 1 . . .        . . . . # * * * * *
*        x x x x _ . . 1 . .        . . . . # * * * * *
*        x x x x _ . . . 1 .        . . . . # * * * * *
*     n  x x x x _ . . . . 1     n  . . . . # * * * * *
*
*             matrix L                   matrix U
*
*  where rows and columns of the matrix U with numbers k, k+1, ..., n
*  form the active submatrix (eliminated elements are marked by '#' and
*  other elements of the active submatrix are marked by '*'). Note that
*  each eliminated non-zero element u[i,k] of the matrix U gives the
*  corresponding element l[i,k] of the matrix L (marked by '_').
*
*  Actually all operations are performed on the matrix V. Should note
*  that the row-wise representation corresponds to the matrix V, but the
*  column-wise representation corresponds to the active submatrix of the
*  matrix V, i.e. elements of the matrix V, which doesn't belong to the
*  active submatrix, are missing from the column linked lists.
*
*  Let u[k,k] = v[p,q] be the pivot. In order to eliminate subdiagonal
*  elements u[i',k] = v[i,q], i' = k+1, k+2, ..., n, the routine applies
*  the following elementary gaussian transformations:
*
*     (i-th row of V) := (i-th row of V) - f[i,p] * (p-th row of V),
*
*  where f[i,p] = v[i,q] / v[p,q] is a gaussian multiplier.
*
*  Additionally, in order to keep the main equality A = F*V, each time
*  when the routine applies the transformation to i-th row of the matrix
*  V, it also adds f[i,p] as a new element to the matrix F.
*
*  IMPORTANT: On entry the working arrays flag and work should contain
*  zeros. This status is provided by the routine on exit. */

static void eliminate(LUX *lux, LUXWKA *wka, LUXELM *piv, int flag[],
      mpq_t work[])
{     DMP *pool = lux->pool;
      LUXELM **F_row = lux->F_row;
      LUXELM **F_col = lux->F_col;
      mpq_t *V_piv = lux->V_piv;
      LUXELM **V_row = lux->V_row;
      LUXELM **V_col = lux->V_col;
      int *R_len = wka->R_len;
      int *R_head = wka->R_head;
      int *R_prev = wka->R_prev;
      int *R_next = wka->R_next;
      int *C_len = wka->C_len;
      int *C_head = wka->C_head;
      int *C_prev = wka->C_prev;
      int *C_next = wka->C_next;
      LUXELM *fip, *vij, *vpj, *viq, *next;
      mpq_t temp;
      int i, j, p, q;
      mpq_init(temp);
      /* determine row and column indices of the pivot v[p,q] */
      xassert(piv != NULL);
      p = piv->i, q = piv->j;
      /* remove p-th (pivot) row from the active set; it will never
         return there */
      if (R_prev[p] == 0)
         R_head[R_len[p]] = R_next[p];
      else
         R_next[R_prev[p]] = R_next[p];
      if (R_next[p] == 0)
         ;
      else
         R_prev[R_next[p]] = R_prev[p];
      /* remove q-th (pivot) column from the active set; it will never
         return there */
      if (C_prev[q] == 0)
         C_head[C_len[q]] = C_next[q];
      else
         C_next[C_prev[q]] = C_next[q];
      if (C_next[q] == 0)
         ;
      else
         C_prev[C_next[q]] = C_prev[q];
      /* store the pivot value in a separate array */
      mpq_set(V_piv[p], piv->val);
      /* remove the pivot from p-th row */
      if (piv->r_prev == NULL)
         V_row[p] = piv->r_next;
      else
         piv->r_prev->r_next = piv->r_next;
      if (piv->r_next == NULL)
         ;
      else
         piv->r_next->r_prev = piv->r_prev;
      R_len[p]--;
      /* remove the pivot from q-th column */
      if (piv->c_prev == NULL)
         V_col[q] = piv->c_next;
      else
         piv->c_prev->c_next = piv->c_next;
      if (piv->c_next == NULL)
         ;
      else
         piv->c_next->c_prev = piv->c_prev;
      C_len[q]--;
      /* free the space occupied by the pivot */
      mpq_clear(piv->val);
      dmp_free_atom(pool, piv, sizeof(LUXELM));
      /* walk through p-th (pivot) row, which already does not contain
         the pivot v[p,q], and do the following... */
      for (vpj = V_row[p]; vpj != NULL; vpj = vpj->r_next)
      {  /* get column index of v[p,j] */
         j = vpj->j;
         /* store v[p,j] in the working array */
         flag[j] = 1;
         mpq_set(work[j], vpj->val);
         /* remove j-th column from the active set; it will return there
            later with a new length */
         if (C_prev[j] == 0)
            C_head[C_len[j]] = C_next[j];
         else
            C_next[C_prev[j]] = C_next[j];
         if (C_next[j] == 0)
            ;
         else
            C_prev[C_next[j]] = C_prev[j];
         /* v[p,j] leaves the active submatrix, so remove it from j-th
            column; however, v[p,j] is kept in p-th row */
         if (vpj->c_prev == NULL)
            V_col[j] = vpj->c_next;
         else
            vpj->c_prev->c_next = vpj->c_next;
         if (vpj->c_next == NULL)
            ;
         else
            vpj->c_next->c_prev = vpj->c_prev;
         C_len[j]--;
      }
      /* now walk through q-th (pivot) column, which already does not
         contain the pivot v[p,q], and perform gaussian elimination */
      while (V_col[q] != NULL)
      {  /* element v[i,q] has to be eliminated */
         viq = V_col[q];
         /* get row index of v[i,q] */
         i = viq->i;
         /* remove i-th row from the active set; later it will return
            there with a new length */
         if (R_prev[i] == 0)
            R_head[R_len[i]] = R_next[i];
         else
            R_next[R_prev[i]] = R_next[i];
         if (R_next[i] == 0)
            ;
         else
            R_prev[R_next[i]] = R_prev[i];
         /* compute gaussian multiplier f[i,p] = v[i,q] / v[p,q] and
            store it in the matrix F */
         fip = dmp_get_atom(pool, sizeof(LUXELM));
         fip->i = i, fip->j = p;
         mpq_init(fip->val);
         mpq_div(fip->val, viq->val, V_piv[p]);
         fip->r_prev = NULL;
         fip->r_next = F_row[i];
         fip->c_prev = NULL;
         fip->c_next = F_col[p];
         if (fip->r_next != NULL) fip->r_next->r_prev = fip;
         if (fip->c_next != NULL) fip->c_next->c_prev = fip;
         F_row[i] = F_col[p] = fip;
         /* v[i,q] has to be eliminated, so remove it from i-th row */
         if (viq->r_prev == NULL)
            V_row[i] = viq->r_next;
         else
            viq->r_prev->r_next = viq->r_next;
         if (viq->r_next == NULL)
            ;
         else
            viq->r_next->r_prev = viq->r_prev;
         R_len[i]--;
         /* and also from q-th column */
         V_col[q] = viq->c_next;
         C_len[q]--;
         /* free the space occupied by v[i,q] */
         mpq_clear(viq->val);
         dmp_free_atom(pool, viq, sizeof(LUXELM));
         /* perform gaussian transformation:
            (i-th row) := (i-th row) - f[i,p] * (p-th row)
            note that now p-th row, which is in the working array,
            does not contain the pivot v[p,q], and i-th row does not
            contain the element v[i,q] to be eliminated */
         /* walk through i-th row and transform existing non-zero
            elements */
         for (vij = V_row[i]; vij != NULL; vij = next)
         {  next = vij->r_next;
            /* get column index of v[i,j] */
            j = vij->j;
            /* v[i,j] := v[i,j] - f[i,p] * v[p,j] */
            if (flag[j])
            {  /* v[p,j] != 0 */
               flag[j] = 0;
               mpq_mul(temp, fip->val, work[j]);
               mpq_sub(vij->val, vij->val, temp);
               if (mpq_sgn(vij->val) == 0)
               {  /* new v[i,j] is zero, so remove it from the active
                     submatrix */
                  /* remove v[i,j] from i-th row */
                  if (vij->r_prev == NULL)
                     V_row[i] = vij->r_next;
                  else
                     vij->r_prev->r_next = vij->r_next;
                  if (vij->r_next == NULL)
                     ;
                  else
                     vij->r_next->r_prev = vij->r_prev;
                  R_len[i]--;
                  /* remove v[i,j] from j-th column */
                  if (vij->c_prev == NULL)
                     V_col[j] = vij->c_next;
                  else
                     vij->c_prev->c_next = vij->c_next;
                  if (vij->c_next == NULL)
                     ;
                  else
                     vij->c_next->c_prev = vij->c_prev;
                  C_len[j]--;
                  /* free the space occupied by v[i,j] */
                  mpq_clear(vij->val);
                  dmp_free_atom(pool, vij, sizeof(LUXELM));
               }
            }
         }
         /* now flag is the pattern of the set v[p,*] \ v[i,*] */
         /* walk through p-th (pivot) row and create new elements in
            i-th row, which appear due to fill-in */
         for (vpj = V_row[p]; vpj != NULL; vpj = vpj->r_next)
         {  j = vpj->j;
            if (flag[j])
            {  /* create new non-zero v[i,j] = 0 - f[i,p] * v[p,j] and
                  add it to i-th row and j-th column */
               vij = dmp_get_atom(pool, sizeof(LUXELM));
               vij->i = i, vij->j = j;
               mpq_init(vij->val);
               mpq_mul(vij->val, fip->val, work[j]);
               mpq_neg(vij->val, vij->val);
               vij->r_prev = NULL;
               vij->r_next = V_row[i];
               vij->c_prev = NULL;
               vij->c_next = V_col[j];
               if (vij->r_next != NULL) vij->r_next->r_prev = vij;
               if (vij->c_next != NULL) vij->c_next->c_prev = vij;
               V_row[i] = V_col[j] = vij;
               R_len[i]++, C_len[j]++;
            }
            else
            {  /* there is no fill-in, because v[i,j] already exists in
                  i-th row; restore the flag, which was reset before */
               flag[j] = 1;
            }
         }
         /* now i-th row has been completely transformed and can return
            to the active set with a new length */
         R_prev[i] = 0;
         R_next[i] = R_head[R_len[i]];
         if (R_next[i] != 0) R_prev[R_next[i]] = i;
         R_head[R_len[i]] = i;
      }
      /* at this point q-th (pivot) column must be empty */
      xassert(C_len[q] == 0);
      /* walk through p-th (pivot) row again and do the following... */
      for (vpj = V_row[p]; vpj != NULL; vpj = vpj->r_next)
      {  /* get column index of v[p,j] */
         j = vpj->j;
         /* erase v[p,j] from the working array */
         flag[j] = 0;
         mpq_set_si(work[j], 0, 1);
         /* now j-th column has been completely transformed, so it can
            return to the active list with a new length */
         C_prev[j] = 0;
         C_next[j] = C_head[C_len[j]];
         if (C_next[j] != 0) C_prev[C_next[j]] = j;
         C_head[C_len[j]] = j;
      }
      mpq_clear(temp);
      /* return to the factorizing routine */
      return;
}

/***********************************************************************
*  lux_decomp - compute LU-factorization
*
*  SYNOPSIS
*
*  #include "lux.h"
*  int lux_decomp(LUX *lux, int (*col)(void *info, int j, int ind[],
*     mpq_t val[]), void *info);
*
*  DESCRIPTION
*
*  The routine lux_decomp computes LU-factorization of a given square
*  matrix A.
*
*  The parameter lux specifies LU-factorization data structure built by
*  means of the routine lux_create.
*
*  The formal routine col specifies the original matrix A. In order to
*  obtain j-th column of the matrix A the routine lux_decomp calls the
*  routine col with the parameter j (1 <= j <= n, where n is the order
*  of A). In response the routine col should store row indices and
*  numerical values of non-zero elements of j-th column of A to the
*  locations ind[1], ..., ind[len] and val[1], ..., val[len], resp.,
*  where len is the number of non-zeros in j-th column, which should be
*  returned on exit. Neiter zero nor duplicate elements are allowed.
*
*  The parameter info is a transit pointer passed to the formal routine
*  col; it can be used for various purposes.
*
*  RETURNS
*
*  The routine lux_decomp returns the singularity flag. Zero flag means
*  that the original matrix A is non-singular while non-zero flag means
*  that A is (exactly!) singular.
*
*  Note that LU-factorization is valid in both cases, however, in case
*  of singularity some rows of the matrix V (including pivot elements)
*  will be empty.
*
*  REPAIRING SINGULAR MATRIX
*
*  If the routine lux_decomp returns non-zero flag, it provides all
*  necessary information that can be used for "repairing" the matrix A,
*  where "repairing" means replacing linearly dependent columns of the
*  matrix A by appropriate columns of the unity matrix. This feature is
*  needed when the routine lux_decomp is used for reinverting the basis
*  matrix within the simplex method procedure.
*
*  On exit linearly dependent columns of the matrix U have the numbers
*  rank+1, rank+2, ..., n, where rank is the exact rank of the matrix A
*  stored by the routine to the member lux->rank. The correspondence
*  between columns of A and U is the same as between columns of V and U.
*  Thus, linearly dependent columns of the matrix A have the numbers
*  Q_col[rank+1], Q_col[rank+2], ..., Q_col[n], where Q_col is an array
*  representing the permutation matrix Q in column-like format. It is
*  understood that each j-th linearly dependent column of the matrix U
*  should be replaced by the unity vector, where all elements are zero
*  except the unity diagonal element u[j,j]. On the other hand j-th row
*  of the matrix U corresponds to the row of the matrix V (and therefore
*  of the matrix A) with the number P_row[j], where P_row is an array
*  representing the permutation matrix P in row-like format. Thus, each
*  j-th linearly dependent column of the matrix U should be replaced by
*  a column of the unity matrix with the number P_row[j].
*
*  The code that repairs the matrix A may look like follows:
*
*     for (j = rank+1; j <= n; j++)
*     {  replace column Q_col[j] of the matrix A by column P_row[j] of
*        the unity matrix;
*     }
*
*  where rank, P_row, and Q_col are members of the structure LUX. */

int lux_decomp(LUX *lux, int (*col)(void *info, int j, int ind[],
      mpq_t val[]), void *info)
{     int n = lux->n;
      LUXELM **V_row = lux->V_row;
      LUXELM **V_col = lux->V_col;
      int *P_row = lux->P_row;
      int *P_col = lux->P_col;
      int *Q_row = lux->Q_row;
      int *Q_col = lux->Q_col;
      LUXELM *piv, *vij;
      LUXWKA *wka;
      int i, j, k, p, q, t, *flag;
      mpq_t *work;
      /* allocate working area */
      wka = xmalloc(sizeof(LUXWKA));
      wka->R_len = xcalloc(1+n, sizeof(int));
      wka->R_head = xcalloc(1+n, sizeof(int));
      wka->R_prev = xcalloc(1+n, sizeof(int));
      wka->R_next = xcalloc(1+n, sizeof(int));
      wka->C_len = xcalloc(1+n, sizeof(int));
      wka->C_head = xcalloc(1+n, sizeof(int));
      wka->C_prev = xcalloc(1+n, sizeof(int));
      wka->C_next = xcalloc(1+n, sizeof(int));
      /* initialize LU-factorization data structures */
      initialize(lux, col, info, wka);
      /* allocate working arrays */
      flag = xcalloc(1+n, sizeof(int));
      work = xcalloc(1+n, sizeof(mpq_t));
      for (k = 1; k <= n; k++)
      {  flag[k] = 0;
         mpq_init(work[k]);
      }
      /* main elimination loop */
      for (k = 1; k <= n; k++)
      {  /* choose a pivot element v[p,q] */
         piv = find_pivot(lux, wka);
         if (piv == NULL)
         {  /* no pivot can be chosen, because the active submatrix is
               empty */
            break;
         }
         /* determine row and column indices of the pivot element */
         p = piv->i, q = piv->j;
         /* let v[p,q] correspond to u[i',j']; permute k-th and i'-th
            rows and k-th and j'-th columns of the matrix U = P*V*Q to
            move the element u[i',j'] to the position u[k,k] */
         i = P_col[p], j = Q_row[q];
         xassert(k <= i && i <= n && k <= j && j <= n);
         /* permute k-th and i-th rows of the matrix U */
         t = P_row[k];
         P_row[i] = t, P_col[t] = i;
         P_row[k] = p, P_col[p] = k;
         /* permute k-th and j-th columns of the matrix U */
         t = Q_col[k];
         Q_col[j] = t, Q_row[t] = j;
         Q_col[k] = q, Q_row[q] = k;
         /* eliminate subdiagonal elements of k-th column of the matrix
            U = P*V*Q using the pivot element u[k,k] = v[p,q] */
         eliminate(lux, wka, piv, flag, work);
      }
      /* determine the rank of A (and V) */
      lux->rank = k - 1;
      /* free working arrays */
      xfree(flag);
      for (k = 1; k <= n; k++) mpq_clear(work[k]);
      xfree(work);
      /* build column lists of the matrix V using its row lists */
      for (j = 1; j <= n; j++)
         xassert(V_col[j] == NULL);
      for (i = 1; i <= n; i++)
      {  for (vij = V_row[i]; vij != NULL; vij = vij->r_next)
         {  j = vij->j;
            vij->c_prev = NULL;
            vij->c_next = V_col[j];
            if (vij->c_next != NULL) vij->c_next->c_prev = vij;
            V_col[j] = vij;
         }
      }
      /* free working area */
      xfree(wka->R_len);
      xfree(wka->R_head);
      xfree(wka->R_prev);
      xfree(wka->R_next);
      xfree(wka->C_len);
      xfree(wka->C_head);
      xfree(wka->C_prev);
      xfree(wka->C_next);
      xfree(wka);
      /* return to the calling program */
      return (lux->rank < n);
}

/***********************************************************************
*  lux_f_solve - solve system F*x = b or F'*x = b
*
*  SYNOPSIS
*
*  #include "lux.h"
*  void lux_f_solve(LUX *lux, int tr, mpq_t x[]);
*
*  DESCRIPTION
*
*  The routine lux_f_solve solves either the system F*x = b (if the
*  flag tr is zero) or the system F'*x = b (if the flag tr is non-zero),
*  where the matrix F is a component of LU-factorization specified by
*  the parameter lux, F' is a matrix transposed to F.
*
*  On entry the array x should contain elements of the right-hand side
*  vector b in locations x[1], ..., x[n], where n is the order of the
*  matrix F. On exit this array will contain elements of the solution
*  vector x in the same locations. */

void lux_f_solve(LUX *lux, int tr, mpq_t x[])
{     int n = lux->n;
      LUXELM **F_row = lux->F_row;
      LUXELM **F_col = lux->F_col;
      int *P_row = lux->P_row;
      LUXELM *fik, *fkj;
      int i, j, k;
      mpq_t temp;
      mpq_init(temp);
      if (!tr)
      {  /* solve the system F*x = b */
         for (j = 1; j <= n; j++)
         {  k = P_row[j];
            if (mpq_sgn(x[k]) != 0)
            {  for (fik = F_col[k]; fik != NULL; fik = fik->c_next)
               {  mpq_mul(temp, fik->val, x[k]);
                  mpq_sub(x[fik->i], x[fik->i], temp);
               }
            }
         }
      }
      else
      {  /* solve the system F'*x = b */
         for (i = n; i >= 1; i--)
         {  k = P_row[i];
            if (mpq_sgn(x[k]) != 0)
            {  for (fkj = F_row[k]; fkj != NULL; fkj = fkj->r_next)
               {  mpq_mul(temp, fkj->val, x[k]);
                  mpq_sub(x[fkj->j], x[fkj->j], temp);
               }
            }
         }
      }
      mpq_clear(temp);
      return;
}

/***********************************************************************
*  lux_v_solve - solve system V*x = b or V'*x = b
*
*  SYNOPSIS
*
*  #include "lux.h"
*  void lux_v_solve(LUX *lux, int tr, double x[]);
*
*  DESCRIPTION
*
*  The routine lux_v_solve solves either the system V*x = b (if the
*  flag tr is zero) or the system V'*x = b (if the flag tr is non-zero),
*  where the matrix V is a component of LU-factorization specified by
*  the parameter lux, V' is a matrix transposed to V.
*
*  On entry the array x should contain elements of the right-hand side
*  vector b in locations x[1], ..., x[n], where n is the order of the
*  matrix V. On exit this array will contain elements of the solution
*  vector x in the same locations. */

void lux_v_solve(LUX *lux, int tr, mpq_t x[])
{     int n = lux->n;
      mpq_t *V_piv = lux->V_piv;
      LUXELM **V_row = lux->V_row;
      LUXELM **V_col = lux->V_col;
      int *P_row = lux->P_row;
      int *Q_col = lux->Q_col;
      LUXELM *vij;
      int i, j, k;
      mpq_t *b, temp;
      b = xcalloc(1+n, sizeof(mpq_t));
      for (k = 1; k <= n; k++)
         mpq_init(b[k]), mpq_set(b[k], x[k]), mpq_set_si(x[k], 0, 1);
      mpq_init(temp);
      if (!tr)
      {  /* solve the system V*x = b */
         for (k = n; k >= 1; k--)
         {  i = P_row[k], j = Q_col[k];
            if (mpq_sgn(b[i]) != 0)
            {  mpq_set(x[j], b[i]);
               mpq_div(x[j], x[j], V_piv[i]);
               for (vij = V_col[j]; vij != NULL; vij = vij->c_next)
               {  mpq_mul(temp, vij->val, x[j]);
                  mpq_sub(b[vij->i], b[vij->i], temp);
               }
            }
         }
      }
      else
      {  /* solve the system V'*x = b */
         for (k = 1; k <= n; k++)
         {  i = P_row[k], j = Q_col[k];
            if (mpq_sgn(b[j]) != 0)
            {  mpq_set(x[i], b[j]);
               mpq_div(x[i], x[i], V_piv[i]);
               for (vij = V_row[i]; vij != NULL; vij = vij->r_next)
               {  mpq_mul(temp, vij->val, x[i]);
                  mpq_sub(b[vij->j], b[vij->j], temp);
               }
            }
         }
      }
      for (k = 1; k <= n; k++) mpq_clear(b[k]);
      mpq_clear(temp);
      xfree(b);
      return;
}

/***********************************************************************
*  lux_solve - solve system A*x = b or A'*x = b
*
*  SYNOPSIS
*
*  #include "lux.h"
*  void lux_solve(LUX *lux, int tr, mpq_t x[]);
*
*  DESCRIPTION
*
*  The routine lux_solve solves either the system A*x = b (if the flag
*  tr is zero) or the system A'*x = b (if the flag tr is non-zero),
*  where the parameter lux specifies LU-factorization of the matrix A,
*  A' is a matrix transposed to A.
*
*  On entry the array x should contain elements of the right-hand side
*  vector b in locations x[1], ..., x[n], where n is the order of the
*  matrix A. On exit this array will contain elements of the solution
*  vector x in the same locations. */

void lux_solve(LUX *lux, int tr, mpq_t x[])
{     if (lux->rank < lux->n)
         xfault("lux_solve: LU-factorization has incomplete rank\n");
      if (!tr)
      {  /* A = F*V, therefore inv(A) = inv(V)*inv(F) */
         lux_f_solve(lux, 0, x);
         lux_v_solve(lux, 0, x);
      }
      else
      {  /* A' = V'*F', therefore inv(A') = inv(F')*inv(V') */
         lux_v_solve(lux, 1, x);
         lux_f_solve(lux, 1, x);
      }
      return;
}

/***********************************************************************
*  lux_delete - delete LU-factorization
*
*  SYNOPSIS
*
*  #include "lux.h"
*  void lux_delete(LUX *lux);
*
*  DESCRIPTION
*
*  The routine lux_delete deletes LU-factorization data structure,
*  which the parameter lux points to, freeing all the memory allocated
*  to this object. */

void lux_delete(LUX *lux)
{     int n = lux->n;
      LUXELM *fij, *vij;
      int i;
      for (i = 1; i <= n; i++)
      {  for (fij = lux->F_row[i]; fij != NULL; fij = fij->r_next)
            mpq_clear(fij->val);
         mpq_clear(lux->V_piv[i]);
         for (vij = lux->V_row[i]; vij != NULL; vij = vij->r_next)
            mpq_clear(vij->val);
      }
      dmp_delete_pool(lux->pool);
      xfree(lux->F_row);
      xfree(lux->F_col);
      xfree(lux->V_piv);
      xfree(lux->V_row);
      xfree(lux->V_col);
      xfree(lux->P_row);
      xfree(lux->P_col);
      xfree(lux->Q_row);
      xfree(lux->Q_col);
      xfree(lux);
      return;
}

/* eof */