/* 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