praat/external/glpk/glpluf.c

/* glpluf.c (LU-factorization) */

/***********************************************************************
*  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 "glpenv.h"
#include "glpluf.h"
#define xfault xerror

/* CAUTION: DO NOT CHANGE THE LIMIT BELOW */

#define N_MAX 100000000 /* = 100*10^6 */
/* maximal order of the original matrix */

/***********************************************************************
*  NAME
*
*  luf_create_it - create LU-factorization
*
*  SYNOPSIS
*
*  #include "glpluf.h"
*  LUF *luf_create_it(void);
*
*  DESCRIPTION
*
*  The routine luf_create_it creates a program object, which represents
*  LU-factorization of a square matrix.
*
*  RETURNS
*
*  The routine luf_create_it returns a pointer to the object created. */

LUF *luf_create_it(void)
{     LUF *luf;
      luf = xmalloc(sizeof(LUF));
      luf->n_max = luf->n = 0;
      luf->valid = 0;
      luf->fr_ptr = luf->fr_len = NULL;
      luf->fc_ptr = luf->fc_len = NULL;
      luf->vr_ptr = luf->vr_len = luf->vr_cap = NULL;
      luf->vr_piv = NULL;
      luf->vc_ptr = luf->vc_len = luf->vc_cap = NULL;
      luf->pp_row = luf->pp_col = NULL;
      luf->qq_row = luf->qq_col = NULL;
      luf->sv_size = 0;
      luf->sv_beg = luf->sv_end = 0;
      luf->sv_ind = NULL;
      luf->sv_val = NULL;
      luf->sv_head = luf->sv_tail = 0;
      luf->sv_prev = luf->sv_next = NULL;
      luf->vr_max = NULL;
      luf->rs_head = luf->rs_prev = luf->rs_next = NULL;
      luf->cs_head = luf->cs_prev = luf->cs_next = NULL;
      luf->flag = NULL;
      luf->work = NULL;
      luf->new_sva = 0;
      luf->piv_tol = 0.10;
      luf->piv_lim = 4;
      luf->suhl = 1;
      luf->eps_tol = 1e-15;
      luf->max_gro = 1e+10;
      luf->nnz_a = luf->nnz_f = luf->nnz_v = 0;
      luf->max_a = luf->big_v = 0.0;
      luf->rank = 0;
      return luf;
}

/***********************************************************************
*  NAME
*
*  luf_defrag_sva - defragment the sparse vector area
*
*  SYNOPSIS
*
*  #include "glpluf.h"
*  void luf_defrag_sva(LUF *luf);
*
*  DESCRIPTION
*
*  The routine luf_defrag_sva defragments the sparse vector area (SVA)
*  gathering all unused locations in one continuous extent. In order to
*  do that the routine moves all unused locations from the left part of
*  SVA (which contains rows and columns of the matrix V) to the middle
*  part (which contains free locations). This is attained by relocating
*  elements of rows and columns of the matrix V toward the beginning of
*  the left part.
*
*  NOTE that this "garbage collection" involves changing row and column
*  pointers of the matrix V. */

void luf_defrag_sva(LUF *luf)
{     int n = luf->n;
      int *vr_ptr = luf->vr_ptr;
      int *vr_len = luf->vr_len;
      int *vr_cap = luf->vr_cap;
      int *vc_ptr = luf->vc_ptr;
      int *vc_len = luf->vc_len;
      int *vc_cap = luf->vc_cap;
      int *sv_ind = luf->sv_ind;
      double *sv_val = luf->sv_val;
      int *sv_next = luf->sv_next;
      int sv_beg = 1;
      int i, j, k;
      /* skip rows and columns, which do not need to be relocated */
      for (k = luf->sv_head; k != 0; k = sv_next[k])
      {  if (k <= n)
         {  /* i-th row of the matrix V */
            i = k;
            if (vr_ptr[i] != sv_beg) break;
            vr_cap[i] = vr_len[i];
            sv_beg += vr_cap[i];
         }
         else
         {  /* j-th column of the matrix V */
            j = k - n;
            if (vc_ptr[j] != sv_beg) break;
            vc_cap[j] = vc_len[j];
            sv_beg += vc_cap[j];
         }
      }
      /* relocate other rows and columns in order to gather all unused
         locations in one continuous extent */
      for (k = k; k != 0; k = sv_next[k])
      {  if (k <= n)
         {  /* i-th row of the matrix V */
            i = k;
            memmove(&sv_ind[sv_beg], &sv_ind[vr_ptr[i]],
               vr_len[i] * sizeof(int));
            memmove(&sv_val[sv_beg], &sv_val[vr_ptr[i]],
               vr_len[i] * sizeof(double));
            vr_ptr[i] = sv_beg;
            vr_cap[i] = vr_len[i];
            sv_beg += vr_cap[i];
         }
         else
         {  /* j-th column of the matrix V */
            j = k - n;
            memmove(&sv_ind[sv_beg], &sv_ind[vc_ptr[j]],
               vc_len[j] * sizeof(int));
            memmove(&sv_val[sv_beg], &sv_val[vc_ptr[j]],
               vc_len[j] * sizeof(double));
            vc_ptr[j] = sv_beg;
            vc_cap[j] = vc_len[j];
            sv_beg += vc_cap[j];
         }
      }
      /* set new pointer to the beginning of the free part */
      luf->sv_beg = sv_beg;
      return;
}

/***********************************************************************
*  NAME
*
*  luf_enlarge_row - enlarge row capacity
*
*  SYNOPSIS
*
*  #include "glpluf.h"
*  int luf_enlarge_row(LUF *luf, int i, int cap);
*
*  DESCRIPTION
*
*  The routine luf_enlarge_row enlarges capacity of the i-th row of the
*  matrix V to cap locations (assuming that its current capacity is less
*  than cap). In order to do that the routine relocates elements of the
*  i-th row to the end of the left part of SVA (which contains rows and
*  columns of the matrix V) and then expands the left part by allocating
*  cap free locations from the free part. If there are less than cap
*  free locations, the routine defragments the sparse vector area.
*
*  Due to "garbage collection" this operation may change row and column
*  pointers of the matrix V.
*
*  RETURNS
*
*  If no error occured, the routine returns zero. Otherwise, in case of
*  overflow of the sparse vector area, the routine returns non-zero. */

int luf_enlarge_row(LUF *luf, int i, int cap)
{     int n = luf->n;
      int *vr_ptr = luf->vr_ptr;
      int *vr_len = luf->vr_len;
      int *vr_cap = luf->vr_cap;
      int *vc_cap = luf->vc_cap;
      int *sv_ind = luf->sv_ind;
      double *sv_val = luf->sv_val;
      int *sv_prev = luf->sv_prev;
      int *sv_next = luf->sv_next;
      int ret = 0;
      int cur, k, kk;
      xassert(1 <= i && i <= n);
      xassert(vr_cap[i] < cap);
      /* if there are less than cap free locations, defragment SVA */
      if (luf->sv_end - luf->sv_beg < cap)
      {  luf_defrag_sva(luf);
         if (luf->sv_end - luf->sv_beg < cap)
         {  ret = 1;
            goto done;
         }
      }
      /* save current capacity of the i-th row */
      cur = vr_cap[i];
      /* copy existing elements to the beginning of the free part */
      memmove(&sv_ind[luf->sv_beg], &sv_ind[vr_ptr[i]],
         vr_len[i] * sizeof(int));
      memmove(&sv_val[luf->sv_beg], &sv_val[vr_ptr[i]],
         vr_len[i] * sizeof(double));
      /* set new pointer and new capacity of the i-th row */
      vr_ptr[i] = luf->sv_beg;
      vr_cap[i] = cap;
      /* set new pointer to the beginning of the free part */
      luf->sv_beg += cap;
      /* now the i-th row starts in the rightmost location among other
         rows and columns of the matrix V, so its node should be moved
         to the end of the row/column linked list */
      k = i;
      /* remove the i-th row node from the linked list */
      if (sv_prev[k] == 0)
         luf->sv_head = sv_next[k];
      else
      {  /* capacity of the previous row/column can be increased at the
            expense of old locations of the i-th row */
         kk = sv_prev[k];
         if (kk <= n) vr_cap[kk] += cur; else vc_cap[kk-n] += cur;
         sv_next[sv_prev[k]] = sv_next[k];
      }
      if (sv_next[k] == 0)
         luf->sv_tail = sv_prev[k];
      else
         sv_prev[sv_next[k]] = sv_prev[k];
      /* insert the i-th row node to the end of the linked list */
      sv_prev[k] = luf->sv_tail;
      sv_next[k] = 0;
      if (sv_prev[k] == 0)
         luf->sv_head = k;
      else
         sv_next[sv_prev[k]] = k;
      luf->sv_tail = k;
done: return ret;
}

/***********************************************************************
*  NAME
*
*  luf_enlarge_col - enlarge column capacity
*
*  SYNOPSIS
*
*  #include "glpluf.h"
*  int luf_enlarge_col(LUF *luf, int j, int cap);
*
*  DESCRIPTION
*
*  The routine luf_enlarge_col enlarges capacity of the j-th column of
*  the matrix V to cap locations (assuming that its current capacity is
*  less than cap). In order to do that the routine relocates elements
*  of the j-th column to the end of the left part of SVA (which contains
*  rows and columns of the matrix V) and then expands the left part by
*  allocating cap free locations from the free part. If there are less
*  than cap free locations, the routine defragments the sparse vector
*  area.
*
*  Due to "garbage collection" this operation may change row and column
*  pointers of the matrix V.
*
*  RETURNS
*
*  If no error occured, the routine returns zero. Otherwise, in case of
*  overflow of the sparse vector area, the routine returns non-zero. */

int luf_enlarge_col(LUF *luf, int j, int cap)
{     int n = luf->n;
      int *vr_cap = luf->vr_cap;
      int *vc_ptr = luf->vc_ptr;
      int *vc_len = luf->vc_len;
      int *vc_cap = luf->vc_cap;
      int *sv_ind = luf->sv_ind;
      double *sv_val = luf->sv_val;
      int *sv_prev = luf->sv_prev;
      int *sv_next = luf->sv_next;
      int ret = 0;
      int cur, k, kk;
      xassert(1 <= j && j <= n);
      xassert(vc_cap[j] < cap);
      /* if there are less than cap free locations, defragment SVA */
      if (luf->sv_end - luf->sv_beg < cap)
      {  luf_defrag_sva(luf);
         if (luf->sv_end - luf->sv_beg < cap)
         {  ret = 1;
            goto done;
         }
      }
      /* save current capacity of the j-th column */
      cur = vc_cap[j];
      /* copy existing elements to the beginning of the free part */
      memmove(&sv_ind[luf->sv_beg], &sv_ind[vc_ptr[j]],
         vc_len[j] * sizeof(int));
      memmove(&sv_val[luf->sv_beg], &sv_val[vc_ptr[j]],
         vc_len[j] * sizeof(double));
      /* set new pointer and new capacity of the j-th column */
      vc_ptr[j] = luf->sv_beg;
      vc_cap[j] = cap;
      /* set new pointer to the beginning of the free part */
      luf->sv_beg += cap;
      /* now the j-th column starts in the rightmost location among
         other rows and columns of the matrix V, so its node should be
         moved to the end of the row/column linked list */
      k = n + j;
      /* remove the j-th column node from the linked list */
      if (sv_prev[k] == 0)
         luf->sv_head = sv_next[k];
      else
      {  /* capacity of the previous row/column can be increased at the
            expense of old locations of the j-th column */
         kk = sv_prev[k];
         if (kk <= n) vr_cap[kk] += cur; else vc_cap[kk-n] += cur;
         sv_next[sv_prev[k]] = sv_next[k];
      }
      if (sv_next[k] == 0)
         luf->sv_tail = sv_prev[k];
      else
         sv_prev[sv_next[k]] = sv_prev[k];
      /* insert the j-th column node to the end of the linked list */
      sv_prev[k] = luf->sv_tail;
      sv_next[k] = 0;
      if (sv_prev[k] == 0)
         luf->sv_head = k;
      else
         sv_next[sv_prev[k]] = k;
      luf->sv_tail = k;
done: return ret;
}

/***********************************************************************
*  reallocate - reallocate LU-factorization arrays
*
*  This routine reallocates arrays, whose size depends of n, the order
*  of the matrix A to be factorized. */

static void reallocate(LUF *luf, int n)
{     int n_max = luf->n_max;
      luf->n = n;
      if (n <= n_max) goto done;
      if (luf->fr_ptr != NULL) xfree(luf->fr_ptr);
      if (luf->fr_len != NULL) xfree(luf->fr_len);
      if (luf->fc_ptr != NULL) xfree(luf->fc_ptr);
      if (luf->fc_len != NULL) xfree(luf->fc_len);
      if (luf->vr_ptr != NULL) xfree(luf->vr_ptr);
      if (luf->vr_len != NULL) xfree(luf->vr_len);
      if (luf->vr_cap != NULL) xfree(luf->vr_cap);
      if (luf->vr_piv != NULL) xfree(luf->vr_piv);
      if (luf->vc_ptr != NULL) xfree(luf->vc_ptr);
      if (luf->vc_len != NULL) xfree(luf->vc_len);
      if (luf->vc_cap != NULL) xfree(luf->vc_cap);
      if (luf->pp_row != NULL) xfree(luf->pp_row);
      if (luf->pp_col != NULL) xfree(luf->pp_col);
      if (luf->qq_row != NULL) xfree(luf->qq_row);
      if (luf->qq_col != NULL) xfree(luf->qq_col);
      if (luf->sv_prev != NULL) xfree(luf->sv_prev);
      if (luf->sv_next != NULL) xfree(luf->sv_next);
      if (luf->vr_max != NULL) xfree(luf->vr_max);
      if (luf->rs_head != NULL) xfree(luf->rs_head);
      if (luf->rs_prev != NULL) xfree(luf->rs_prev);
      if (luf->rs_next != NULL) xfree(luf->rs_next);
      if (luf->cs_head != NULL) xfree(luf->cs_head);
      if (luf->cs_prev != NULL) xfree(luf->cs_prev);
      if (luf->cs_next != NULL) xfree(luf->cs_next);
      if (luf->flag != NULL) xfree(luf->flag);
      if (luf->work != NULL) xfree(luf->work);
      luf->n_max = n_max = n + 100;
      luf->fr_ptr = xcalloc(1+n_max, sizeof(int));
      luf->fr_len = xcalloc(1+n_max, sizeof(int));
      luf->fc_ptr = xcalloc(1+n_max, sizeof(int));
      luf->fc_len = xcalloc(1+n_max, sizeof(int));
      luf->vr_ptr = xcalloc(1+n_max, sizeof(int));
      luf->vr_len = xcalloc(1+n_max, sizeof(int));
      luf->vr_cap = xcalloc(1+n_max, sizeof(int));
      luf->vr_piv = xcalloc(1+n_max, sizeof(double));
      luf->vc_ptr = xcalloc(1+n_max, sizeof(int));
      luf->vc_len = xcalloc(1+n_max, sizeof(int));
      luf->vc_cap = xcalloc(1+n_max, sizeof(int));
      luf->pp_row = xcalloc(1+n_max, sizeof(int));
      luf->pp_col = xcalloc(1+n_max, sizeof(int));
      luf->qq_row = xcalloc(1+n_max, sizeof(int));
      luf->qq_col = xcalloc(1+n_max, sizeof(int));
      luf->sv_prev = xcalloc(1+n_max+n_max, sizeof(int));
      luf->sv_next = xcalloc(1+n_max+n_max, sizeof(int));
      luf->vr_max = xcalloc(1+n_max, sizeof(double));
      luf->rs_head = xcalloc(1+n_max, sizeof(int));
      luf->rs_prev = xcalloc(1+n_max, sizeof(int));
      luf->rs_next = xcalloc(1+n_max, sizeof(int));
      luf->cs_head = xcalloc(1+n_max, sizeof(int));
      luf->cs_prev = xcalloc(1+n_max, sizeof(int));
      luf->cs_next = xcalloc(1+n_max, sizeof(int));
      luf->flag = xcalloc(1+n_max, sizeof(int));
      luf->work = xcalloc(1+n_max, sizeof(double));
done: return;
}

/***********************************************************************
*  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. (Row-wise representation of the matrix F is not used
*  at the factorization stage and therefore is not initialized.)
*
*  If no error occured, the routine returns zero. Otherwise, in case of
*  overflow of the sparse vector area, the routine returns non-zero. */

static int initialize(LUF *luf, int (*col)(void *info, int j, int rn[],
      double aj[]), void *info)
{     int n = luf->n;
      int *fc_ptr = luf->fc_ptr;
      int *fc_len = luf->fc_len;
      int *vr_ptr = luf->vr_ptr;
      int *vr_len = luf->vr_len;
      int *vr_cap = luf->vr_cap;
      int *vc_ptr = luf->vc_ptr;
      int *vc_len = luf->vc_len;
      int *vc_cap = luf->vc_cap;
      int *pp_row = luf->pp_row;
      int *pp_col = luf->pp_col;
      int *qq_row = luf->qq_row;
      int *qq_col = luf->qq_col;
      int *sv_ind = luf->sv_ind;
      double *sv_val = luf->sv_val;
      int *sv_prev = luf->sv_prev;
      int *sv_next = luf->sv_next;
      double *vr_max = luf->vr_max;
      int *rs_head = luf->rs_head;
      int *rs_prev = luf->rs_prev;
      int *rs_next = luf->rs_next;
      int *cs_head = luf->cs_head;
      int *cs_prev = luf->cs_prev;
      int *cs_next = luf->cs_next;
      int *flag = luf->flag;
      double *work = luf->work;
      int ret = 0;
      int i, i_ptr, j, j_beg, j_end, k, len, nnz, sv_beg, sv_end, ptr;
      double big, val;
      /* free all locations of the sparse vector area */
      sv_beg = 1;
      sv_end = luf->sv_size + 1;
      /* (row-wise representation of the matrix F is not initialized,
         because it is not used at the factorization stage) */
      /* build the matrix F in column-wise format (initially F = I) */
      for (j = 1; j <= n; j++)
      {  fc_ptr[j] = sv_end;
         fc_len[j] = 0;
      }
      /* clear rows of the matrix V; clear the flag array */
      for (i = 1; i <= n; i++)
         vr_len[i] = vr_cap[i] = 0, flag[i] = 0;
      /* build the matrix V in column-wise format (initially V = A);
         count non-zeros in rows of this matrix; count total number of
         non-zeros; compute largest of absolute values of elements */
      nnz = 0;
      big = 0.0;
      for (j = 1; j <= n; j++)
      {  int *rn = pp_row;
         double *aj = work;
         /* obtain j-th column of the matrix A */
         len = col(info, j, rn, aj);
         if (!(0 <= len && len <= n))
            xfault("luf_factorize: j = %d; len = %d; invalid column len"
               "gth\n", j, len);
         /* check for free locations */
         if (sv_end - sv_beg < len)
         {  /* overflow of the sparse vector area */
            ret = 1;
            goto done;
         }
         /* set pointer to the j-th column */
         vc_ptr[j] = sv_beg;
         /* set length of the j-th column */
         vc_len[j] = vc_cap[j] = len;
         /* count total number of non-zeros */
         nnz += len;
         /* walk through elements of the j-th column */
         for (ptr = 1; ptr <= len; ptr++)
         {  /* get row index and numerical value of a[i,j] */
            i = rn[ptr];
            val = aj[ptr];
            if (!(1 <= i && i <= n))
               xfault("luf_factorize: i = %d; j = %d; invalid row index"
                  "\n", i, j);
            if (flag[i])
               xfault("luf_factorize: i = %d; j = %d; duplicate element"
                  " not allowed\n", i, j);
            if (val == 0.0)
               xfault("luf_factorize: i = %d; j = %d; zero element not "
                  "allowed\n", i, j);
            /* add new element v[i,j] = a[i,j] to j-th column */
            sv_ind[sv_beg] = i;
            sv_val[sv_beg] = val;
            sv_beg++;
            /* big := max(big, |a[i,j]|) */
            if (val < 0.0) val = - val;
            if (big < val) big = val;
            /* mark non-zero in the i-th position of the j-th column */
            flag[i] = 1;
            /* increase length of the i-th row */
            vr_cap[i]++;
         }
         /* reset all non-zero marks */
         for (ptr = 1; ptr <= len; ptr++) flag[rn[ptr]] = 0;
      }
      /* allocate rows of the matrix V */
      for (i = 1; i <= n; i++)
      {  /* get length of the i-th row */
         len = vr_cap[i];
         /* check for free locations */
         if (sv_end - sv_beg < len)
         {  /* overflow of the sparse vector area */
            ret = 1;
            goto done;
         }
         /* set pointer to the i-th row */
         vr_ptr[i] = sv_beg;
         /* reserve locations for the i-th row */
         sv_beg += len;
      }
      /* build the matrix V in row-wise format using representation of
         this matrix in column-wise format */
      for (j = 1; j <= n; j++)
      {  /* walk through elements of the j-th column */
         j_beg = vc_ptr[j];
         j_end = j_beg + vc_len[j] - 1;
         for (k = j_beg; k <= j_end; k++)
         {  /* get row index and numerical value of v[i,j] */
            i = sv_ind[k];
            val = sv_val[k];
            /* store element in the i-th row */
            i_ptr = vr_ptr[i] + vr_len[i];
            sv_ind[i_ptr] = j;
            sv_val[i_ptr] = val;
            /* increase count of the i-th row */
            vr_len[i]++;
         }
      }
      /* initialize the matrices P and Q (initially P = Q = I) */
      for (k = 1; k <= n; k++)
         pp_row[k] = pp_col[k] = qq_row[k] = qq_col[k] = k;
      /* set sva partitioning pointers */
      luf->sv_beg = sv_beg;
      luf->sv_end = sv_end;
      /* the initial physical order of rows and columns of the matrix V
         is n+1, ..., n+n, 1, ..., n (firstly columns, then rows) */
      luf->sv_head = n+1;
      luf->sv_tail = n;
      for (i = 1; i <= n; i++)
      {  sv_prev[i] = i-1;
         sv_next[i] = i+1;
      }
      sv_prev[1] = n+n;
      sv_next[n] = 0;
      for (j = 1; j <= n; j++)
      {  sv_prev[n+j] = n+j-1;
         sv_next[n+j] = n+j+1;
      }
      sv_prev[n+1] = 0;
      sv_next[n+n] = 1;
      /* clear working arrays */
      for (k = 1; k <= n; k++)
      {  flag[k] = 0;
         work[k] = 0.0;
      }
      /* initialize some statistics */
      luf->nnz_a = nnz;
      luf->nnz_f = 0;
      luf->nnz_v = nnz;
      luf->max_a = big;
      luf->big_v = big;
      luf->rank = -1;
      /* initially the active submatrix is the entire matrix V */
      /* largest of absolute values of elements in each active row is
         unknown yet */
      for (i = 1; i <= n; i++) vr_max[i] = -1.0;
      /* build linked lists of active rows */
      for (len = 0; len <= n; len++) rs_head[len] = 0;
      for (i = 1; i <= n; i++)
      {  len = vr_len[i];
         rs_prev[i] = 0;
         rs_next[i] = rs_head[len];
         if (rs_next[i] != 0) rs_prev[rs_next[i]] = i;
         rs_head[len] = i;
      }
      /* build linked lists of active columns */
      for (len = 0; len <= n; len++) cs_head[len] = 0;
      for (j = 1; j <= n; j++)
      {  len = vc_len[j];
         cs_prev[j] = 0;
         cs_next[j] = cs_head[len];
         if (cs_next[j] != 0) cs_prev[cs_next[j]] = j;
         cs_head[len] = j;
      }
done: /* return to the factorizing routine */
      return ret;
}

/***********************************************************************
*  find_pivot - choose a pivot element
*
*  This routine chooses a pivot element in the active submatrix of the
*  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 doesn't 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).
* 
*  For the reason of numerical stability the routine applies so called
*  threshold pivoting proposed by J.Reid. It is assumed that an element
*  v[i,j] can be selected as a pivot candidate if it is not very small
*  (in absolute value) among other elements in the same row, i.e. if it
*  satisfies to the stability condition |v[i,j]| >= tol * max|v[i,*]|,
*  where 0 < tol < 1 is a given tolerance.
* 
*  In order to keep sparsity of the matrix V the routine uses Markowitz
*  strategy, trying to choose such element v[p,q], which satisfies to
*  the stability condition (see above) and has smallest Markowitz cost
*  (nr[p]-1) * (nc[q]-1), where nr[p] and nc[q] are numbers of non-zero
*  elements, respectively, in the p-th row and in the 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.
* 
*  If the pivot element v[p,q] has been chosen, the routine stores its
*  indices to the locations *p and *q and returns zero. Otherwise, if
*  the active submatrix is empty and therefore the pivot element can't
*  be chosen, the routine returns non-zero. */

static int find_pivot(LUF *luf, int *_p, int *_q)
{     int n = luf->n;
      int *vr_ptr = luf->vr_ptr;
      int *vr_len = luf->vr_len;
      int *vc_ptr = luf->vc_ptr;
      int *vc_len = luf->vc_len;
      int *sv_ind = luf->sv_ind;
      double *sv_val = luf->sv_val;
      double *vr_max = luf->vr_max;
      int *rs_head = luf->rs_head;
      int *rs_next = luf->rs_next;
      int *cs_head = luf->cs_head;
      int *cs_prev = luf->cs_prev;
      int *cs_next = luf->cs_next;
      double piv_tol = luf->piv_tol;
      int piv_lim = luf->piv_lim;
      int suhl = luf->suhl;
      int p, q, len, i, i_beg, i_end, i_ptr, j, j_beg, j_end, j_ptr,
         ncand, next_j, min_p, min_q, min_len;
      double best, cost, big, temp;
      /* initially no pivot candidates have been found so far */
      p = q = 0, 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 pivot */
      j = cs_head[1];
      if (j != 0)
      {  xassert(vc_len[j] == 1);
         p = sv_ind[vc_ptr[j]], q = j;
         goto done;
      }
      /* if in the active submatrix there is a row that has the only
         non-zero (row singleton), choose it as pivot */
      i = rs_head[1];
      if (i != 0)
      {  xassert(vr_len[i] == 1);
         p = i, q = sv_ind[vr_ptr[i]];
         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 that have len non-zeros */
         for (j = cs_head[len]; j != 0; j = next_j)
         {  /* the j-th column has len non-zeros */
            j_beg = vc_ptr[j];
            j_end = j_beg + vc_len[j] - 1;
            /* save pointer to the next column with the same length */
            next_j = cs_next[j];
            /* find an element in the j-th column, which is placed in a
               row with minimal number of non-zeros and satisfies to the
               stability condition (such element may not exist) */
            min_p = min_q = 0, min_len = INT_MAX;
            for (j_ptr = j_beg; j_ptr <= j_end; j_ptr++)
            {  /* get row index of v[i,j] */
               i = sv_ind[j_ptr];
               i_beg = vr_ptr[i];
               i_end = i_beg + vr_len[i] - 1;
               /* if the i-th row is not shorter than that one, where
                  minimal element is currently placed, skip v[i,j] */
               if (vr_len[i] >= min_len) continue;
               /* determine the largest of absolute values of elements
                  in the i-th row */
               big = vr_max[i];
               if (big < 0.0)
               {  /* the largest value is unknown yet; compute it */
                  for (i_ptr = i_beg; i_ptr <= i_end; i_ptr++)
                  {  temp = sv_val[i_ptr];
                     if (temp < 0.0) temp = - temp;
                     if (big < temp) big = temp;
                  }
                  vr_max[i] = big;
               }
               /* find v[i,j] in the i-th row */
               for (i_ptr = vr_ptr[i]; sv_ind[i_ptr] != j; i_ptr++);
               xassert(i_ptr <= i_end);
               /* if v[i,j] doesn't satisfy to the stability condition,
                  skip it */
               temp = sv_val[i_ptr];
               if (temp < 0.0) temp = - temp;
               if (temp < piv_tol * big) continue;
               /* v[i,j] is better than the current minimal element */
               min_p = i, min_q = j, min_len = vr_len[i];
               /* if Markowitz cost of the current minimal 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)
               {  p = min_p, q = min_q;
                  goto done;
               }
            }
            /* the j-th column has been scanned */
            if (min_p != 0)
            {  /* the minimal element is a next pivot candidate */
               ncand++;
               /* compute its Markowitz cost */
               cost = (double)(min_len - 1) * (double)(len - 1);
               /* choose between the minimal element and the current
                  candidate */
               if (cost < best) p = min_p, q = min_q, best = cost;
               /* if piv_lim candidates have been considered, there are
                  doubts that a much better candidate exists; therefore
                  it's time to terminate the search */
               if (ncand == piv_lim) goto done;
            }
            else
            {  /* the j-th column has no elements, which satisfy to the
                  stability condition; Uwe Suhl suggests to exclude such
                  column from the further consideration until it becomes
                  a column singleton; in hard cases this significantly
                  reduces a time needed for pivot searching */
               if (suhl)
               {  /* remove the j-th column from the active set */
                  if (cs_prev[j] == 0)
                     cs_head[len] = cs_next[j];
                  else
                     cs_next[cs_prev[j]] = cs_next[j];
                  if (cs_next[j] == 0)
                     /* nop */;
                  else
                     cs_prev[cs_next[j]] = cs_prev[j];
                  /* the following assignment is used to avoid an error
                     when the routine eliminate (see below) will try to
                     remove the j-th column from the active set */
                  cs_prev[j] = cs_next[j] = j;
               }
            }
         }
         /* consider active rows that have len non-zeros */
         for (i = rs_head[len]; i != 0; i = rs_next[i])
         {  /* the i-th row has len non-zeros */
            i_beg = vr_ptr[i];
            i_end = i_beg + vr_len[i] - 1;
            /* determine the largest of absolute values of elements in
               the i-th row */
            big = vr_max[i];
            if (big < 0.0)
            {  /* the largest value is unknown yet; compute it */
               for (i_ptr = i_beg; i_ptr <= i_end; i_ptr++)
               {  temp = sv_val[i_ptr];
                  if (temp < 0.0) temp = - temp;
                  if (big < temp) big = temp;
               }
               vr_max[i] = big;
            }
            /* find an element in the i-th row, which is placed in a
               column with minimal number of non-zeros and satisfies to
               the stability condition (such element always exists) */
            min_p = min_q = 0, min_len = INT_MAX;
            for (i_ptr = i_beg; i_ptr <= i_end; i_ptr++)
            {  /* get column index of v[i,j] */
               j = sv_ind[i_ptr];
               /* if the j-th column is not shorter than that one, where
                  minimal element is currently placed, skip v[i,j] */
               if (vc_len[j] >= min_len) continue;
               /* if v[i,j] doesn't satisfy to the stability condition,
                  skip it */
               temp = sv_val[i_ptr];
               if (temp < 0.0) temp = - temp;
               if (temp < piv_tol * big) continue;
               /* v[i,j] is better than the current minimal element */
               min_p = i, min_q = j, min_len = vc_len[j];
               /* if Markowitz cost of the current minimal 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)
               {  p = min_p, q = min_q;
                  goto done;
               }
            }
            /* the i-th row has been scanned */
            if (min_p != 0)
            {  /* the minimal element is a next pivot candidate */
               ncand++;
               /* compute its Markowitz cost */
               cost = (double)(len - 1) * (double)(min_len - 1);
               /* choose between the minimal element and the current
                  candidate */
               if (cost < best) p = min_p, q = min_q, best = cost;
               /* if piv_lim candidates have been considered, there are
                  doubts that a much better candidate exists; therefore
                  it's time to terminate the search */
               if (ncand == piv_lim) goto done;
            }
            else
            {  /* this can't be because this can never be */
               xassert(min_p != min_p);
            }
         }
      }
done: /* bring the pivot to the factorizing routine */
      *_p = p, *_q = q;
      return (p == 0);
}

/***********************************************************************
*  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 parameters p and q are, respectively, row and column indices of
*  the element v[p,q], which corresponds to the element 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.
*
*  If no error occured, the routine returns zero. Otherwise, in case of
*  overflow of the sparse vector area, the routine returns non-zero. */

static int eliminate(LUF *luf, int p, int q)
{     int n = luf->n;
      int *fc_ptr = luf->fc_ptr;
      int *fc_len = luf->fc_len;
      int *vr_ptr = luf->vr_ptr;
      int *vr_len = luf->vr_len;
      int *vr_cap = luf->vr_cap;
      double *vr_piv = luf->vr_piv;
      int *vc_ptr = luf->vc_ptr;
      int *vc_len = luf->vc_len;
      int *vc_cap = luf->vc_cap;
      int *sv_ind = luf->sv_ind;
      double *sv_val = luf->sv_val;
      int *sv_prev = luf->sv_prev;
      int *sv_next = luf->sv_next;
      double *vr_max = luf->vr_max;
      int *rs_head = luf->rs_head;
      int *rs_prev = luf->rs_prev;
      int *rs_next = luf->rs_next;
      int *cs_head = luf->cs_head;
      int *cs_prev = luf->cs_prev;
      int *cs_next = luf->cs_next;
      int *flag = luf->flag;
      double *work = luf->work;
      double eps_tol = luf->eps_tol;
      /* at this stage the row-wise representation of the matrix F is
         not used, so fr_len can be used as a working array */
      int *ndx = luf->fr_len;
      int ret = 0;
      int len, fill, i, i_beg, i_end, i_ptr, j, j_beg, j_end, j_ptr, k,
         p_beg, p_end, p_ptr, q_beg, q_end, q_ptr;
      double fip, val, vpq, temp;
      xassert(1 <= p && p <= n);
      xassert(1 <= q && q <= n);
      /* remove the p-th (pivot) row from the active set; this row will
         never return there */
      if (rs_prev[p] == 0)
         rs_head[vr_len[p]] = rs_next[p];
      else
         rs_next[rs_prev[p]] = rs_next[p];
      if (rs_next[p] == 0)
         ;
      else
         rs_prev[rs_next[p]] = rs_prev[p];
      /* remove the q-th (pivot) column from the active set; this column
         will never return there */
      if (cs_prev[q] == 0)
         cs_head[vc_len[q]] = cs_next[q];
      else
         cs_next[cs_prev[q]] = cs_next[q];
      if (cs_next[q] == 0)
         ;
      else
         cs_prev[cs_next[q]] = cs_prev[q];
      /* find the pivot v[p,q] = u[k,k] in the p-th row */
      p_beg = vr_ptr[p];
      p_end = p_beg + vr_len[p] - 1;
      for (p_ptr = p_beg; sv_ind[p_ptr] != q; p_ptr++) /* nop */;
      xassert(p_ptr <= p_end);
      /* store value of the pivot */
      vpq = (vr_piv[p] = sv_val[p_ptr]);
      /* remove the pivot from the p-th row */
      sv_ind[p_ptr] = sv_ind[p_end];
      sv_val[p_ptr] = sv_val[p_end];
      vr_len[p]--;
      p_end--;
      /* find the pivot v[p,q] = u[k,k] in the q-th column */
      q_beg = vc_ptr[q];
      q_end = q_beg + vc_len[q] - 1;
      for (q_ptr = q_beg; sv_ind[q_ptr] != p; q_ptr++) /* nop */;
      xassert(q_ptr <= q_end);
      /* remove the pivot from the q-th column */
      sv_ind[q_ptr] = sv_ind[q_end];
      vc_len[q]--;
      q_end--;
      /* walk through the p-th (pivot) row, which doesn't contain the
         pivot v[p,q] already, and do the following... */
      for (p_ptr = p_beg; p_ptr <= p_end; p_ptr++)
      {  /* get column index of v[p,j] */
         j = sv_ind[p_ptr];
         /* store v[p,j] to the working array */
         flag[j] = 1;
         work[j] = sv_val[p_ptr];
         /* remove the j-th column from the active set; this column will
            return there later with new length */
         if (cs_prev[j] == 0)
            cs_head[vc_len[j]] = cs_next[j];
         else
            cs_next[cs_prev[j]] = cs_next[j];
         if (cs_next[j] == 0)
            ;
         else
            cs_prev[cs_next[j]] = cs_prev[j];
         /* find v[p,j] in the j-th column */
         j_beg = vc_ptr[j];
         j_end = j_beg + vc_len[j] - 1;
         for (j_ptr = j_beg; sv_ind[j_ptr] != p; j_ptr++) /* nop */;
         xassert(j_ptr <= j_end);
         /* since v[p,j] leaves the active submatrix, remove it from the
            j-th column; however, v[p,j] is kept in the p-th row */
         sv_ind[j_ptr] = sv_ind[j_end];
         vc_len[j]--;
      }
      /* walk through the q-th (pivot) column, which doesn't contain the
         pivot v[p,q] already, and perform gaussian elimination */
      while (q_beg <= q_end)
      {  /* element v[i,q] should be eliminated */
         /* get row index of v[i,q] */
         i = sv_ind[q_beg];
         /* remove the i-th row from the active set; later this row will
            return there with new length */
         if (rs_prev[i] == 0)
            rs_head[vr_len[i]] = rs_next[i];
         else
            rs_next[rs_prev[i]] = rs_next[i];
         if (rs_next[i] == 0)
            ;
         else
            rs_prev[rs_next[i]] = rs_prev[i];
         /* find v[i,q] in the i-th row */
         i_beg = vr_ptr[i];
         i_end = i_beg + vr_len[i] - 1;
         for (i_ptr = i_beg; sv_ind[i_ptr] != q; i_ptr++) /* nop */;
         xassert(i_ptr <= i_end);
         /* compute gaussian multiplier f[i,p] = v[i,q] / v[p,q] */
         fip = sv_val[i_ptr] / vpq;
         /* since v[i,q] should be eliminated, remove it from the i-th
            row */
         sv_ind[i_ptr] = sv_ind[i_end];
         sv_val[i_ptr] = sv_val[i_end];
         vr_len[i]--;
         i_end--;
         /* and from the q-th column */
         sv_ind[q_beg] = sv_ind[q_end];
         vc_len[q]--;
         q_end--;
         /* perform gaussian transformation:
            (i-th row) := (i-th row) - f[i,p] * (p-th row)
            note that now the p-th row, which is in the working array,
            doesn't contain the pivot v[p,q], and the i-th row doesn't
            contain the eliminated element v[i,q] */
         /* walk through the i-th row and transform existing non-zero
            elements */
         fill = vr_len[p];
         for (i_ptr = i_beg; i_ptr <= i_end; i_ptr++)
         {  /* get column index of v[i,j] */
            j = sv_ind[i_ptr];
            /* v[i,j] := v[i,j] - f[i,p] * v[p,j] */
            if (flag[j])
            {  /* v[p,j] != 0 */
               temp = (sv_val[i_ptr] -= fip * work[j]);
               if (temp < 0.0) temp = - temp;
               flag[j] = 0;
               fill--; /* since both v[i,j] and v[p,j] exist */
               if (temp == 0.0 || temp < eps_tol)
               {  /* new v[i,j] is closer to zero; replace it by exact
                     zero, i.e. remove it from the active submatrix */
                  /* remove v[i,j] from the i-th row */
                  sv_ind[i_ptr] = sv_ind[i_end];
                  sv_val[i_ptr] = sv_val[i_end];
                  vr_len[i]--;
                  i_ptr--;
                  i_end--;
                  /* find v[i,j] in the j-th column */
                  j_beg = vc_ptr[j];
                  j_end = j_beg + vc_len[j] - 1;
                  for (j_ptr = j_beg; sv_ind[j_ptr] != i; j_ptr++);
                  xassert(j_ptr <= j_end);
                  /* remove v[i,j] from the j-th column */
                  sv_ind[j_ptr] = sv_ind[j_end];
                  vc_len[j]--;
               }
               else
               {  /* v_big := max(v_big, |v[i,j]|) */
                  if (luf->big_v < temp) luf->big_v = temp;
               }
            }
         }
         /* now flag is the pattern of the set v[p,*] \ v[i,*], and fill
            is number of non-zeros in this set; therefore up to fill new
            non-zeros may appear in the i-th row */
         if (vr_len[i] + fill > vr_cap[i])
         {  /* enlarge the i-th row */
            if (luf_enlarge_row(luf, i, vr_len[i] + fill))
            {  /* overflow of the sparse vector area */
               ret = 1;
               goto done;
            }
            /* defragmentation may change row and column pointers of the
               matrix V */
            p_beg = vr_ptr[p];
            p_end = p_beg + vr_len[p] - 1;
            q_beg = vc_ptr[q];
            q_end = q_beg + vc_len[q] - 1;
         }
         /* walk through the p-th (pivot) row and create new elements
            of the i-th row that appear due to fill-in; column indices
            of these new elements are accumulated in the array ndx */
         len = 0;
         for (p_ptr = p_beg; p_ptr <= p_end; p_ptr++)
         {  /* get column index of v[p,j], which may cause fill-in */
            j = sv_ind[p_ptr];
            if (flag[j])
            {  /* compute new non-zero v[i,j] = 0 - f[i,p] * v[p,j] */
               temp = (val = - fip * work[j]);
               if (temp < 0.0) temp = - temp;
               if (temp == 0.0 || temp < eps_tol)
                  /* if v[i,j] is closer to zero; just ignore it */;
               else
               {  /* add v[i,j] to the i-th row */
                  i_ptr = vr_ptr[i] + vr_len[i];
                  sv_ind[i_ptr] = j;
                  sv_val[i_ptr] = val;
                  vr_len[i]++;
                  /* remember column index of v[i,j] */
                  ndx[++len] = j;
                  /* big_v := max(big_v, |v[i,j]|) */
                  if (luf->big_v < temp) luf->big_v = temp;
               }
            }
            else
            {  /* there is no fill-in, because v[i,j] already exists in
                  the i-th row; restore the flag of the element v[p,j],
                  which was reset before */
               flag[j] = 1;
            }
         }
         /* add new non-zeros v[i,j] to the corresponding columns */
         for (k = 1; k <= len; k++)
         {  /* get column index of new non-zero v[i,j] */
            j = ndx[k];
            /* one free location is needed in the j-th column */
            if (vc_len[j] + 1 > vc_cap[j])
            {  /* enlarge the j-th column */
               if (luf_enlarge_col(luf, j, vc_len[j] + 10))
               {  /* overflow of the sparse vector area */
                  ret = 1;
                  goto done;
               }
               /* defragmentation may change row and column pointers of
                  the matrix V */
               p_beg = vr_ptr[p];
               p_end = p_beg + vr_len[p] - 1;
               q_beg = vc_ptr[q];
               q_end = q_beg + vc_len[q] - 1;
            }
            /* add new non-zero v[i,j] to the j-th column */
            j_ptr = vc_ptr[j] + vc_len[j];
            sv_ind[j_ptr] = i;
            vc_len[j]++;
         }
         /* now the i-th row has been completely transformed, therefore
            it can return to the active set with new length */
         rs_prev[i] = 0;
         rs_next[i] = rs_head[vr_len[i]];
         if (rs_next[i] != 0) rs_prev[rs_next[i]] = i;
         rs_head[vr_len[i]] = i;
         /* the largest of absolute values of elements in the i-th row
            is currently unknown */
         vr_max[i] = -1.0;
         /* at least one free location is needed to store the gaussian
            multiplier */
         if (luf->sv_end - luf->sv_beg < 1)
         {  /* there are no free locations at all; defragment SVA */
            luf_defrag_sva(luf);
            if (luf->sv_end - luf->sv_beg < 1)
            {  /* overflow of the sparse vector area */
               ret = 1;
               goto done;
            }
            /* defragmentation may change row and column pointers of the
               matrix V */
            p_beg = vr_ptr[p];
            p_end = p_beg + vr_len[p] - 1;
            q_beg = vc_ptr[q];
            q_end = q_beg + vc_len[q] - 1;
         }
         /* add the element f[i,p], which is the gaussian multiplier,
            to the matrix F */
         luf->sv_end--;
         sv_ind[luf->sv_end] = i;
         sv_val[luf->sv_end] = fip;
         fc_len[p]++;
         /* end of elimination loop */
      }
      /* at this point the q-th (pivot) column should be empty */
      xassert(vc_len[q] == 0);
      /* reset capacity of the q-th column */
      vc_cap[q] = 0;
      /* remove node of the q-th column from the addressing list */
      k = n + q;
      if (sv_prev[k] == 0)
         luf->sv_head = sv_next[k];
      else
         sv_next[sv_prev[k]] = sv_next[k];
      if (sv_next[k] == 0)
         luf->sv_tail = sv_prev[k];
      else
         sv_prev[sv_next[k]] = sv_prev[k];
      /* the p-th column of the matrix F has been completely built; set
         its pointer */
      fc_ptr[p] = luf->sv_end;
      /* walk through the p-th (pivot) row and do the following... */
      for (p_ptr = p_beg; p_ptr <= p_end; p_ptr++)
      {  /* get column index of v[p,j] */
         j = sv_ind[p_ptr];
         /* erase v[p,j] from the working array */
         flag[j] = 0;
         work[j] = 0.0;
         /* the j-th column has been completely transformed, therefore
            it can return to the active set with new length; however
            the special case c_prev[j] = c_next[j] = j means that the
            routine find_pivot excluded the j-th column from the active
            set due to Uwe Suhl's rule, and therefore in this case the
            column can return to the active set only if it is a column
            singleton */
         if (!(vc_len[j] != 1 && cs_prev[j] == j && cs_next[j] == j))
         {  cs_prev[j] = 0;
            cs_next[j] = cs_head[vc_len[j]];
            if (cs_next[j] != 0) cs_prev[cs_next[j]] = j;
            cs_head[vc_len[j]] = j;
         }
      }
done: /* return to the factorizing routine */
      return ret;
}

/***********************************************************************
*  build_v_cols - build the matrix V in column-wise format
*
*  This routine builds the column-wise representation of the matrix V
*  using its row-wise representation.
*
*  If no error occured, the routine returns zero. Otherwise, in case of
*  overflow of the sparse vector area, the routine returns non-zero. */

static int build_v_cols(LUF *luf)
{     int n = luf->n;
      int *vr_ptr = luf->vr_ptr;
      int *vr_len = luf->vr_len;
      int *vc_ptr = luf->vc_ptr;
      int *vc_len = luf->vc_len;
      int *vc_cap = luf->vc_cap;
      int *sv_ind = luf->sv_ind;
      double *sv_val = luf->sv_val;
      int *sv_prev = luf->sv_prev;
      int *sv_next = luf->sv_next;
      int ret = 0;
      int i, i_beg, i_end, i_ptr, j, j_ptr, k, nnz;
      /* it is assumed that on entry all columns of the matrix V are
         empty, i.e. vc_len[j] = vc_cap[j] = 0 for all j = 1, ..., n,
         and have been removed from the addressing list */
      /* count non-zeros in columns of the matrix V; count total number
         of non-zeros in this matrix */
      nnz = 0;
      for (i = 1; i <= n; i++)
      {  /* walk through elements of the i-th row and count non-zeros
            in the corresponding columns */
         i_beg = vr_ptr[i];
         i_end = i_beg + vr_len[i] - 1;
         for (i_ptr = i_beg; i_ptr <= i_end; i_ptr++)
            vc_cap[sv_ind[i_ptr]]++;
         /* count total number of non-zeros */
         nnz += vr_len[i];
      }
      /* store total number of non-zeros */
      luf->nnz_v = nnz;
      /* check for free locations */
      if (luf->sv_end - luf->sv_beg < nnz)
      {  /* overflow of the sparse vector area */
         ret = 1;
         goto done;
      }
      /* allocate columns of the matrix V */
      for (j = 1; j <= n; j++)
      {  /* set pointer to the j-th column */
         vc_ptr[j] = luf->sv_beg;
         /* reserve locations for the j-th column */
         luf->sv_beg += vc_cap[j];
      }
      /* build the matrix V in column-wise format using this matrix in
         row-wise format */
      for (i = 1; i <= n; i++)
      {  /* walk through elements of the i-th row */
         i_beg = vr_ptr[i];
         i_end = i_beg + vr_len[i] - 1;
         for (i_ptr = i_beg; i_ptr <= i_end; i_ptr++)
         {  /* get column index */
            j = sv_ind[i_ptr];
            /* store element in the j-th column */
            j_ptr = vc_ptr[j] + vc_len[j];
            sv_ind[j_ptr] = i;
            sv_val[j_ptr] = sv_val[i_ptr];
            /* increase length of the j-th column */
            vc_len[j]++;
         }
      }
      /* now columns are placed in the sparse vector area behind rows
         in the order n+1, n+2, ..., n+n; so insert column nodes in the
         addressing list using this order */
      for (k = n+1; k <= n+n; k++)
      {  sv_prev[k] = k-1;
         sv_next[k] = k+1;
      }
      sv_prev[n+1] = luf->sv_tail;
      sv_next[luf->sv_tail] = n+1;
      sv_next[n+n] = 0;
      luf->sv_tail = n+n;
done: /* return to the factorizing routine */
      return ret;
}

/***********************************************************************
*  build_f_rows - build the matrix F in row-wise format
*
*  This routine builds the row-wise representation of the matrix F using
*  its column-wise representation.
*
*  If no error occured, the routine returns zero. Otherwise, in case of
*  overflow of the sparse vector area, the routine returns non-zero. */

static int build_f_rows(LUF *luf)
{     int n = luf->n;
      int *fr_ptr = luf->fr_ptr;
      int *fr_len = luf->fr_len;
      int *fc_ptr = luf->fc_ptr;
      int *fc_len = luf->fc_len;
      int *sv_ind = luf->sv_ind;
      double *sv_val = luf->sv_val;
      int ret = 0;
      int i, j, j_beg, j_end, j_ptr, ptr, nnz;
      /* clear rows of the matrix F */
      for (i = 1; i <= n; i++) fr_len[i] = 0;
      /* count non-zeros in rows of the matrix F; count total number of
         non-zeros in this matrix */
      nnz = 0;
      for (j = 1; j <= n; j++)
      {  /* walk through elements of the j-th column and count non-zeros
            in the corresponding rows */
         j_beg = fc_ptr[j];
         j_end = j_beg + fc_len[j] - 1;
         for (j_ptr = j_beg; j_ptr <= j_end; j_ptr++)
            fr_len[sv_ind[j_ptr]]++;
         /* increase total number of non-zeros */
         nnz += fc_len[j];
      }
      /* store total number of non-zeros */
      luf->nnz_f = nnz;
      /* check for free locations */
      if (luf->sv_end - luf->sv_beg < nnz)
      {  /* overflow of the sparse vector area */
         ret = 1;
         goto done;
      }
      /* allocate rows of the matrix F */
      for (i = 1; i <= n; i++)
      {  /* set pointer to the end of the i-th row; later this pointer
            will be set to the beginning of the i-th row */
         fr_ptr[i] = luf->sv_end;
         /* reserve locations for the i-th row */
         luf->sv_end -= fr_len[i];
      }
      /* build the matrix F in row-wise format using this matrix in
         column-wise format */
      for (j = 1; j <= n; j++)
      {  /* walk through elements of the j-th column */
         j_beg = fc_ptr[j];
         j_end = j_beg + fc_len[j] - 1;
         for (j_ptr = j_beg; j_ptr <= j_end; j_ptr++)
         {  /* get row index */
            i = sv_ind[j_ptr];
            /* store element in the i-th row */
            ptr = --fr_ptr[i];
            sv_ind[ptr] = j;
            sv_val[ptr] = sv_val[j_ptr];
         }
      }
done: /* return to the factorizing routine */
      return ret;
}

/***********************************************************************
*  NAME
*
*  luf_factorize - compute LU-factorization
*
*  SYNOPSIS
*
*  #include "glpluf.h"
*  int luf_factorize(LUF *luf, int n, int (*col)(void *info, int j,
*     int ind[], double val[]), void *info);
*
*  DESCRIPTION
*
*  The routine luf_factorize computes LU-factorization of a specified
*  square matrix A.
*
*  The parameter luf specifies LU-factorization program object created
*  by the routine luf_create_it.
*
*  The parameter n specifies the order of A, n > 0.
*
*  The formal routine col specifies the matrix A to be factorized. To
*  obtain j-th column of A the routine luf_factorize calls the routine
*  col with the parameter j (1 <= j <= n). In response the routine col
*  should store row indices and numerical values of non-zero elements
*  of j-th column of A to locations ind[1,...,len] and val[1,...,len],
*  respectively, where len is the number of non-zeros in j-th column
*  returned on exit. Neither zero nor duplicate elements are allowed.
*
*  The parameter info is a transit pointer passed to the routine col.
*
*  RETURNS
*
*  0  LU-factorization has been successfully computed.
*
*  LUF_ESING
*     The specified matrix is singular within the working precision.
*     (On some elimination step the active submatrix is exactly zero,
*     so no pivot can be chosen.)
*
*  LUF_ECOND
*     The specified matrix is ill-conditioned.
*     (On some elimination step too intensive growth of elements of the
*     active submatix has been detected.)
*
*  If matrix A is well scaled, the return code LUF_ECOND may also mean
*  that the threshold pivoting tolerance piv_tol should be increased.
*
*  In case of non-zero return code the factorization becomes invalid.
*  It should not be used in other operations until the cause of failure
*  has been eliminated and the factorization has been recomputed again
*  with the routine luf_factorize.
*
*  REPAIRING SINGULAR MATRIX
*
*  If the routine luf_factorize returns non-zero code, 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 this routine is used for factorizing the basis matrix
*  within the simplex method procedure.
*
*  On exit linearly dependent columns of the (partially transformed)
*  matrix U have numbers rank+1, rank+2, ..., n, where rank is estimated
*  rank of the matrix A stored by the routine to the member luf->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 numbers qq_col[rank+1], qq_col[rank+2], ..., qq_col[n], where
*  qq_col is the column-like representation of the permutation matrix Q.
*  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 pp_row[j], where pp_row is
*  the row-like representation of the permutation matrix P. Thus, each
*  j-th linearly dependent column of the matrix U should be replaced by
*  column of the unity matrix with the number pp_row[j].
*
*  The code that repairs the matrix A may look like follows:
*
*     for (j = rank+1; j <= n; j++)
*     {  replace the column qq_col[j] of the matrix A by the column
*        pp_row[j] of the unity matrix;
*     }
*
*  where rank, pp_row, and qq_col are members of the structure LUF. */

int luf_factorize(LUF *luf, int n, int (*col)(void *info, int j,
      int ind[], double val[]), void *info)
{     int *pp_row, *pp_col, *qq_row, *qq_col;
      double max_gro = luf->max_gro;
      int i, j, k, p, q, t, ret;
      if (n < 1)
         xfault("luf_factorize: n = %d; invalid parameter\n", n);
      if (n > N_MAX)
         xfault("luf_factorize: n = %d; matrix too big\n", n);
      /* invalidate the factorization */
      luf->valid = 0;
      /* reallocate arrays, if necessary */
      reallocate(luf, n);
      pp_row = luf->pp_row;
      pp_col = luf->pp_col;
      qq_row = luf->qq_row;
      qq_col = luf->qq_col;
      /* estimate initial size of the SVA, if not specified */
      if (luf->sv_size == 0 && luf->new_sva == 0)
         luf->new_sva = 5 * (n + 10);
more: /* reallocate the sparse vector area, if required */
      if (luf->new_sva > 0)
      {  if (luf->sv_ind != NULL) xfree(luf->sv_ind);
         if (luf->sv_val != NULL) xfree(luf->sv_val);
         luf->sv_size = luf->new_sva;
         luf->sv_ind = xcalloc(1+luf->sv_size, sizeof(int));
         luf->sv_val = xcalloc(1+luf->sv_size, sizeof(double));
         luf->new_sva = 0;
      }
      /* initialize LU-factorization data structures */
      if (initialize(luf, col, info))
      {  /* overflow of the sparse vector area */
         luf->new_sva = luf->sv_size + luf->sv_size;
         xassert(luf->new_sva > luf->sv_size);
         goto more;
      }
      /* main elimination loop */
      for (k = 1; k <= n; k++)
      {  /* choose a pivot element v[p,q] */
         if (find_pivot(luf, &p, &q))
         {  /* no pivot can be chosen, because the active submatrix is
               exactly zero */
            luf->rank = k - 1;
            ret = LUF_ESING;
            goto done;
         }
         /* 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 = pp_col[p], j = qq_row[q];
         xassert(k <= i && i <= n && k <= j && j <= n);
         /* permute k-th and i-th rows of the matrix U */
         t = pp_row[k];
         pp_row[i] = t, pp_col[t] = i;
         pp_row[k] = p, pp_col[p] = k;
         /* permute k-th and j-th columns of the matrix U */
         t = qq_col[k];
         qq_col[j] = t, qq_row[t] = j;
         qq_col[k] = q, qq_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] */
         if (eliminate(luf, p, q))
         {  /* overflow of the sparse vector area */
            luf->new_sva = luf->sv_size + luf->sv_size;
            xassert(luf->new_sva > luf->sv_size);
            goto more;
         }
         /* check relative growth of elements of the matrix V */
         if (luf->big_v > max_gro * luf->max_a)
         {  /* the growth is too intensive, therefore most probably the
               matrix A is ill-conditioned */
            luf->rank = k - 1;
            ret = LUF_ECOND;
            goto done;
         }
      }
      /* now the matrix U = P*V*Q is upper triangular, the matrix V has
         been built in row-wise format, and the matrix F has been built
         in column-wise format */
      /* defragment the sparse vector area in order to merge all free
         locations in one continuous extent */
      luf_defrag_sva(luf);
      /* build the matrix V in column-wise format */
      if (build_v_cols(luf))
      {  /* overflow of the sparse vector area */
         luf->new_sva = luf->sv_size + luf->sv_size;
         xassert(luf->new_sva > luf->sv_size);
         goto more;
      }
      /* build the matrix F in row-wise format */
      if (build_f_rows(luf))
      {  /* overflow of the sparse vector area */
         luf->new_sva = luf->sv_size + luf->sv_size;
         xassert(luf->new_sva > luf->sv_size);
         goto more;
      }
      /* the LU-factorization has been successfully computed */
      luf->valid = 1;
      luf->rank = n;
      ret = 0;
      /* if there are few free locations in the sparse vector area, try
         increasing its size in the future */
      t = 3 * (n + luf->nnz_v) + 2 * luf->nnz_f;
      if (luf->sv_size < t)
      {  luf->new_sva = luf->sv_size;
         while (luf->new_sva < t)
         {  k = luf->new_sva;
            luf->new_sva = k + k;
            xassert(luf->new_sva > k);
         }
      }
done: /* return to the calling program */
      return ret;
}

/***********************************************************************
*  NAME
*
*  luf_f_solve - solve system F*x = b or F'*x = b
*
*  SYNOPSIS
*
*  #include "glpluf.h"
*  void luf_f_solve(LUF *luf, int tr, double x[]);
*
*  DESCRIPTION
*
*  The routine luf_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 luf, 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 luf_f_solve(LUF *luf, int tr, double x[])
{     int n = luf->n;
      int *fr_ptr = luf->fr_ptr;
      int *fr_len = luf->fr_len;
      int *fc_ptr = luf->fc_ptr;
      int *fc_len = luf->fc_len;
      int *pp_row = luf->pp_row;
      int *sv_ind = luf->sv_ind;
      double *sv_val = luf->sv_val;
      int i, j, k, beg, end, ptr;
      double xk;
      if (!luf->valid)
         xfault("luf_f_solve: LU-factorization is not valid\n");
      if (!tr)
      {  /* solve the system F*x = b */
         for (j = 1; j <= n; j++)
         {  k = pp_row[j];
            xk = x[k];
            if (xk != 0.0)
            {  beg = fc_ptr[k];
               end = beg + fc_len[k] - 1;
               for (ptr = beg; ptr <= end; ptr++)
                  x[sv_ind[ptr]] -= sv_val[ptr] * xk;
            }
         }
      }
      else
      {  /* solve the system F'*x = b */
         for (i = n; i >= 1; i--)
         {  k = pp_row[i];
            xk = x[k];
            if (xk != 0.0)
            {  beg = fr_ptr[k];
               end = beg + fr_len[k] - 1;
               for (ptr = beg; ptr <= end; ptr++)
                  x[sv_ind[ptr]] -= sv_val[ptr] * xk;
            }
         }
      }
      return;
}

/***********************************************************************
*  NAME
*
*  luf_v_solve - solve system V*x = b or V'*x = b
*
*  SYNOPSIS
*
*  #include "glpluf.h"
*  void luf_v_solve(LUF *luf, int tr, double x[]);
*
*  DESCRIPTION
*
*  The routine luf_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 luf, 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 luf_v_solve(LUF *luf, int tr, double x[])
{     int n = luf->n;
      int *vr_ptr = luf->vr_ptr;
      int *vr_len = luf->vr_len;
      double *vr_piv = luf->vr_piv;
      int *vc_ptr = luf->vc_ptr;
      int *vc_len = luf->vc_len;
      int *pp_row = luf->pp_row;
      int *qq_col = luf->qq_col;
      int *sv_ind = luf->sv_ind;
      double *sv_val = luf->sv_val;
      double *b = luf->work;
      int i, j, k, beg, end, ptr;
      double temp;
      if (!luf->valid)
         xfault("luf_v_solve: LU-factorization is not valid\n");
      for (k = 1; k <= n; k++) b[k] = x[k], x[k] = 0.0;
      if (!tr)
      {  /* solve the system V*x = b */
         for (k = n; k >= 1; k--)
         {  i = pp_row[k], j = qq_col[k];
            temp = b[i];
            if (temp != 0.0)
            {  x[j] = (temp /= vr_piv[i]);
               beg = vc_ptr[j];
               end = beg + vc_len[j] - 1;
               for (ptr = beg; ptr <= end; ptr++)
                  b[sv_ind[ptr]] -= sv_val[ptr] * temp;
            }
         }
      }
      else
      {  /* solve the system V'*x = b */
         for (k = 1; k <= n; k++)
         {  i = pp_row[k], j = qq_col[k];
            temp = b[j];
            if (temp != 0.0)
            {  x[i] = (temp /= vr_piv[i]);
               beg = vr_ptr[i];
               end = beg + vr_len[i] - 1;
               for (ptr = beg; ptr <= end; ptr++)
                  b[sv_ind[ptr]] -= sv_val[ptr] * temp;
            }
         }
      }
      return;
}

/***********************************************************************
*  NAME
*
*  luf_a_solve - solve system A*x = b or A'*x = b
*
*  SYNOPSIS
*
*  #include "glpluf.h"
*  void luf_a_solve(LUF *luf, int tr, double x[]);
*
*  DESCRIPTION
*
*  The routine luf_a_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 luf 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 luf_a_solve(LUF *luf, int tr, double x[])
{     if (!luf->valid)
         xfault("luf_a_solve: LU-factorization is not valid\n");
      if (!tr)
      {  /* A = F*V, therefore inv(A) = inv(V)*inv(F) */
         luf_f_solve(luf, 0, x);
         luf_v_solve(luf, 0, x);
      }
      else
      {  /* A' = V'*F', therefore inv(A') = inv(F')*inv(V') */
         luf_v_solve(luf, 1, x);
         luf_f_solve(luf, 1, x);
      }
      return;
}

/***********************************************************************
*  NAME
*
*  luf_delete_it - delete LU-factorization
*
*  SYNOPSIS
*
*  #include "glpluf.h"
*  void luf_delete_it(LUF *luf);
*
*  DESCRIPTION
*
*  The routine luf_delete deletes LU-factorization specified by the
*  parameter luf and frees all the memory allocated to this program
*  object. */

void luf_delete_it(LUF *luf)
{     if (luf->fr_ptr != NULL) xfree(luf->fr_ptr);
      if (luf->fr_len != NULL) xfree(luf->fr_len);
      if (luf->fc_ptr != NULL) xfree(luf->fc_ptr);
      if (luf->fc_len != NULL) xfree(luf->fc_len);
      if (luf->vr_ptr != NULL) xfree(luf->vr_ptr);
      if (luf->vr_len != NULL) xfree(luf->vr_len);
      if (luf->vr_cap != NULL) xfree(luf->vr_cap);
      if (luf->vr_piv != NULL) xfree(luf->vr_piv);
      if (luf->vc_ptr != NULL) xfree(luf->vc_ptr);
      if (luf->vc_len != NULL) xfree(luf->vc_len);
      if (luf->vc_cap != NULL) xfree(luf->vc_cap);
      if (luf->pp_row != NULL) xfree(luf->pp_row);
      if (luf->pp_col != NULL) xfree(luf->pp_col);
      if (luf->qq_row != NULL) xfree(luf->qq_row);
      if (luf->qq_col != NULL) xfree(luf->qq_col);
      if (luf->sv_ind != NULL) xfree(luf->sv_ind);
      if (luf->sv_val != NULL) xfree(luf->sv_val);
      if (luf->sv_prev != NULL) xfree(luf->sv_prev);
      if (luf->sv_next != NULL) xfree(luf->sv_next);
      if (luf->vr_max != NULL) xfree(luf->vr_max);
      if (luf->rs_head != NULL) xfree(luf->rs_head);
      if (luf->rs_prev != NULL) xfree(luf->rs_prev);
      if (luf->rs_next != NULL) xfree(luf->rs_next);
      if (luf->cs_head != NULL) xfree(luf->cs_head);
      if (luf->cs_prev != NULL) xfree(luf->cs_prev);
      if (luf->cs_next != NULL) xfree(luf->cs_next);
      if (luf->flag != NULL) xfree(luf->flag);
      if (luf->work != NULL) xfree(luf->work);
      xfree(luf);
      return;
}

/* eof */