ascend4/solver/linsolqr.c

/*
 *  linsol II: Ascend Sparse Linear Solvers
 *  by Benjamin Allan
 *  Created: 2/6/90
 *  Version: $Revision: 1.26 $
 *  Version control file: $RCSfile: linsolqr.c,v $
 *  Date last modified: $Date: 2000/01/25 02:26:58 $
 *  Last modified by: $Author: ballan $
 *
 *  This file is part of the SLV solver.
 *
 *  Copyright (C) 1994 Benjamin Andrew Allan
 *
 *  Based on linsol:
 *  Copyright (C) 1990 Karl Michael Westerberg
 *  Copyright (C) 1993 Joseph Zaher
 *  Copyright (C) 1994 Joseph Zaher, Benjamin Andrew Allan
 *  and on sqr.pas v1.5:
 *  Copyright (C) 1994 Boyd Safrit
 *
 *  The SLV solver 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 2 of the
 *  License, or (at your option) any later version.
 *
 *  The SLV solver is distributed in 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 the program; if not, write to the Free Software Foundation,
 *  Inc., 675 Mass Ave, Cambridge, MA 02139 USA.  Check the file named
 *  COPYING.  COPYING is found in ../compiler.
 */

/*
 *  Implementation Notes:
 *   9/95: there are many places in linsolqr.c where lower level matrix
 *   operations have been replaced by higher level ones. These used to
 *   be flagged as such so that people who wanted to construct newer
 *   operations from the lowest level mtx operators could see how to
 *   use the lowlevel operators. Since it has been decided that the old
 *   linsol will NOT be discarded completely (we like having easy
 *   benchmarks to look better than), the user wishing to code new
 *   operations is directed back to linsol to see how lowlevel operations
 *   can be expensively done.
 *
 *   10/95: (BAA) I attempted to optimize number_drag through use of
 *   pointer arithmetic. Testing the resulting code without optimization
 *   turned on gave about 20% improvement over the original on an alpha
 *   using the native compiler. However, with optimization on the array
 *   subscripted version which is left here won by ~40%.
 *   The results for gcc were similar: pointer 10% better unoptimized,
 *   subscripted 30% better when -O2.
 *   For both compilers -O3 was worse! gcc was ~1.8x slower than the
 *   native cc both optimized and not.
 *   I have not tested a loop unrolled drag, but I expect unrolling is
 *   exactly what the optimizer is doing, so why bother?
 */

#include <math.h>
#include <stdarg.h>
#include "utilities/ascConfig.h"
#include "utilities/ascMalloc.h"
#include "utilities/ascPanic.h"
#include "utilities/mem.h"
#include "utilities/set.h"
#include "general/tm_time.h"
#include "solver/mtx.h"
#include "solver/linsolqr.h"

#ifdef __WIN32__
#define copysign _copysign
#endif /* __WIN32__ */

#define LINSOLQR_DEBUG FALSE
/* qr code debugging flag */
#define LINSOL_DEBUG FALSE
/* linsol debugging flag */

#define D_ZERO (real64)0.0
#define ZERO (int32)0
#define D_ONE (real64)1.0

/* housekeeping flags */
#define BUILD_DEAD_CODE 0
#undef BUILD_KIRK_CODE


#define RBADEBUG 0
/* turns on spew that will generate output files
 * comparable to and in the same location as for rankiba2
 * user specific hardcoded pathnames are involved.
 * don't try this at home.
 */
#if RBADEBUG
extern FILE *gscr;
#endif

/*
 * timing messages disable flag.
 */
int g_linsolqr_timing = 1;

/***************************************************************************\
  Data structures for linsol
  Any time you modify these, please verify linsolqr_size is still correct.
\***************************************************************************/
struct rhs_list {
   real64 *rhs;            /* Vector of rhs values */
   real64 *varvalue;       /* Solution of the linear system */
   boolean solved;               /* ? has the rhs been solved */
   boolean transpose;            /* ? transpose the coef matrix */
   struct rhs_list *next;
};

struct qr_fill_t {
  int32 cnt;     /* number of real nonzeros in a column */
  int32 col;     /* column number prior to being sorted */
  int32 fill;    /* fill if this (or if col?) were the pivot column  */
  int32 overlap; /* scratch for fillable calculation. */
};
/* fill counting structure for sparse qr */

struct lu_auxdata {       /* all indexed by cur row or cur col number */
  real64 *pivlist;  /* vector of pivots (diagonal of L for ranki) */
  real64 *tmp;      /* row elimination, dependency buffer */
  int32 cap;        /* current capacity of the vectors if not null */
};
/* structure for lu algorithms */

struct qr_auxdata {       /* all indexed by cur row or cur col number */
  mtx_matrix_t hhvects;   /* slave matrix of sys->factors for holding u's */
  mtx_sparse_t sp;        /* sparse data space */
  real64 *alpha;    /* column norms for the Q\R factor matrix */
  real64 *sigma;    /* column norms for the S (R inverse) matrix */
  real64 *tau;      /* column tau of the Householder transforms */
  real64 *hhcol;    /* Householder transform scratch vector */
  real64 *hhrow;    /* Householder transform scratch vector */
  real64 nu;        /* current condition number of inverse */
  real64 anu;       /* condition number of coef */
  real64 asubnu;    /* condition number of active region in coef */
  struct qr_fill_t *fill; /* fill counting structure */
  mtx_region_t facreg;    /* region over which qr is done if nonsquare */
  int32 cap;        /* current capacity of these vectors if not null */
};
/* structure for CondQR algorithms */

struct linsolqr_header {
   int integrity;
   enum factor_class fclass;     /* Type of factoring expected */
   enum factor_method fmethod;   /* Method factoring expected */
   enum reorder_method rmethod;  /* Type of most recent reordering */
   int32 capacity;         /* Capacity of arrays */
   mtx_matrix_t coef;            /* Coefficient matrix */
   struct rhs_list *rl;          /* List of rhs vectors */
   int rlength;                  /* Length of rhs list */
   real64 pivot_zero;      /* Smallest acceptable pivot */
   real64 ptol;            /* Pivot selection tolerance */
   real64 ctol;            /* Condition selection tolerance */
   real64 dtol;            /* Matrix entry drop tolerance */
   mtx_matrix_t factors;         /* Matrix with UL and dependence info (ranki),
                                    or L and dependence info  (ranki2),
                                    or R and dependence info (qr methods)  */
   mtx_matrix_t inverse;         /* Matrix with U multipliers (ranki2) or R^-1
                                    (condqr), in both cases slave of factors */
   boolean factored;             /* ? has the matrix been factored */
   boolean rowdeps;              /* ? have row dependencies been calc'ed */
   boolean coldeps;              /* ? have col dependencies been calc'ed */
   mtx_range_t rng;              /* Pivot range */
   mtx_region_t reg;             /* Bounding region given */
   int32 rank;             /* Rank of the matrix */
   real64 smallest_pivot;  /* Smallest pivot accepted */
   struct qr_auxdata *qrdata;    /* Data vectors for qr methods */
   struct lu_auxdata *ludata;    /* Data vectors for lu methods */
};
/* linsol main structure */

/***************************************************************************\
  string functions for linsol io
\***************************************************************************/

char *linsolqr_rmethods() {
  static char names[]="Natural, SPK1, TSPK1";
  return names;
}

char *linsolqr_fmethods() {
  static char names[] =
    "SPK1/RANKI,SPK1/RANKI+ROW,Fast-SPK1/RANKI,Fast-SPK1/RANKI+ROW,Fastest-SPK1/MR-RANKI,CondQR,CPQR";
  return names;
}

/** **/
enum reorder_method linsolqr_rmethod_to_enum(char *name) {
  if (strcmp(name,"SPK1")==0) return spk1;
  if (strcmp(name,"TSPK1")==0) return tspk1;
  if (strcmp(name,"Natural")==0) return natural;
  return unknown_r;
}

enum factor_method linsolqr_fmethod_to_enum(char *name) {
  if (strcmp(name,"KIRK-STUFF")==0) return ranki_ka;
  if (strcmp(name,"SPK1/RANKI")==0) return ranki_kw;
  if (strcmp(name,"SPK1/RANKI+ROW")==0) return ranki_jz;
  if (strcmp(name,"Fast-SPK1/RANKI")==0) return ranki_kw2;
  if (strcmp(name,"Fast-SPK1/RANKI+ROW")==0) return ranki_jz2;
  if (strcmp(name,"Fastest-SPK1/MR-RANKI")==0) return ranki_ba2;
  if (strcmp(name,"CondQR")==0) return cond_qr;
  if (strcmp(name,"CPQR")==0) return plain_qr;
  return unknown_f;
}

enum factor_class linsolqr_fmethod_to_fclass(enum factor_method fm)
{
  switch (fm) {
    /* ranki-like things */
    case ranki_kw:
    case ranki_jz:
    case ranki_ba2:
    case ranki_kw2:
    case ranki_jz2:
    case ranki_ka:
      return ranki;
    /* implemented qr things */
    case plain_qr:
      return s_qr;
    /* other stuff unimplemented. all fall through */
    case cond_qr:
    case opt_qr:
    case ls_qr:
    case gauss_ba2:
    case symmetric_lu:
    case unknown_f:
    default:
      return unknown_c;
  }
}

/** **/
char *linsolqr_enum_to_rmethod(enum reorder_method meth) {
  switch (meth) {
  case spk1: return "SPK1";
  case tspk1: return "TSPK1";
  case natural: return "Natural";
  default: return "Unknown reordering method.";
  }
}

char *linsolqr_enum_to_fmethod(enum factor_method meth) {
  switch (meth) {
  case ranki_kw: return "SPK1/RANKI";
  case ranki_jz: return "SPK1/RANKI+ROW";
  case ranki_ba2: return "Fastest-SPK1/MR-RANKI";
  case ranki_kw2: return "Fast-SPK1/RANKI";
  case ranki_jz2: return "Fast-SPK1/RANKI+ROW";
  case ranki_ka: return "KIRK-STUFF";
  case cond_qr: return "CondQR";
  case plain_qr: return "CPQR";
  default: return "Unknown factorization method.";
  }
}

/** **/
char *linsolqr_rmethod_description(enum reorder_method meth) {
  static char sspk1[]="SPK1 reordering ala Stadtherr.";
  static char stspk1[]="SPK1 reordering column wise ala Stadtherr.";
  static char nat[]="Ordering as received from user";
  switch (meth) {
  case spk1: return sspk1;
  case tspk1: return stspk1;
  case natural: return nat;
  default: return "Unknown reordering method.";
  }
}

char *linsolqr_fmethod_description(enum factor_method meth) {
  static char ex_ranki[] =
    "SPK1 reordering with RANKI LU factoring ala Stadtherr.";
  static char ex_rankib[] =
    "SPK1 reordering with RANKI LU factoring ala Stadtherr, but fast.";
  static char ex_rankij[]="SPK1/RANKI LU with pseudo-complete pivoting";
  static char ex_rankik[]="KIRK-STUFF/RANKI LU with pseudo-complete pivoting";
  static char ex_condqr[]="Sparse QR with condition controlled pivoting.";
  static char ex_plainqr[]="Sparse QR with column pivoting.";
  switch (meth) {
  case ranki_kw: return ex_ranki;
  case ranki_kw2: return ex_ranki;
  case ranki_ba2: return ex_rankib;
  case ranki_jz: return ex_rankij;
  case ranki_jz2: return ex_rankij;
  case ranki_ka: return ex_rankik;
  case cond_qr: return ex_condqr;
  case plain_qr: return ex_plainqr;
  default: return "Unknown factorization method.";
  }
}

/***************************************************************************\
  internal check routines
\***************************************************************************/
#define OK        ((int)439828676)
#define DESTROYED ((int)839276847)
#define CHECK_SYSTEM(a) check_system((a),__FILE__,__LINE__)
static int check_system(linsolqr_system_t sys, char *file, int line)
/**
 ***  Checks the system handle. Return 0 if ok, 1 otherwise.
 **/
{
   if( ISNULL(sys) ) {
      FPRINTF(stderr,"ERROR:  (linsolqr) check_system\n");
      FPRINTF(stderr,"        NULL system handle.\n");
      return 1;
   }

   switch( sys->integrity ) {
   case OK:
      return 0;
   case DESTROYED:
      FPRINTF(stderr,"ERROR:  (linsolqr) check_system, line %d\n",line);
      FPRINTF(stderr,"        System was recently destroyed.\n");
      break;
   default:
      FPRINTF(stderr,"ERROR:  (linsolqr) check_system, line %d\n",line);
      FPRINTF(stderr,"        System was reused or never created.\n");
      break;
   }
   return 1;
}

/***************************************************************************\
  external calls (and some of their internals)
\***************************************************************************/
static void destroy_rhs_list(struct rhs_list *rl)
/**
 ***  Destroys rhs list.
 **/
{
   while( NOTNULL(rl) ) {
      struct rhs_list *p;
      p = rl;
      rl = rl->next;
      if( NOTNULL(p->varvalue)  )
     ascfree( (POINTER)(p->varvalue) );
      ascfree( (POINTER)p );
   }
}

static struct rhs_list *find_rhs(struct rhs_list *rl,real64 *rhs)
/**
 ***  Searches for rhs in rhs list, returning it or NULL if not found.
 **/
{
   for( ; NOTNULL(rl) ; rl = rl->next )
      if( rl->rhs == rhs )
         break;
   return(rl);
}

static struct lu_auxdata *create_ludata()
/**
 *** Creates an empty zeroed ludata struct.
 *** Filled up by insure_lu_capacity.
 **/
{
  struct lu_auxdata *d;
  d=(struct lu_auxdata *)asccalloc(1,sizeof(struct lu_auxdata));
  return d;
}

static struct qr_auxdata *create_qrdata()
/**
 *** Creates an empty zeroed qrdata struct.
 *** Filled up by insure_qr_capacity.
 **/
{
  struct qr_auxdata *d;
  d=(struct qr_auxdata *)asccalloc(1,sizeof(struct qr_auxdata));
  return d;
}

static void destroy_ludata(struct lu_auxdata *d) {
/**
 *** Conditionally destroys a ludata struct and its parts.
 **/
  if (NOTNULL(d)) {
    if (NOTNULL(d->pivlist))  ascfree(d->pivlist);
    d->pivlist=NULL;
    if (NOTNULL(d->tmp))  ascfree(d->tmp);
    d->tmp=NULL;
    ascfree(d);
  }
}

static void destroy_qrdata(struct qr_auxdata *d) {
/**
 *** Conditionally destroys a qrdata struct and its parts.
 **/
  if (NOTNULL(d)) {
    if (NOTNULL(d->alpha)) ascfree(d->alpha);
    if (NOTNULL(d->sigma)) ascfree(d->sigma);
    if (NOTNULL(d->tau))   ascfree(d->tau);
    if (NOTNULL(d->hhcol)) ascfree(d->hhcol);
    if (NOTNULL(d->hhrow)) ascfree(d->hhrow);
    if (NOTNULL(d->fill))  ascfree(d->fill);

    ascfree(d);
  }
}

linsolqr_system_t linsolqr_create()
{
   linsolqr_system_t sys;

   sys = (linsolqr_system_t)ascmalloc( sizeof(struct linsolqr_header) );
   sys->fclass = unknown_c;
   sys->fmethod = unknown_f;
   sys->rmethod = unknown_r;
   sys->integrity = OK;
   sys->capacity = 0;
   sys->coef = NULL;
   sys->rl = NULL;
   sys->rlength = 0;
   sys->pivot_zero = 1e-12;           /* default value */
   sys->ptol = 0.1;                   /* default value */
   sys->ctol = 0.1;                   /* default value */
   sys->dtol = LINQR_DROP_TOLERANCE;  /* default value */
   sys->factors = NULL;
   sys->inverse = NULL;
   sys->factored = FALSE;
   sys->rowdeps = FALSE;
   sys->coldeps = FALSE;
   sys->rng.low = sys->rng.high = -1;
   sys->reg.row.low = sys->reg.row.high = -1;
   sys->reg.col.low = sys->reg.col.high = -1;
   sys->rank = -1;
   sys->smallest_pivot = MAXDOUBLE;
   sys->qrdata = NULL;
   sys->ludata = NULL;
   return(sys);
}

void linsolqr_destroy(linsolqr_system_t sys)
{
   if(CHECK_SYSTEM(sys)) {
     FPRINTF(stderr,"Bad linsolqr_system_t found. Not destroyed.\n");
     return;
   }
   if( NOTNULL(sys->coef) ) {
     FPRINTF(stderr,"Notice:  (linsolqr) non-NULL coefficient matrix\n");
     FPRINTF(stderr,"                    in linsolqr system not being\n");
     FPRINTF(stderr,"                    destroyed with linsolqr_system_t.\n");
   }
   if( NOTNULL(sys->inverse) )
      mtx_destroy(sys->inverse);
   if( NOTNULL(sys->factors) )
      mtx_destroy(sys->factors);
   destroy_rhs_list(sys->rl);
   destroy_qrdata(sys->qrdata);
   destroy_ludata(sys->ludata);
   sys->integrity = DESTROYED;
   ascfree( (POINTER)sys );
}

void linsolqr_set_matrix(linsolqr_system_t sys,mtx_matrix_t mtx)
{
   if(CHECK_SYSTEM(sys)) {
     FPRINTF(stderr,"Bad linsolqr_system_t found. coef mtx not set.\n");
     return;
   }
   sys->coef = mtx;
}


void linsolqr_set_region(linsolqr_system_t sys,mtx_region_t region)
{
   if(CHECK_SYSTEM(sys)) {
     FPRINTF(stderr,"Bad linsolqr_system_t found. coef mtx not set.\n");
     return;
   }
   sys->reg = region;
}


mtx_matrix_t linsolqr_get_matrix(linsolqr_system_t sys)
{
   CHECK_SYSTEM(sys);
   return( sys->coef );
}

mtx_region_t *linsolqr_get_region(linsolqr_system_t sys)
{
    return &(sys->reg);
}

void linsolqr_add_rhs(linsolqr_system_t sys,
                      real64 *rhs,
                      boolean transpose)
{
   struct rhs_list *rl;

   if(CHECK_SYSTEM(sys)) {
     FPRINTF(stderr,"Bad linsolqr_system_t found. rhs not added.\n");
     return;
   }
   rl = find_rhs(sys->rl,rhs);
   if( NOTNULL(rl) ) {
      return;
   } else {
     if( NOTNULL(rhs) ) {
       rl = (struct rhs_list *)ascmalloc( sizeof(struct rhs_list) );
       rl->rhs = rhs;
       rl->varvalue = NULL;
       rl->solved = FALSE;
       rl->transpose = transpose;
       rl->next = sys->rl;
       sys->rl = rl;
       ++(sys->rlength);
       if (sys->capacity>0)  {
         rl->varvalue =
          (real64 *)ascmalloc(sys->capacity*sizeof(real64));
       }
     }
   }
}

void linsolqr_remove_rhs(linsolqr_system_t sys,real64 *rhs)
{
   struct rhs_list **q;

   if(CHECK_SYSTEM(sys)) {
     FPRINTF(stderr,"Bad linsolqr_system_t found. rhs not removed.\n");
     return;
   }
   for( q = &(sys->rl) ; NOTNULL(*q) ; q = &((*q)->next) )
      if( (*q)->rhs == rhs )
     break;
   if( NOTNULL(*q) ) {
      struct rhs_list *p;
      p = *q;
      *q = p->next;
      if( NOTNULL(p->varvalue) )
     ascfree( (POINTER)(p->varvalue) );
      ascfree( (POINTER)p );
      --(sys->rlength);
   } else if( NOTNULL(rhs) ) {
      FPRINTF(stderr,"ERROR:  (linsolqr) linsolqr_remove_rhs\n");
      FPRINTF(stderr,"        Rhs does not exist.\n");
   }
}

int32 linsolqr_number_of_rhs(linsolqr_system_t sys)
{
   CHECK_SYSTEM(sys);
   return( sys->rlength );
}

real64 *linsolqr_get_rhs(linsolqr_system_t sys,int n)
{
   struct rhs_list *rl;
   int count;

   CHECK_SYSTEM(sys);

   count = sys->rlength - 1 - n;
   if( count < 0 )  return(NULL);
   for( rl = sys->rl ; count > 0 && NOTNULL(rl) ; rl = rl->next )
      --count;
   return( ISNULL(rl) ? NULL : rl->rhs );
}

void linsolqr_matrix_was_changed(linsolqr_system_t sys)
{
   if(CHECK_SYSTEM(sys)) {
     FPRINTF(stderr,
       "Bad linsolqr_system_t found. matrix change message ignored.\n");
     return;
   }
   sys->rowdeps = sys->coldeps = sys->factored = FALSE;
}

void linsolqr_rhs_was_changed(linsolqr_system_t sys, real64 *rhs)
{
   struct rhs_list *rl;

   if(CHECK_SYSTEM(sys)) {
     FPRINTF(stderr,"Bad linsolqr_system_t found. rhs change ignored.\n");
     return;
   }
   rl = find_rhs(sys->rl,rhs);
   if( NOTNULL(rl) ) {
      rl->solved = FALSE;
   } else if( NOTNULL(rhs) ) {
      FPRINTF(stderr,"ERROR:  (linsolqr) linsolqr_rhs_was_modified\n");
      FPRINTF(stderr,"        Rhs does not exist.\n");
   }
}

void linsolqr_set_pivot_zero(linsolqr_system_t sys,real64 pivot_zero)
{
   if(CHECK_SYSTEM(sys)) {
     FPRINTF(stderr,"Bad linsolqr_system_t found. set_pivot_zero ignored.\n");
     return;
   }
   if( pivot_zero < D_ZERO ) {
      FPRINTF(stderr,"ERROR:  (linsolqr) linsolqr_set_pivot_zero\n");
      FPRINTF(stderr,"        Invalid pivot zero of %g\n",pivot_zero);
   } else {
      sys->pivot_zero = pivot_zero;
      linsolqr_matrix_was_changed(sys);
   }
}

real64 linsolqr_pivot_zero(linsolqr_system_t sys)
{
   CHECK_SYSTEM(sys);
   return( sys->pivot_zero );
}

void linsolqr_set_pivot_tolerance(linsolqr_system_t sys, real64 ptol)
{
   if(CHECK_SYSTEM(sys)) {
     FPRINTF(stderr,"Bad linsolqr_system_t found. set_pivot_tol ignored.\n");
     return;
   }
   if( ptol < D_ZERO || ptol > D_ONE ) {
      FPRINTF(stderr,"ERROR:  (linsolqr) linsolqr_set_pivot_tolerance\n");
      FPRINTF(stderr,"        Invalid pivot tolerance of %g\n",ptol);
   } else {
      sys->ptol = ptol;
      linsolqr_matrix_was_changed(sys);
   }
}

void linsolqr_set_condition_tolerance(linsolqr_system_t sys, real64 ctol)
{
   if(CHECK_SYSTEM(sys)) {
     FPRINTF(stderr,
       "Bad linsolqr_system_t found. set_condition_tolerance ignored.\n");
     return;
   }
   if( ctol < D_ZERO || ctol > D_ONE ) {
      FPRINTF(stderr,"ERROR:  (linsolqr) linsolqr_set_condition_tolerance\n");
      FPRINTF(stderr,"        Invalid condition tolerance of %g\n",ctol);
   } else {
      sys->ctol = ctol;
      linsolqr_matrix_was_changed(sys);
   }
}

void linsolqr_set_drop_tolerance(linsolqr_system_t sys, real64 dtol)
{
   if(CHECK_SYSTEM(sys)) {
     FPRINTF(stderr,
       "Bad linsolqr_system_t found. set_drop_tolerance ignored.\n");
     return;
   }
   if( dtol < D_ZERO || dtol > D_ONE ) {
      FPRINTF(stderr,"ERROR:  (linsolqr) linsolqr_set_drop_tolerance\n");
      FPRINTF(stderr,"        Invalid drop tolerance of %g\n",dtol);
   } else {
      sys->dtol = dtol;
      linsolqr_matrix_was_changed(sys);
   }
}

real64 linsolqr_pivot_tolerance(linsolqr_system_t sys)
{
   CHECK_SYSTEM(sys);
   return( sys->ptol );
}

real64 linsolqr_condition_tolerance(linsolqr_system_t sys)
{
   CHECK_SYSTEM(sys);
   return( sys->ctol );
}

real64 linsolqr_drop_tolerance(linsolqr_system_t sys)
{
   CHECK_SYSTEM(sys);
   return( sys->dtol );
}

extern enum factor_class linsolqr_fclass(linsolqr_system_t sys)
{
   CHECK_SYSTEM(sys);
   return( sys->fclass );
}
extern enum factor_method linsolqr_fmethod(linsolqr_system_t sys)
{
   CHECK_SYSTEM(sys);
   return( sys->fmethod );
}
extern enum reorder_method linsolqr_rmethod(linsolqr_system_t sys)
{
   CHECK_SYSTEM(sys);
   return( sys->rmethod );
}

int32 linsolqr_rank(linsolqr_system_t sys)
{
   CHECK_SYSTEM(sys);
   if( !sys->factored ) {
      FPRINTF(stderr,"ERROR:  (linsolqr) linsolqr_rank\n");
      FPRINTF(stderr,"        System not factored yet.\n");
   }
   return(sys->rank);
}

real64 linsolqr_smallest_pivot(linsolqr_system_t sys)
{
   CHECK_SYSTEM(sys);
#if LINSOL_DEBUG
   if( !sys->factored ) {
      FPRINTF(stderr,"ERROR:  (linsolqr) linsolqr_smallest_pivot\n");
      FPRINTF(stderr,"        System not factored yet.\n");
   }
#endif
   return(sys->smallest_pivot);
}

/***************************************************************************\
  Commonly used internal functions for sparse linear solvers based on mtx.
   void insure_capacity(sys)
   void insure_lu_capacity(sys)
   void insure_qr_capacity(sys)
   void forward_substitute(sys,rvec,transpose)
   void backward_substitute(sys,rvec,transpose)
   macro SQR(x)
   int find_pivot_number(vec,len,tol,eps,ivec,rvec,maxi)
\***************************************************************************/

static real64 *raise_capacity(real64 *vec,
                                    int32 oldcap,
                                    int32 newcap)
/**
 ***  Raises the capacity of the array and returns a new array.
 ***  It is assumed that oldcap < newcap.  vec is destroyed or
 ***  returned as appropriate.
 ***  If !NDEBUG, the vector expanded is also set to 0.
 **/
{
  real64 *newvec=NULL;
  int i;
  if (newcap < oldcap) {
#ifndef NDEBUG
    for (i = 0; i < newcap; i++) {
      newvec[i] = 0.0;
    }
#endif
    return vec;
  }
  if (NOTNULL(vec)) {
    /* don't call realloc on null with newcap 0 or it frees */
    newvec=(real64 *)ascrealloc(vec,(newcap * sizeof(real64)));
  } else {
    newvec=(newcap > 0 ?
      (real64 *)ascmalloc( newcap * sizeof(real64) ) : NULL);
  }
#ifndef NDEBUG
  for (i = 0; i < newcap; i++) {
    newvec[i] = 0.0;
  }
#endif
  return newvec;
}

static int32 *raise_qr_int_capacity(int32 *vec,
                                          int32 oldcap,
                                          int32 newcap)
/**
 ***  Raises the capacity of the index array and returns a new array.
 ***  It is assumed that oldcap < newcap.  vec is destroyed or
 ***  returned as appropriate. vec returned is zeroed.
 ***  calling with newcap=0 does not force deallocation.
 **/
{
  int32 *newvec=NULL;
  if (newcap < oldcap)
    return vec;
  if (NOTNULL(vec)) {/* don't call realloc on null with newcap 0 or it frees */
    newvec=(int32 *)ascrealloc(vec,(newcap * sizeof(int32)));
    mtx_zero_int32(vec,newcap);
  } else {
    newvec= (newcap > 0 ?
      (int32 *)asccalloc( newcap , sizeof(int32) ) : NULL);
  }
  return newvec;
}

static real64 *raise_qr_real_capacity(real64 *vec,
                                       int32 oldcap,
                                       int32 newcap)
/**
 ***  Raises the capacity of the real array and returns a new array.
 ***  It is assumed that oldcap < newcap.  vec is destroyed or
 ***  returned as appropriate. vec returned is zeroed.
 ***  calling with newcap=0 does not force deallocation.
 **/
{
  real64 *newvec=NULL;
  if (newcap < oldcap)
    return vec;
  if (NOTNULL(vec)) {/* don't call realloc on null with newcap 0 or it frees */
    newvec=(real64 *)ascrealloc(vec,(newcap * sizeof(real64)));
    mtx_zero_real64(vec,newcap);
  } else {
    newvec= (newcap > 0 ?
      (real64 *)asccalloc( newcap , sizeof(real64) ) : NULL);
  }
  return newvec;
}

static struct qr_fill_t *raise_qr_fill_capacity(struct qr_fill_t *vec,
                                                int32 oldcap,
                                                int32 newcap)
/**
 ***  Raises the capacity of the fill array and returns a new array.
 ***  It is assumed that oldcap < newcap.  vec is destroyed or
 ***  returned as appropriate. vec returned is zeroed.
 ***  calling with newcap=0 does not force deallocation.
 **/
{
  struct qr_fill_t *newvec=NULL;
  if (newcap < oldcap)
    return vec;
  if (NOTNULL(vec)) {/* don't call realloc on null with newcap 0 or it frees */
    newvec=(struct qr_fill_t *)
      ascrealloc(vec,(newcap * sizeof(struct qr_fill_t)));
    mtx_zero_char((char *)vec,newcap*sizeof(struct qr_fill_t));
  } else {
    newvec= (newcap > 0 ?
      (struct qr_fill_t *)asccalloc( newcap , sizeof(struct qr_fill_t) ) :
      NULL);
  }
  return newvec;
}

static void insure_capacity(linsolqr_system_t sys)
/**
 ***  Insures that the capacity of all of the solution vectors
 ***  for each rhs is large enough.
 ***  The above implies a malloc if varvalue is null.
 ***  Assumes varvalue are at sys->capacity already.
 **/
{
   int32 req_cap;
   req_cap = mtx_capacity(sys->coef);

   if( req_cap > sys->capacity ) {
      struct rhs_list *rl;

      for( rl = sys->rl ; NOTNULL(rl) ; rl = rl->next )
         rl->varvalue = raise_capacity(rl->varvalue,sys->capacity,req_cap);
      sys->capacity = req_cap;
   }
}

static void insure_lu_capacity(linsolqr_system_t sys)
/**
 ***  Insures that the capacity of all of the ludata vectors.
 ***  The above implies a malloc if vector is null.
 ***  If not null, implies an extension if needed.
 **/
{
  int32 req_cap;
  if (!sys || !(sys->coef) || !(sys->ludata)) {
    FPRINTF(stderr,"linsolqr: insure_lu_capacity called with NULL pointer");
    return;
  }
  req_cap = mtx_capacity(sys->coef);

  if( req_cap > sys->ludata->cap ) {
    sys->ludata->pivlist =
      raise_capacity(sys->ludata->pivlist,sys->ludata->cap,req_cap);
    sys->ludata->tmp =
      raise_capacity(sys->ludata->tmp,sys->ludata->cap,req_cap);
    sys->ludata->cap = req_cap;
  }
}

static void insure_qr_capacity(linsolqr_system_t sys)
/**
 ***  Insures that the capacity of all of the qrdata vectors
 ***  is large enough.
 ***  The above implies a calloc if vector is null.
 ***  If not null, implies an extension and zeroing of the vector.
 ***  Also zeroes the simple elements of the qrdata.
 **/
{
  int32 req_cap;
  if (!sys || !(sys->coef) || !(sys->qrdata)) {
    FPRINTF(stderr,"linsolqr: insure_qr_capacity called with NULL pointer");
    return;
  }
  req_cap = mtx_capacity(sys->coef);

  if( req_cap > sys->qrdata->cap ) {
    sys->qrdata->alpha =
      raise_qr_real_capacity(sys->qrdata->alpha,sys->qrdata->cap,req_cap);
    sys->qrdata->sigma =
      raise_qr_real_capacity(sys->qrdata->sigma,sys->qrdata->cap,req_cap);
    sys->qrdata->tau =
      raise_qr_real_capacity(sys->qrdata->tau,sys->qrdata->cap,req_cap);
    sys->qrdata->hhcol =
      raise_qr_real_capacity(sys->qrdata->hhcol,sys->qrdata->cap,req_cap);
    sys->qrdata->hhrow =
      raise_qr_real_capacity(sys->qrdata->hhrow,sys->qrdata->cap,req_cap);

    sys->qrdata->sp.data =
      raise_qr_real_capacity(sys->qrdata->sp.data,sys->qrdata->cap,req_cap);
    sys->qrdata->sp.idata =
      raise_qr_int_capacity(sys->qrdata->sp.idata,sys->qrdata->cap,req_cap);

    sys->qrdata->fill =
      raise_qr_fill_capacity(sys->qrdata->fill,sys->qrdata->cap,req_cap);

    sys->qrdata->cap = req_cap;
    sys->qrdata->sp.cap = req_cap;
    sys->qrdata->sp.len = 0;
  }
  sys->qrdata->nu=D_ZERO;
  sys->qrdata->anu=D_ZERO;
  sys->qrdata->asubnu=D_ZERO;
}

static void forward_substitute(linsolqr_system_t sys,
                               real64 *arr,
                               boolean transpose)
/**
 ***  Forward substitute.  It is assumed that the L (or U) part of
 ***  sys->factors is computed.  This function converts c to x in place.  The
 ***  values are stored in arr indexed by original row number (or original
 ***  column number).
 ***
 ***     transpose = FALSE:                  transpose = TRUE:
 ***                                          T
 ***     L x = c                             U x = c
 ***
 ***  The following formula holds:
 ***     0<=k<r ==> x(k) = [c(k) - L(k,(0..k-1)) dot x(0..k-1)] / L(k,k)
 ***  or
 ***     0<=k<r ==> x(k) = [c(k) - U((0..k-1),k) dot x(0..k-1)] / U(k,k)
 ***
 ***  U(k,k) is 1 by construction. L(k,k) is stored in sys->ludata->pivlist[k]
 ***  and in matrix.
 **/
{
   mtx_range_t dot_rng;
   mtx_coord_t nz;
   real64 sum, *pivlist;
   mtx_matrix_t mtx;
   int32 dotlim;
   boolean nonzero_found=FALSE;

   mtx=sys->factors;
   pivlist=sys->ludata->pivlist;
   dot_rng.low = sys->rng.low;
   dotlim=dot_rng.low+sys->rank;
   if (transpose) {     /* arr is indexed by original column number */
     for( nz.col=dot_rng.low; nz.col < dotlim; ++(nz.col) ) {
       register int32 org_col;

       dot_rng.high = nz.col - 1;
       org_col = mtx_col_to_org(mtx,nz.col);
       if (arr[org_col]!=D_ZERO) nonzero_found=TRUE;
       if (nonzero_found) {
         sum=mtx_col_dot_full_org_vec(mtx,nz.col,arr,&dot_rng,TRUE);
         /* arr[org_col] = (arr[org_col] - sum)  / D_ONE */;
         arr[org_col] -=sum;
       }
     }

   } else {             /* arr is indexed by original row number */
     for( nz.row=dot_rng.low; nz.row < dotlim; ++(nz.row) ) {
       register int32 org_row;

       dot_rng.high = nz.row - 1;
       org_row = mtx_row_to_org(mtx,nz.row);
       if (arr[org_row]!=D_ZERO) nonzero_found=TRUE;
       if (nonzero_found) {
         sum = mtx_row_dot_full_org_vec(mtx,nz.row,arr,&dot_rng,TRUE);
/*
         nz.col = nz.row;
         arr[org_row] = (arr[org_row] - sum) / mtx_value(mtx,&nz);
*/
         arr[org_row] = (arr[org_row] - sum) / pivlist[nz.row];
       }
     }
   }
}

static void backward_substitute(linsolqr_system_t sys,
                                real64 *arr,
                                boolean transpose)
/**
 ***  Backward substitute.  It is assumed that the U (or L) part of
 ***  sys->factors is computed.  This function converts rhs to c in place.  The
 ***  values are stored in arr indexed by original row number (or original
 ***  column number).
 ***
 ***    transpose = FALSE:                  transpose = TRUE:
 ***                                         T
 ***    U c = rhs                           L c = rhs
 ***
 ***  The following formula holds:
 ***  transpose=FALSE: (the usual for J.dx=-f where rhs is -f)
 ***     0<=k<r ==> c(k) = [rhs(k) - U(k,(k+1..r-1)) dot c(k+1..r-1)] / U(k,k)
 ***     working up from the bottom.
 ***  or
 ***  transpose=TRUE:
 ***     0<=k<r ==> c(k) = [rhs(k) - L((k+1..r-1),k) dot c(k+1..r-1)] / L(k,k)
 ***     working right to left.
 ***
 ***  U(k,k) is 1 by construction. L(k,k) is stored in sys->ludata->pivlist[k]
 ***  and in matrix.
 **/
{
   mtx_range_t dot_rng;
   mtx_coord_t nz;
   real64 sum, *pivlist;
   mtx_matrix_t mtx;
   int32 dotlim;
   boolean nonzero_found=FALSE; /* once TRUE, substitution must be done
                                   over remaining rows/cols */

   dot_rng.high = sys->rng.low + sys->rank - 1;
   dotlim=sys->rng.low;
   mtx=sys->factors;
   pivlist=sys->ludata->pivlist;
   if (transpose) {     /* arr is indexed by original column number */
     for( nz.col = dot_rng.high ; nz.col >= dotlim ; --(nz.col) ) {
       register int32 org_col;

       dot_rng.low = nz.col + 1;

       org_col = mtx_col_to_org(mtx,nz.col);
       if (arr[org_col]!=D_ZERO) nonzero_found=TRUE;
       if (nonzero_found) {
         sum= mtx_col_dot_full_org_vec(mtx,nz.col,arr,&dot_rng,TRUE);
/*
  reminders:
         nz.row = nz.col;
         arr[org_col] = (arr[org_col] - sum) / mtx_value(mtx,&nz);
*/
         arr[org_col] = (arr[org_col] - sum) / pivlist[nz.col];
       }
     }
   } else {             /* arr is indexed by original row number */
     /* we are working from the bottom up */
     for( nz.row = dot_rng.high ; nz.row >= dotlim ; --(nz.row) ) {
       register int32 org_row;

       dot_rng.low = nz.row + 1;
       org_row = mtx_row_to_org(mtx,nz.row);

       if (arr[org_row]!=D_ZERO) nonzero_found=TRUE;
       if (nonzero_found) {
         sum= mtx_row_dot_full_org_vec(mtx,nz.row,arr,&dot_rng,TRUE);
/*
  reminders:
         nz.row = nz.col;
         arr[org_row] = (arr[org_row] - sum) / D_ONE;
*/
         arr[org_row] -= sum;
       }
     }
   }
}

#ifdef SQR
#undef SQR
#endif
#define SQR(x) ((x)*(x))
/*
 *   int32 find_pivot_number(vec,len,tol,eps,ivec,rvec,maxi);
 *   real64 *vec,*rvec,tol;
 *   int32 len, *ivec, *maxi;
 *
 *   Search array vec of positive numbers for the first entry, k, which passes
 *   (vec[k] >= tol*vecmax) where vecmax is the largest value in vec.
 *   vec is an array len long. rvec, ivec are (worst case) len long.
 *   *maxi will be set to the index of the first occurence of the largest
 *   number in vec which is >= eps.
 *   0.0 <= tol <= 1.
 *   If tol <=0, no search is done (*maxi = 0, k = 0).
 *   Eps is an absolute number which vec values must be >= to before
 *   being eligible for the tol test.
 *
 *   Best case, highly nonlinear feature:
 *     if on entry *maxi = len, we will do a simplified search for the
 *     first k which satisfies tol and eps tests, based on the value
 *     stored at vec[len-1].
 *
 *   We could allocate and deallocate rvec/ivec since they are
 *   temporaries, but it is more efficient for caller to manage that.
 *   If tol>=1.0, rvec and ivec are ignored and may be NULL.
 *   NOTE: we are going on value vec[i], not fabs(vec[i]).
 *   Returns k.
 * Remember: GIGO.
 * Failure modes- if all numbers are <= 0.0, k and *maxi will be returned 0.
 * Failure modes- if all numbers are < eps, k and *maxi will be returned 0.
 *
 * find_pivot_index, should we implement, would search an int vector.
 *
 * Hint: if you know you have not increased any values in vec and have not
 * changed the value at vec[maxi] before a subsequent search, then call
 * with len = previous maxi+1 to reduce the search time.
 */
static int32 find_pivot_number(const real64 *vec,
                                     const int32 len,
                                     const real64 tol,
                                     const real64 eps,
                                     int32 * const ivec,
                                     real64 * const rvec,
                                     int32 * const maxi)
{
  int32 j,k;
  if (tol <= D_ZERO ) {
    *maxi = 0;
    return 0;
  }
  if (tol >= D_ONE) {
    register real64 biggest = MAX(eps,D_ZERO);
    /* get the biggest, period */
    k = 0;
    for (j=0; j < len; j++) {
      if (vec[j] > biggest) {
        biggest = vec[j];
        k = j;
      }
    }
    *maxi = k;
  } else {
    int32 bigi;
    bigi = 0;
    rvec[0] = eps;
    ivec[0] = 0;
    if (*maxi==len) {
      /* cheap search against eps and tol, no list. */
      register real64 thold;
      thold = tol * vec[len-1];
      if (thold >= eps) {
        for (k = 0; k < len && vec[k] < thold; k++);
      } else {
        for (k = 0; k < len && vec[k] < eps; k++);
      }
      /* adjust maxi to still point at last location, as indicated by len */
      (*maxi)--;
    } else {
      real64 rlast;
      rlast = eps;
      /* create short list */
      for (j=0; j < len; j++) {
        while (j < len && vec[j] <= rlast) j++; /* skip to next max */
        if (j < len) {
          ivec[bigi] = j;
          rvec[bigi] = vec[j];
          bigi++;  /* bigi ends up being len of rvec or, if vec all<= eps, 0 */
          rlast = rvec[bigi-1];
        }
      }
      if (bigi == 0) {
        *maxi = k = 0;
      } else {
        register real64 thold;
        /* Note: if bigi is enormous, we should do a binary search,
           not linear. */
        *maxi = ivec[bigi-1];
        thold = tol * rvec[bigi-1];
        /* get pivot from shortlist, searching backward */
        if (thold >= eps) {
          for (k = bigi-1; k >= 0 && rvec[k] >= thold; k--);
        } else {
          for (k = bigi-1; k >= 0 && rvec[k] >= eps; k--);
        }
        /* translate pivot to vec index */
        k = ivec[k+1];
      }
    }
  }
  return k;
}

static void calc_dependent_rows_ranki1(linsolqr_system_t sys)
{
  mtx_coord_t nz;
  real64 value;
  mtx_range_t colrange;
  mtx_range_t rowrange;
  real64 *lc;
  mtx_matrix_t mtx;

  sys->rowdeps = TRUE;
  if( ( (sys->reg.row.low == sys->rng.low) &&
        ( sys->reg.row.high == sys->rng.low+sys->rank-1 )
      ) || sys->rank==0 )
    return;

  lc = sys->ludata->tmp;
  colrange.low = sys->rng.low;
  colrange.high = colrange.low + sys->rank - 1;
  rowrange.low = sys->rng.high;
  rowrange.high = sys->rng.low+sys->rank;
  mtx=sys->factors;

  nz.row = sys->reg.row.low;
  for( ; nz.row <= sys->reg.row.high; nz.row++ ) {
    if( nz.row == sys->rng.low ) {
      nz.row = rowrange.high-1;
      continue;
    }
    mtx_zero_real64(lc,(sys->capacity));
    mtx_org_row_vec(mtx,nz.row,lc,&colrange);
    if( nz.row < rowrange.high || nz.row > rowrange.low )
      backward_substitute(sys,lc,TRUE);
    forward_substitute(sys,lc,TRUE);
    mtx_clear_row(mtx,nz.row,&colrange);
    for( nz.col=colrange.low; nz.col <= colrange.high; nz.col++ ) {
      value = lc[mtx_col_to_org(mtx,nz.col)];
      if( value != D_ZERO ) mtx_fill_value(mtx,&nz,value);
    }
  }
}


static void calc_dependent_cols_ranki1(linsolqr_system_t sys)
{
  mtx_coord_t nz;
  real64 value;
  mtx_range_t rowrange;
  mtx_range_t colrange;
  real64 *lc;
  mtx_matrix_t mtx;

  sys->coldeps = TRUE;
  if( ( (sys->reg.col.low == sys->rng.low) &&
        ( sys->reg.col.high == sys->rng.low+sys->rank-1 )
      ) || sys->rank==0 )
    return;

  lc = sys->ludata->tmp;
  rowrange.low = sys->rng.low;
  rowrange.high = rowrange.low + sys->rank - 1;
  colrange.high = sys->rng.low+sys->rank;
  colrange.low = sys->rng.high;
  mtx=sys->factors;

  nz.col = sys->reg.col.low;
  for( ; nz.col <= sys->reg.col.high; nz.col++ ) {
    if( nz.col == sys->rng.low ) {
      nz.col = colrange.high-1;
      continue;
    }
    mtx_zero_real64(lc,sys->capacity);
    mtx_org_col_vec(mtx,nz.col,lc,&rowrange);
    if( nz.col < colrange.high || nz.col > colrange.low )
      backward_substitute(sys,lc,FALSE);
    forward_substitute(sys,lc,FALSE);
    mtx_clear_col(mtx,nz.col,&rowrange);
    for( nz.row=rowrange.low; nz.row <= rowrange.high; nz.row++ ) {
        value = lc[mtx_row_to_org(mtx,nz.row)];
      if( value != D_ZERO ) mtx_fill_value(mtx,&nz,value);
    }
  }
}

/**
 *** Given a matrix and a diagonal range, sets the diagonal elements to 0
 *** but does not delete them in the range low to high inclusive.
 *** worst cost= k*(nonzeros in rows of range treated).
 *** likely cost= k*(nonzeros in column pivoted rows of range treated)
 ***              since you can smartly put in the diagonal last and find
 ***              it first on many occasions.
 *** This is a good bit cheaper than deleting the elements unless one
 *** expects thousands of solves.
 **/
static void zero_diagonal_elements(mtx_matrix_t mtx,
                                   int32 low, int32 high)
{
  mtx_coord_t nz;
  for (nz.row = low; nz.row <= high; nz.row++) {
    nz.col = nz.row;
    mtx_set_value(mtx,&nz,0.0);
  }
}

static void zero_unpivoted_vars(linsolqr_system_t sys,
                real64 *varvalues,
                boolean transpose)
/**
 ***  Sets the values of unpivoted variables to zero.
 **/
{
   int32 ndx,order;
   mtx_matrix_t mtx;

   mtx = sys->factors;
   order = mtx_order(mtx);
   for( ndx = 0 ; ndx < sys->rng.low ; ++ndx )
     if (transpose)
       varvalues[mtx_col_to_org(mtx,ndx)] = D_ZERO;
     else
       varvalues[mtx_row_to_org(mtx,ndx)] = D_ZERO;

   for( ndx = sys->rng.low + sys->rank ; ndx < order ; ++ndx )
     if (transpose)
       varvalues[mtx_col_to_org(mtx,ndx)] = D_ZERO;
     else
       varvalues[mtx_row_to_org(mtx,ndx)] = D_ZERO;
}

#if LINSOL_DEBUG
static void debug_out_factors(FILE *fp,linsolqr_system_t sys)
/**
 ***  Outputs permutation and values of the nonzero elements in the
 ***  factor matrix square region given by sys->rng.
 **/
{
   mtx_region_t reg;
   reg.row.low=reg.col.low=sys->rng.low;
   reg.row.high=reg.col.high=sys->rng.high;
   mtx_write_region(fp,sys->factors,&reg);
}
#endif /* LINSOL_DEBUG */

/***************************************************************************\
  Reordering functions for SPK1, and possibly for other schemes to be
  implemented later.
  The stuff here is almost, but not quite, black magic. Don't tinker with it.
\***************************************************************************/

/*********************************
 begin of spk1 stuff
*********************************/
struct reorder_list {            /* List of rows/columns and their counts. */
   int32 ndx;
   int32 count;
   struct reorder_list *next;
};

struct reorder_vars {
   mtx_matrix_t mtx;
   mtx_region_t reg;             /* Current active region */
   int32 colhigh;          /* Original value of reg.col.high */
   struct reorder_list *tlist;   /* Array of (enough) list elements */
   int32 *rowcount;        /* rowcount[reg.row.low .. reg.row.high] */
};

static void adjust_row_count(struct reorder_vars *vars,int32 removed_col)
/**
 ***  Adjusts the row counts to account for the (to be) removed column.
 **/
{
   mtx_coord_t nz;
   real64 value;
   nz.col = removed_col;
   nz.row = mtx_FIRST;
   while( value = mtx_next_in_col(vars->mtx,&nz,&(vars->reg.row)),
     nz.row != mtx_LAST )
      --(vars->rowcount[nz.row]);
}

static void assign_row_and_col(struct reorder_vars *vars,
                          int32 row,
                          int32 col)
/**
 ***  Assigns the given row to the given column, moving the row and column
 ***  to the beginning of the active region and removing them (readjusting
 ***  the region).  The row counts are NOT adjusted.  If col == mtx_NONE,
 ***  then no column is assigned and the row is moved to the end of the
 ***  active block instead.  Otherwise, it is assumed that the column
 ***  is active.
 **/
{
   if( col == mtx_NONE ) {
      mtx_swap_rows(vars->mtx,row,vars->reg.row.high);
      vars->rowcount[row] = vars->rowcount[vars->reg.row.high];
      --(vars->reg.row.high);
   } else {
      mtx_swap_rows(vars->mtx,row,vars->reg.row.low);
      vars->rowcount[row] = vars->rowcount[vars->reg.row.low];
      mtx_swap_cols(vars->mtx,col,vars->reg.col.low);
      ++(vars->reg.row.low);
      ++(vars->reg.col.low);
   }
}

static void push_column_on_stack(struct reorder_vars *vars,int32 col)
/**
 ***  Pushes the given column onto the stack.  It is assumed that the
 ***  column is active.  Row counts are adjusted.
 **/
{
   adjust_row_count(vars,col);
   mtx_swap_cols(vars->mtx,col,vars->reg.col.high);
   --(vars->reg.col.high);
}

static int32 pop_column_from_stack(struct reorder_vars *vars)
/**
 ***  Pops the column on the "top" of the stack off of the stack and
 ***  returns the column index, where it now lies in the active region.
 ***  If the stack is empty, mtx_NONE is returned.  Row counts are NOT
 ***  adjusted (this together with a subsequent assignment of this column
 ***  ==> no row count adjustment necessary).
 **/
{
   if( vars->reg.col.high < vars->colhigh )
      return(++(vars->reg.col.high));
   else
      return( mtx_NONE );
}

static void assign_null_rows(struct reorder_vars *vars)
/**
 ***  Assigns empty rows, moving them to the assigned region.  It is
 ***  assumed that row counts are correct.  Columns are assigned off the
 ***  stack.
 **/
{
   int32 row;

   for( row = vars->reg.row.low ; row <= vars->reg.row.high ; ++row )
      if( vars->rowcount[row] == 0 )
         assign_row_and_col(vars , row , pop_column_from_stack(vars));
}

static void forward_triangularize(struct reorder_vars *vars)
/**
 ***  Forward triangularizes the region, assigning singleton rows with their
 ***  one and only incident column until there are no more.  The row counts
 ***  must be correct, and they are updated.
 **/
{
   boolean change;
   mtx_coord_t nz;
   real64 value;

   do {
      change = FALSE;
      for( nz.row = vars->reg.row.low ;
          nz.row <= vars->reg.row.high && vars->rowcount[nz.row] != 1;
          ++nz.row ) ;
      if( nz.row <= vars->reg.row.high ) {
         /* found singleton row */
         nz.col = mtx_FIRST; /* this is somehow coming back with nz.col -1 */
         value = mtx_next_in_row(vars->mtx,&nz,&(vars->reg.col));
         adjust_row_count(vars,nz.col);
         assign_row_and_col(vars,nz.row,nz.col);
         change = TRUE;
      }
   } while( change );
}

static int32 select_row(struct reorder_vars *vars)
/**
 ***  Selects a row and returns its index.  It is assumed that there is a
 ***  row.  Row counts must be correct.  vars->tlist will be used.
 **/
{
   int32 min_row_count;
   int32 max_col_count;
   int32 row;
   int32 i, nties=-2;  /* # elements currently defined in vars->tlist */
   int32 sum;
   mtx_coord_t nz;
   real64 value;
   mtx_matrix_t mtx;
   mtx_range_t *colrng, *rowrng;

   /* Set to something > any possible value */
   min_row_count = vars->reg.col.high-vars->reg.col.low+2;
   for( row = vars->reg.row.low ; row <= vars->reg.row.high ; ++row )
      if( vars->rowcount[row] <= min_row_count ) {
         if( vars->rowcount[row] < min_row_count ) {
            min_row_count = vars->rowcount[row];
            nties = 0;
         }
         vars->tlist[nties++].ndx = row;
      }
   /**
    ***  vars->tlist[0..nties-1] is a list of row numbers which tie for
    ***   minimum row count.
    **/

   max_col_count = -1;   /* < any possible value */
   i = 0;
   mtx = vars->mtx;
   colrng=&(vars->reg.col);
   rowrng=&(vars->reg.row);
   while( i < nties ) {

      sum = 0;
      nz.row = vars->tlist[i].ndx;
      nz.col = mtx_FIRST;
      while( value = mtx_next_in_row(mtx,&nz,colrng),
        nz.col != mtx_LAST )
         sum += mtx_nonzeros_in_col(mtx,nz.col,rowrng);
      if( sum > max_col_count ) {
         max_col_count = sum;
         row = nz.row;
      }
      i++;
   }
   /**
    ***  Now row contains the row with the minimum row count, which has the
    ***  greatest total column count of incident columns among all rows with
    ***  the same (minimum) row count.  Select it.
    **/
   return(row);
}

static void spk1_reorder(struct reorder_vars *vars)
/**
 ***  Reorders the assigned matrix vars->mtx within the specified bounding
 ***  block region vars->reg.  The region is split into 6 subregions during
 ***  reordering:  the rows are divided in two, and the columns divided in
 ***  three.  Initially everything is in the active subregion.  Ultimately,
 ***  everything will be assigned.
 ***
 ***          <-- assigned -->|<-- active-->|<-- on stack -->|
 ***     ----+----------------+-------------+----------------+
 ***       a |                |             |                |
 ***       s |                |             |                |
 ***       s |                |             |                |
 ***       i |    (SQUARE)    |             |                |
 ***       g |                |             |                |
 ***       n |                |             |                |
 ***       e |                |             |                |
 ***       d |                |             |                |
 ***     ----+----------------+-------------+----------------+
 ***       a |                |             |                |
 ***       c |                |    ACTIVE   |                |
 ***       t |                |    REGION   |                |
 ***       i |                |             |                |
 ***       v |                |             |                |
 ***       e |                |             |                |
 ***     ----+----------------+-------------+----------------+
 ***
 ***  The algorithm is roughly as follows:
 ***    (1) select a row (by some criterion).
 ***    (2) push columns incident on that row onto the stack in decreasing
 ***        order of their length.
 ***    (3) pop first column off the stack and assign it to the selected
 ***        row.
 ***    (4) forward-triangularize (assign singleton rows with their one
 ***        and only incident column, until no longer possible).
 ***
 ***    (1)-(4) should be repeated until the active subregion becomes empty.
 ***
 ***  Everything above was written as though the entire matrix is
 ***  involved.  In reality, only the relevant square region is involved.
 **/
{
   int32 row, size;
   int32 *rowcount_array_origin;
   mtx_matrix_t mtx;

   size = MAX(vars->reg.row.high,vars->reg.col.high) + 1;
   vars->tlist = size > 0 ? (struct reorder_list *)
      ascmalloc( size*sizeof(struct reorder_list) ) : NULL;
   vars->rowcount = rowcount_array_origin = size > 0 ? (int32 *)
      ascmalloc( size*sizeof(int32) ) : NULL;
   mtx = vars->mtx;

   vars->colhigh = vars->reg.col.high;
   /* Establish row counts */
   for( row = vars->reg.row.low ; row <= vars->reg.row.high ; ++row )
      vars->rowcount[row] =
     mtx_nonzeros_in_row(mtx,row,&(vars->reg.col));

   while(TRUE) {
      struct reorder_list *head;
      int32 nelts;   /* # elements "allocated" from vars->tlist */
      mtx_coord_t nz;
      real64 value;

      forward_triangularize(vars);
      assign_null_rows(vars);
      if( vars->reg.row.low>vars->reg.row.high ||
     vars->reg.col.low>vars->reg.col.high ) {
     /* Active region is now empty, done */
     if( NOTNULL(vars->tlist) )
        ascfree( vars->tlist );
     if( NOTNULL(rowcount_array_origin) )
        ascfree( rowcount_array_origin );
         return;
      }

      head = NULL;
      nelts = 0;
      nz.row = select_row(vars);
      nz.col = mtx_FIRST;
      while( value = mtx_next_in_row(mtx,&nz,&(vars->reg.col)),
        nz.col != mtx_LAST ) {
         struct reorder_list **q,*p;

         p = &(vars->tlist[nelts++]);
         p->ndx = mtx_col_to_org(mtx,nz.col);
         p->count = mtx_nonzeros_in_col(mtx,nz.col,&(vars->reg.row));
         for( q = &head; *q && (*q)->count > p->count; q = &((*q)->next) );
         p->next = *q;
         *q = p;
      }
      /**
       ***  We now have a list of columns which intersect the selected row.
       ***  The list is in order of decreasing column count.
       **/

      /* Push incident columns on stack */
      for( ; NOTNULL(head) ; head = head->next )
         push_column_on_stack(vars,mtx_org_to_col(mtx,head->ndx));

      /* Pop column off stack and assign to selected row */
      assign_row_and_col(vars , nz.row , pop_column_from_stack(vars));
   }
   /*  Not reached. */
}

/*********************************
 end of spk1 stuff
*********************************/
/*********************************
 begin of tspk1 stuff
*********************************/
struct creorder_list {            /* List of columns/rows and their counts. */
   int32 ndx;
   int32 count;
   struct creorder_list *next;
};

struct creorder_vars {
   mtx_matrix_t mtx;
   mtx_region_t reg;             /* Current active region */
   int32 rowhigh;          /* Original value of reg.row.high */
   struct creorder_list *tlist;  /* Array of (enough) list elements */
   int32 *colcount;        /* colcount[reg.col.low .. reg.col.high] */
};

static void adjust_col_count(struct creorder_vars *vars,int32 removed_row)
/**
 ***  Adjusts the column counts to account for the (to be) removed row.
 **/
{
   mtx_coord_t nz;
   real64 value;
   nz.row = removed_row;
   nz.col = mtx_FIRST;
   while( value = mtx_next_in_row(vars->mtx,&nz,&(vars->reg.col)),
     nz.col != mtx_LAST )
      --(vars->colcount[nz.col]);
}

static void assign_col_and_row(struct creorder_vars *vars,
                          int32 col,
                          int32 row)
/**
 ***  Assigns the given row to the given column, moving the row and column
 ***  to the beginning of the active region and removing them (readjusting
 ***  the region).  The col counts are NOT adjusted.  If col == mtx_NONE,
 ***  then no column is assigned and the col is moved to the end of the
 ***  active block instead.  Otherwise, it is assumed that the row
 ***  is active.
 **/
{
   if( row == mtx_NONE ) {
      mtx_swap_cols(vars->mtx,col,vars->reg.col.high);
      vars->colcount[col] = vars->colcount[vars->reg.col.high];
      --(vars->reg.col.high);
   } else {
      mtx_swap_cols(vars->mtx,col,vars->reg.col.low);
      vars->colcount[col] = vars->colcount[vars->reg.col.low];
      mtx_swap_rows(vars->mtx,row,vars->reg.row.low);
      ++(vars->reg.col.low);
      ++(vars->reg.row.low);
   }
}

static void push_row_on_stack(struct creorder_vars *vars,int32 row)
/**
 ***  Pushes the given row onto the stack.  It is assumed that the
 ***  row is active.  Col counts are adjusted.
 **/
{
   adjust_col_count(vars,row);
   mtx_swap_rows(vars->mtx,row,vars->reg.row.high);
   --(vars->reg.row.high);
}

static int32 pop_row_from_stack(struct creorder_vars *vars)
/**
 ***  Pops the row on the "top" of the stack off of the stack and
 ***  returns the row index, where it now lies in the active region.
 ***  If the stack is empty, mtx_NONE is returned.  Col counts are NOT
 ***  adjusted (this together with a subsequent assignment of this row
 ***  ==> no col count adjustment necessary).
 **/
{
   if( vars->reg.row.high < vars->rowhigh )
      return(++(vars->reg.row.high));
   else
      return( mtx_NONE );
}

static void assign_null_cols(struct creorder_vars *vars)
/**
 ***  Assigns empty cols, moving them to the assigned region.  It is
 ***  assumed that col counts are correct.  Rows are assigned off the
 ***  stack.
 **/
{
   int32 col;

   for( col = vars->reg.col.low ; col <= vars->reg.col.high ; ++col )
      if( vars->colcount[col] == 0 )
         assign_col_and_row(vars , col , pop_row_from_stack(vars));
}

static void cforward_triangularize(struct creorder_vars *vars)
/**
 ***  Forward triangularizes the region, assigning singleton columns with their
 ***  one and only incident row until there are no more.  The column counts
 ***  must be correct, and they are updated.
 **/
{
   boolean change;

   do {
      mtx_coord_t nz;
      real64 value;
      change = FALSE;
      for( nz.col = vars->reg.col.low ;
          nz.col <= vars->reg.col.high && vars->colcount[nz.col] != 1;
          ++nz.col ) ;
      if( nz.col <= vars->reg.col.high ) {
         /* found singleton col */
         nz.row = mtx_FIRST;
         value = mtx_next_in_col(vars->mtx,&nz,&(vars->reg.row));
         adjust_col_count(vars,nz.row);
         assign_col_and_row(vars,nz.col,nz.row);
         change = TRUE;
      }
   } while( change );
}

static int32 select_col(struct creorder_vars *vars)
/**
 ***  Selects a col and returns its index.  It is assumed that there is a
 ***  col.  Col counts must be correct.  vars->tlist will be used.
 **/
{
   int32 min_col_count;
   int32 max_row_count;
   int32 col;
   int32 i, nties=-2;  /* # elements currently defined in vars->tlist */
   int32 sum;
   mtx_coord_t nz;
   real64 value;
   mtx_matrix_t mtx;
   mtx_range_t *colrng, *rowrng;

   /* Set to something > any possible value */
   min_col_count = vars->reg.row.high-vars->reg.row.low+2;
   for( col = vars->reg.col.low ; col <= vars->reg.col.high ; ++col )
      if( vars->colcount[col] <= min_col_count ) {
         if( vars->colcount[col] < min_col_count ) {
            min_col_count = vars->colcount[col];
            nties = 0;
         }
         vars->tlist[nties++].ndx = col;
      }
   /**
    ***  vars->tlist[0..nties-1] is a list of row numbers which tie for
    ***   minimum col count.
    **/

   max_row_count = -1;   /* < any possible value */
   i = 0;
   mtx = vars->mtx;
   rowrng=&(vars->reg.row);
   colrng=&(vars->reg.col);
   while( i < nties ) {

      sum = 0;
      nz.row = mtx_FIRST;
      nz.col = vars->tlist[i].ndx;
      while( value = mtx_next_in_col(mtx,&nz,rowrng),
        nz.row != mtx_LAST )
         sum += mtx_nonzeros_in_row(mtx,nz.row,colrng);
      if( sum > max_row_count ) {
         max_row_count = sum;
         col = nz.col;
      }
      i++;
   }
   /**
    ***  Now col contains the col with the minimum col count, which has the
    ***  greatest total row count of incident rows among all cols with
    ***  the same (minimum) col count.  Select it.
    **/
   return(col);
}

static void tspk1_reorder(struct creorder_vars *vars)
/**
 ***  Transpose the picture and explanation that follows:
 ***  Reorders the assigned matrix vars->mtx within the specified bounding
 ***  block region vars->reg.  The region is split into 6 subregions during
 ***  reordering:  the rows are divided in two, and the columns divided in
 ***  three.  Initially everything is in the active subregion.  Ultimately,
 ***  everything will be assigned.
 ***
 ***          <-- assigned -->|<-- active-->|<-- on stack -->|
 ***     ----+----------------+-------------+----------------+
 ***       a |                |             |                |
 ***       s |                |             |                |
 ***       s |                |             |                |
 ***       i |    (SQUARE)    |             |                |
 ***       g |                |             |                |
 ***       n |                |             |                |
 ***       e |                |             |                |
 ***       d |                |             |                |
 ***     ----+----------------+-------------+----------------+
 ***       a |                |             |                |
 ***       c |                |    ACTIVE   |                |
 ***       t |                |    REGION   |                |
 ***       i |                |             |                |
 ***       v |                |             |                |
 ***       e |                |             |                |
 ***     ----+----------------+-------------+----------------+
 ***
 ***  The algorithm is roughly as follows:
 ***    (1) select a row (by some criterion).
 ***    (2) push columns incident on that row onto the stack in decreasing
 ***        order of their length.
 ***    (3) pop first column off the stack and assign it to the selected
 ***        row.
 ***    (4) forward-triangularize (assign singleton rows with their one
 ***        and only incident column, until no longer possible).
 ***
 ***    (1)-(4) should be repeated until the active subregion becomes empty.
 ***
 ***  Everything above was written as though the entire matrix is
 ***  involved.  In reality, only the relevant square region is involved.
 **/
{
   int32 col, size;
   int32 *colcount_array_origin;
   mtx_matrix_t mtx;

   size = vars->reg.col.high - vars->reg.col.low + 1;
   size = MAX(size,vars->reg.row.high - vars->reg.row.low + 1);
   vars->tlist = size > 0 ? (struct creorder_list *)
      ascmalloc( size*sizeof(struct creorder_list) ) : NULL;
   colcount_array_origin = size > 0 ? (int32 *)
      ascmalloc( size*sizeof(int32) ) : NULL;
   vars->colcount =
     NOTNULL(colcount_array_origin) ?
       colcount_array_origin - vars->reg.col.low : NULL;
   mtx = vars->mtx;

   vars->rowhigh = vars->reg.row.high;
   /* Establish col counts */
   for( col = vars->reg.col.low ; col <= vars->reg.col.high ; ++col )
      vars->colcount[col] =
     mtx_nonzeros_in_col(mtx,col,&(vars->reg.row));

   while(TRUE) {
      struct creorder_list *head;
      int32 nelts;   /* # elements "allocated" from vars->tlist */
      mtx_coord_t nz;
      real64 value;

      cforward_triangularize(vars);
      assign_null_cols(vars);
      if( vars->reg.col.low > vars->reg.col.high ||
     vars->reg.row.low > vars->reg.row.high ) {
     /* Active region is now empty, done */
     if( NOTNULL(vars->tlist) )
        ascfree( vars->tlist );
     if( NOTNULL(colcount_array_origin) )
        ascfree( colcount_array_origin );
         return;
      }

      head = NULL;
      nelts = 0;
      nz.row = mtx_FIRST;
      nz.col = select_col(vars);
      while( value = mtx_next_in_col(mtx,&nz,&(vars->reg.row)),
        nz.row != mtx_LAST ) {
         struct creorder_list **q,*p;

         p = &(vars->tlist[nelts++]);
         p->ndx = mtx_row_to_org(mtx,nz.row);
         p->count = mtx_nonzeros_in_row(mtx,nz.row,&(vars->reg.col));
         for( q = &head; *q && (*q)->count > p->count; q = &((*q)->next) );
         p->next = *q;
         *q = p;
      }
      /**
       ***  We now have a list of columns which intersect the selected row.
       ***  The list is in order of decreasing column count.
       **/

      /* Push incident rows on stack */
      for( ; NOTNULL(head) ; head = head->next )
         push_row_on_stack(vars,mtx_org_to_row(mtx,head->ndx));

      /* Pop row off stack and assign to selected col */
      assign_col_and_row(vars , nz.col , pop_row_from_stack(vars));
   }
   /*  Not reached. */
}

/*********************************
 end of tspk1 stuff
*********************************/

static boolean nonempty_row(mtx_matrix_t mtx,int32 row)
/**
 ***  ? row not empty in mtx.
 **/
{
   mtx_coord_t nz;
   real64 value;
   value = mtx_next_in_row(mtx, mtx_coord(&nz,row,mtx_FIRST), mtx_ALL_COLS);
   return( nz.col != mtx_LAST );
}

static boolean nonempty_col(mtx_matrix_t mtx,int32 col)
/**
 ***  ? column not empty in mtx.
 **/
{
   mtx_coord_t nz;
   real64 value;
   value = mtx_next_in_col(mtx, mtx_coord(&nz,mtx_FIRST,col), mtx_ALL_ROWS);
   return( nz.row != mtx_LAST );
}

static void determine_pivot_range(linsolqr_system_t sys)
/**
 ***  Calculates sys->rng from sys->coef.
 **/
{
   sys->reg.row.low = sys->reg.row.high = -1;
   sys->reg.col.low = sys->reg.col.high = -1;

   for( sys->rng.high=mtx_order(sys->coef) ; --(sys->rng.high) >= 0 ; ) {
      if( nonempty_row(sys->coef,sys->rng.high) &&
     sys->reg.row.high < 0 )
     sys->reg.row.high = sys->rng.high;
      if( nonempty_col(sys->coef,sys->rng.high) &&
     sys->reg.col.high < 0 )
     sys->reg.col.high = sys->rng.high;
      if( nonempty_row(sys->coef,sys->rng.high) &&
         nonempty_col(sys->coef,sys->rng.high) )
         break;
   }

   for( sys->rng.low=0 ; sys->rng.low <= sys->rng.high ; ++(sys->rng.low) ) {
      if( nonempty_row(sys->coef,sys->rng.low) &&
     sys->reg.row.low < 0 )
     sys->reg.row.low = sys->rng.low;
      if( nonempty_col(sys->coef,sys->rng.low) &&
     sys->reg.col.low < 0 )
     sys->reg.col.low = sys->rng.low;
      if( nonempty_row(sys->coef,sys->rng.low) &&
         nonempty_col(sys->coef,sys->rng.low) )
         break;
   }
}

static void square_region(linsolqr_system_t sys,mtx_region_t *region)
/**
 *** Get the largest square confined to the diagonal within the region given
 *** and set sys->rng accordingly.
 **/
{
      sys->reg = *region;
      sys->rng.low = MAX(region->row.low,region->col.low);
      sys->rng.high = MIN(region->row.high,region->col.high);
}

static int ranki_reorder(linsolqr_system_t sys,mtx_region_t *region)
/**
 ***  The region to reorder is first isolated by truncating the region
 ***  provided to the largest square region confined to the matrix diagonal.
 ***  It is presumed it will contain no empty rows or columns and will
 ***  provide the basis of candidate pivots when factoring.
 **/
{
   struct reorder_vars vars;
   CHECK_SYSTEM(sys);
   if (sys->fclass != ranki) {
     FPRINTF(stderr,"linsolqr_reorder:  spk1 reorder called on system\n");
     FPRINTF(stderr,"                   with inappropriate factor class\n");
     return 1;
   }
   if( region == mtx_ENTIRE_MATRIX ) determine_pivot_range(sys);
   else square_region(sys,region);

   vars.mtx = sys->coef;
   vars.reg.row.low = vars.reg.col.low = sys->rng.low;
   vars.reg.row.high = vars.reg.col.high = sys->rng.high;
   spk1_reorder(&vars);
   return 0;
}

static int tranki_reorder(linsolqr_system_t sys,mtx_region_t *region)
/**
 ***  The region to reorder is first isolated by truncating the region
 ***  provided to the largest square region confined to the matrix diagonal.
 ***  It is presumed it will contain no empty rows or columns and will
 ***  provide the basis of candidate pivots when factoring.
 **/
{
   struct creorder_vars vars;
   CHECK_SYSTEM(sys);
   if ( !(sys->fclass==ranki || sys->fclass==s_qr) ) {
     FPRINTF(stderr,"linsolqr_reorder:  tspk1 reorder called on system\n");
     FPRINTF(stderr,"                   with inappropriate factor method\n");
     return 1;
   }
   if( region == mtx_ENTIRE_MATRIX ) determine_pivot_range(sys);
   else square_region(sys,region);

   vars.mtx = sys->coef;
   vars.reg.row.low = vars.reg.col.low = sys->rng.low;
   vars.reg.row.high = vars.reg.col.high = sys->rng.high;
   tspk1_reorder(&vars);
   return 0;
}

/***************************************************************************\
  End of reordering functions for SPK1.
\***************************************************************************/

/***************************************************************************\
  RANKI implementation functions.
\***************************************************************************/

static void eliminate_row(mtx_matrix_t mtx,
                          mtx_range_t *rng,
                          int32 row,      /* row to eliminate */
                          real64 *tmp,    /* temporary array */
                          real64 *pivots) /* prior pivots array */
/**
 ***  Eliminates the given row to the left of the diagonal element, assuming
 ***  valid pivots for all of the diagonal elements above it (the elements
 ***  above those diagonal elements, if any exist, are assumed to be U
 ***  elements and therefore ignored).  The temporary array is used by this
 ***  function to do its job.  tmp[k], where rng->low <= k <= rng->high must
 ***  be defined (allocated) but need not be initialized.
 ***  pivots[k] where rng->low <= k <=rng->high must be allocated and
 ***  pivots[rng->low] to pivots[row-1] must contain the previous pivots.
 **/
{
   mtx_coord_t nz;
   mtx_range_t left_to_eliminate,high_cols,tmp_cols;
   boolean do_series=FALSE;

   high_cols.low = row;
   high_cols.high = rng->high;
   left_to_eliminate.low = rng->low;
   left_to_eliminate.high = row - 1;

   /* Move non-zeros left of pivot from matrix to full array */
   mtx_zero_real64(tmp+rng->low,(row-rng->low));
   nz.row = row;
   mtx_cur_row_vec(mtx,row,tmp,&left_to_eliminate);
   mtx_clear_row(mtx,row,&left_to_eliminate);

   /* eliminates nzs from pivot, one by one, filling tmp with multipliers */
   do_series= ( (row-1) >= (rng->low) );
   if (do_series) {
     mtx_add_row_series_init(mtx,row,FALSE);
     /* do nothing between here and the END call to _init which removes
        elements from row */
   }
   for( nz.row = row-1 ; nz.row >= rng->low ; --(nz.row) ) {
      if( tmp[nz.row] == D_ZERO )
         continue;   /* Nothing to do for this row */

      nz.col = nz.row;
      tmp[nz.row] /= pivots[nz.row];
      /* tmp[nz.row] now equals multiplier */

      /* Perform "mtx_add_row" for full array part of the row */
      left_to_eliminate.high = nz.row - 1;
      mtx_cur_vec_add_row(mtx,tmp,nz.row,-tmp[nz.row],
                          &left_to_eliminate,FALSE);

      /* Perform "mtx_add_row" on remaining part of row */
      mtx_add_row_series(nz.row,-tmp[nz.row],&high_cols);
   }
   if (do_series) {
      mtx_add_row_series_init(mtx,mtx_NONE,TRUE);
   }

   mtx_fill_cur_row_vec(mtx,row,tmp,mtx_range(&tmp_cols,rng->low,row-1));

}


/* not in use it seems */
#if BUILD_DEAD_CODE /* yes, baa */
static int32 top_of_spike(linsolqr_system_t sys,int32 col)
/**
 ***  Determines the top row (row of lowest index) in a possible spike
 ***  above the diagonal element in the given column.  If there is no spike,
 ***  then (row = ) col is returned.
 **/
{
   mtx_range_t above_diagonal;
   mtx_coord_t nz;
   real64 value;
   int32 top_row;

   above_diagonal.low = sys->rng.low;
   above_diagonal.high = col-1;
   top_row = nz.col = col;
   nz.row = mtx_FIRST;
   while( value = mtx_next_in_col(sys->factors,&nz,&above_diagonal),
     nz.row != mtx_LAST )
      if( nz.row < top_row )
         top_row = nz.row;
   return( top_row );
}
#endif /* BUILD_DEAD_CODE */

static boolean col_is_a_spike(linsolqr_system_t sys,int32 col)
/**
 ***  Determines if the col is a spike, characterized by having any
 ***  nonzeros above the diagonal. By construction there shouldn't
 ***  be any numeric 0 nonzeros above the diagonal, so we don't prune
 ***  them out here.
 **/
{
  mtx_range_t above_diagonal;
  mtx_coord_t nz;

  above_diagonal.low = sys->rng.low;
  above_diagonal.high = col-1;
  nz.col = col;
  nz.row = mtx_FIRST;
  /* this loop is cheaper than it looks */
  while( (void)mtx_next_in_col(sys->factors,&nz,&above_diagonal),
         nz.row != mtx_LAST ) {
    if( nz.row < col ) return TRUE;
  }
  return( FALSE );
}

/* see implementation notes before modifying this function! */
static void number_drag(real64 *vec, int32 rfrom, int32 rto)
{
  int32 i;
  real64 tmp;
  if (rto <rfrom) {
    tmp=vec[rfrom];
    for (i=rfrom; i> rto; i--) {
      vec[i]=vec[i-1]; /* rotate right */
    }
    vec[rto]=tmp;
    return;
  }
  if (rto > rfrom) {
    tmp=vec[rfrom];
    for (i=rfrom; i< rto; i++) {
      vec[i]=vec[i+1]; /* rotate left */
    }
    vec[rto]=tmp;
  }
}

static void rankikw_factor(linsolqr_system_t sys)
/**
 ***  (the following description also applies to rankijz_factor which is
 ***  different only in pivot selection strategy.)
 ***  This is the heart of the linear equation solver.  This function
 ***  factorizes the matrix into a lower (L) and upper (U) triangular
 ***  matrix.  sys->smallest_pivot and sys->rank are calculated.  The
 ***  RANKI method is utilized.  At the end of elimination, the matrix A
 ***  is factored into A = U L, where U and L are stored as follows:
 ***
 ***      <----- r ----> <- n-r ->
 ***     +--------------+---------+
 ***     |              |         |
 ***     |         U    |         |
 ***     |              |         |
 ***     |   L          |         | r
 ***     |              |         |
 ***     +--------------+---------+
 ***     |              |         |
 ***     |              |   0     | n-r
 ***     |              |         |
 ***     +--------------+---------+
 ***
 ***  The rows and columns have been permuted so that all of the pivoted
 ***  original rows and columns are found in the first r rows and columns
 ***  of the region.  The last n-r rows and columns are unpivoted.  U has
 ***  1's on its diagonal, whereas L's diagonal is stored in the matrix.
 ***
 ***  Everything above was written as though the entire matrix is
 ***  involved.  In reality, only the relevant square region is involved.
 **/
{
   mtx_coord_t nz;
   int32 last_row;
   mtx_range_t pivot_candidates;
   real64 *tmp;
   real64 pivot, *pivots;
   int32 length;
   mtx_matrix_t mtx;

   length = sys->rng.high - sys->rng.low + 1;
   tmp = sys->ludata->tmp;
   /* eliminate row takes care of zeroing the relevant region and won't
      change values outside of it. */
   pivots=sys->ludata->pivlist;
   mtx=sys->factors;

   sys->smallest_pivot = MAXDOUBLE;
   last_row = pivot_candidates.high = sys->rng.high;
   for( nz.row = sys->rng.low ; nz.row <= last_row ; ) {

      pivot_candidates.low = nz.col = nz.row;
      pivots[nz.row]=pivot = mtx_value(mtx,&nz);
      pivot = fabs(pivot);
      if( pivot > sys->pivot_zero &&
     pivot >= sys->ptol * mtx_row_max(mtx,&nz,&pivot_candidates,NULL) &&
     !col_is_a_spike(sys,nz.row) ) {
     /* Good pivot and not a spike: continue with next row */
         if( pivot < sys->smallest_pivot )
            sys->smallest_pivot = pivot;
         ++(nz.row);
         continue;
      }
      /* pivots for rows nz.row back to sys->rng->low are stored in pivots */
      /**
       ***  Row is a spike row or will
       ***  be when a necessary column
       ***  exchange occurs.
       **/
      eliminate_row(mtx,&(sys->rng),nz.row,tmp,pivots);
      /* pivot will be largest of those available. get size and value */
      pivot=mtx_row_max(mtx,&nz,&pivot_candidates,&(pivots[nz.row]));
      if( pivot <= sys->pivot_zero ) { /* pivot is an epsilon */
         /* Dependent row, drag to the end */
         mtx_drag(mtx,nz.row,last_row);
         number_drag(pivots,nz.row,last_row);
         --last_row;
#undef KAA_DEBUG
#ifdef KAA_DEBUG
     FPRINTF(stderr,"Warning: Row %d is dependent with pivot %20.8g\n",
         nz.row,pivot);
#endif /* KAA_DEBUG */
      } else {
         /* Independent row: nz contains best pivot */
     /* Move pivot to diagonal */
         mtx_swap_cols(mtx,nz.row,nz.col);
     mtx_drag( mtx , nz.row , sys->rng.low );
         number_drag(pivots,nz.row,sys->rng.low);
         if( pivot < sys->smallest_pivot )
            sys->smallest_pivot = pivot;
         ++(nz.row);
      }
   }

   sys->rank = last_row - sys->rng.low + 1;
}


static void rankijz_factor(linsolqr_system_t sys)
/**
 ***  This is an alternate pivoting strategy introduced by Joe Zaher.
 ***  it looks down and across for good pivots rather than just across.
 ***
 **/
{
   real64 biggest;
   mtx_coord_t nz, best;
   mtx_region_t candidates;
   real64 *tmp;
   real64 pivot, *pivots;
   int32 length;
   mtx_matrix_t mtx;

   length = sys->rng.high - sys->rng.low + 1;
   tmp = sys->ludata->tmp;
   /* eliminate row takes care of zeroing the relevant region and won't
      change values outside of it. */
   pivots=sys->ludata->pivlist;
   mtx=sys->factors;

   sys->smallest_pivot = MAXDOUBLE;
   candidates.row.high = sys->rng.high;
   candidates.col.high = sys->rng.high;
   for( nz.row = sys->rng.low ; nz.row <= candidates.row.high ; ) {
      nz.col = nz.row;
      pivots[nz.row] = pivot = mtx_value(mtx,&nz);
      pivot = fabs(pivot);
      candidates.row.low = nz.row;
      candidates.col.low = nz.row;
      if( !col_is_a_spike(sys,nz.row) ) {
         best.col = nz.row;
         mtx_col_max(mtx,&best,&(candidates.row),&biggest);
         if( fabs(biggest) >= sys->pivot_zero ) {
            if( pivot < sys->pivot_zero || pivot < sys->ptol*fabs(biggest) ) {
               mtx_swap_rows(mtx,nz.row,best.row);
               pivots[nz.row] = biggest;
               pivot = fabs(biggest);
            }
            if( pivot < sys->smallest_pivot ) sys->smallest_pivot = pivot;
            ++(nz.row);
            continue;
         }
      }
      /* pivots for rows nz.row back to sys->rng->low are stored in pivots */
      /**
       ***  Row is a spike row or will
       ***  be when a necessary column
       ***  exchange occurs.
       **/
      eliminate_row(mtx,&(sys->rng),nz.row,tmp,pivots);
      /* pivot will be largest of those available. get size and value */
      pivot=mtx_row_max(mtx,&nz,&(candidates.col),&(pivots[nz.row]));
      /* Move pivot to diagonal */
      if( pivot < sys->pivot_zero ) { /* pivot is an epsilon */
         /* Dependent row, drag nz to lower right */
         mtx_drag(mtx,nz.row,candidates.row.high);
         number_drag(pivots,nz.row,candidates.row.high);
         --(candidates.row.high);
      } else {
         /* Independent row, drag nz to upper left */
         mtx_swap_cols(mtx,nz.row,nz.col);
         mtx_drag( mtx , nz.row , sys->rng.low );
         number_drag(pivots,nz.row,sys->rng.low);
         if( pivot < sys->smallest_pivot )
            sys->smallest_pivot = pivot;
         ++(nz.row);
      }
   }

   sys->rank = candidates.row.high - sys->rng.low + 1;
}


#define KAA_DEBUG 0
#if KAA_DEBUG
int ranki_entry(linsolqr_system_t sys,mtx_region_t *region)
#else
static int ranki_entry(linsolqr_system_t sys,mtx_region_t *region)
#endif
/**
 ***  The region to factor is first isolated by truncating the region
 ***  provided to the largest square region confined to the matrix diagonal.
 ***  It is presumed it will contain no empty rows or columns and that it has
 ***  been previously reordered using linsolqr_reorder(sys,region,spk1)
 ***  or a sane variant of spk1.
 ***  This is the entry point for all ranki based strategies, regardless of
 ***  pivot selection subtleties.
 **/
{
   struct rhs_list *rl;
   double time;

   CHECK_SYSTEM(sys);
   if( sys->factored )
      return 0;
   switch(sys->fmethod) {
   case ranki_kw:
   case ranki_jz:
   case ranki_ka: /* add new methods here */
     break;
   default:
     return 1;
   }

   if(ISNULL(sys->ludata))
      return 1;

   if( NOTNULL(sys->inverse) )
      mtx_destroy(sys->inverse);
   sys->inverse=NULL;
   if( NOTNULL(sys->factors) )
      mtx_destroy(sys->factors);
   if( region == mtx_ENTIRE_MATRIX ) determine_pivot_range(sys);
   else square_region(sys,region);

   sys->factors = mtx_copy_region(sys->coef,region);
   sys->rank = -1;
   sys->smallest_pivot = MAXDOUBLE;
   for( rl = sys->rl ; NOTNULL(rl) ; rl = rl->next )
      rl->solved = FALSE;
   insure_capacity(sys);
   insure_lu_capacity(sys);

   time = tm_cpu_time();
   switch(sys->fmethod) {
   case ranki_ka:
   case ranki_kw:
     rankikw_factor(sys);
     break;
   case ranki_jz:
   default:
     rankijz_factor(sys);
   }
 /* no longer done automatically as noone usually cares.
   calc_dependent_rows_ranki1(sys);
   calc_dependent_cols_ranki1(sys);
 */
   sys->factored = TRUE;

#undef KAA_DEBUG
#if KAA_DEBUG
   time = tm_cpu_time() - time;
   FPRINTF(stderr,"Time for Inversion = %f\n",time);
   FPRINTF(stderr,"Non-zeros in Inverse = %d\n",
       mtx_nonzeros_in_region(sys->factors,region));
#endif
   return 0;
}

static int ranki_solve(linsolqr_system_t sys, struct rhs_list *rl)
{
   backward_substitute(sys,rl->varvalue,rl->transpose);
   forward_substitute(sys,rl->varvalue,rl->transpose);
   zero_unpivoted_vars(sys,rl->varvalue,rl->transpose);
   return 0;
}


/*
 *********************************************************************
 * Start of the 2 bodied factorization schemes.
 *
 * ranki2_entry will be used to access these routines.
 * there is now code called:
 * rankikw2_factor
 * rankijz2_factor
 * eliminate_row2
 * ranki2_entry
 * forward_substitute2
 * backward_substitute2
 * ranki2_solve
 *
 * We have fixed the calc_col_dependent and calc_row_dependent
 * routines.
 *********************************************************************
 */

extern int32 mtx_number_ops;

/* this function does no permutation of any sort */
static
void eliminate_row2(mtx_matrix_t mtx,
           mtx_matrix_t upper_mtx,
           mtx_range_t *rng,
           int32 row,       /* row to eliminate */
           real64 *tmp,     /* temporary array */
           real64 *pivots,  /* prior pivots array */
                   real64 dtol      /* drop tolerance */
                   )
{
  mtx_coord_t nz;
  int j,low,high;
  double tmpval;

  /*
   * Move all non-zeros from current row to full array.
   * The full array would have been initialized before,
   * we must put it back in the clean state when we leave.
   * All operations are done over mtx_ALL_COLS.
   *
   * Note: because of an addrow over mtx_ALL_COLS, entries
   * of tmp outside rng may have nonzeros put in them if
   * sys->factors has nonzeros in the outlying columns.
   * If enough of these pile up out there during elimination
   * of a block and the compiler treats double overflow as
   * a signalable error, we will have a floating point
   * exception.
   * Currently sys->factors is a copy of a region from
   * the sys->coef matrix, so we should not have anything
   * out there to pile up. If we make a variant which
   * does not copy the coef matrix but uses it directly,
   * this code needs to be revisited.
   */
  mtx_steal_cur_row_vec(mtx,row,tmp,mtx_ALL_COLS);

  /*
   * Eliminates nzs from pivot, one by one, filling tmp with multipliers
   * We now operate on the entire row, since we are not writing the
   * multipliers back to the matrix.
   */

  low = rng->low;
  high = rng->high;
  for (j = row-1 ; j >= low ; --(j) ) {
    if (tmp[j] == D_ZERO)
      continue;         /* Nothing to do for this row */
    tmpval = tmp[j]/pivots[j];  /* Compute multiplier */

    /*
     * tmpval is now the multiplier. We use it to eliminate the row,
     * but backpatch it to tmp, *after* we do the elimination, as
     * mtx_cur_vec_add_row over all columns will stomp on tmp[j]
     */
    mtx_cur_vec_add_row(mtx,tmp,j,-tmpval,mtx_ALL_COLS,FALSE);
    tmp[j] = tmpval;        /* patch the diagonal */
  }

  /*
   * Fill up the upper triangular matrix.
   * refill the range [low .. row-1]. Remember that we have
   * to zero all nnz's in tmp.
   */
  nz.row = row;
  for (j=low;j<row;j++) {
    if (tmp[j] != D_ZERO) {
      nz.col = j;
      mtx_fill_value(upper_mtx,&nz,tmp[j]);
#if RBADEBUG
      FPRINTF(gscr,"fillingU: or %d oc %d (cc %d) %24.18g\n",
        mtx_row_to_org(mtx,nz.row), mtx_col_to_org(mtx,nz.col), nz.col, tmp[j]);
#endif
      tmp[j] = D_ZERO;
    }
  }
  /*
   * refill the range [row .. high]. We *wont* keep
   * the diagonal as this is sorted out by the pivot
   * array, but here we must keep it, sigh.
   * At this stage we don't know that the diagonal is the
   * pivot, as that is determined after elimination done.
   * Outer loop must delete the ultimate pivot element.
   * Odds are good, however, that it is the current diagonal,
   * so put in the diagonal last.
   */
  