/* glplux.c */

 /***********************************************************************

 *  This code is part of GLPK (GNU Linear Programming Kit).

 *

 *  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,

 *  2009, 2010 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 "glplux.h"

 #define xfault xerror

 #define dmp_create_poolx(size) dmp_create_pool()

 /*----------------------------------------------------------------------

 // lux_create - create LU-factorization.

 //

 // SYNOPSIS

 //

 // #include "glplux.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 "glplux.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 "glplux.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 "glplux.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 "glplux.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 "glplux.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 */