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.
   */
  low = row;
  for (j=high;j>=low;j--) {
    if (tmp[j] != D_ZERO) {
      if (fabs(tmp[j]) > dtol) {
    nz.col = j;
    mtx_fill_value(mtx,&nz,tmp[j]);
#if RBADEBUG
        FPRINTF(gscr,"fillingL: or %d oc %d %24.18g\n",
          mtx_row_to_org(mtx,nz.row), mtx_col_to_org(mtx,nz.col), tmp[j]);
#endif
      }
      tmp[j] = D_ZERO;      /* zero element regardless */
    }
  }
  tmp[row] = D_ZERO;
}

/*
 * This function is the same as rankikw_factor except
 * that it uses 2 matrices; one to store L and one for U.
 * As such it uses eliminate_row2 rather than eliminate_row.
 * Note there has been a change in the way ptol is used,
 * now selecting nearest passing, rather than taking max.
 * We assume there is NO incidence outside the region to be factored
 * in sys->factors.
 */
static void rankikw2_factor(linsolqr_system_t sys)
{
  mtx_coord_t nz;
  int32 last_row, defect;
  mtx_range_t pivot_candidates;
  real64 *tmp;
  real64 pivot, *pivots;
  mtx_matrix_t mtx, upper_mtx;


  /**
   ** tmp and pivots are cur indexed, tmp by cols, pivots by rows.
   ** rather than allocating tmp every time, which can get expensive,
   ** we allocate it once *elsewhere* and zero it all here.
   ** Because we don't know the state of it on entry (due to exceptions,
   ** etc, the rezeroing is unavoidable.
   ** Eliminate_row2 is responsible for returning the working region
   ** zeroed.
   **/

  tmp = sys->ludata->tmp;
  mtx_zero_real64(tmp,sys->capacity);
  pivots = sys->ludata->pivlist;
  mtx = sys->factors;
  upper_mtx = sys->inverse;
  defect = 0;
#if RBADEBUG
  gscr = fopen("/usr1/ballan/tmp/lu/rkw","w+");
#endif

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

#if (RBADEBUG && 0)
    /* remove the && 0 and you better have a lot for free
     * disk space
     */
    mtx_write_region_human_orgrows(gscr,mtx,mtx_ENTIRE_MATRIX);
#endif

    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;
      }
#if RBADEBUG
      FPRINTF(gscr,"Cheap pivot col %d (org %d)\n",
        nz.row, mtx_col_to_org(mtx,nz.row));
      if (nz.row ==179) {
        FPRINTF(stderr,"ben, stop here\n");
      }
#endif
      ++(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_row2(mtx,upper_mtx,&(sys->rng),nz.row,tmp,pivots,sys->dtol);
    /* pivot will be leftmost of those that pass ptol and eps, or 0.0. */
    if (defect) {
      /* unfortunately, can't use mtx_ALL_COLS in search now because we must
         ignore columns that were dragged out with previous singular rows. */
      pivot = mtx_get_pivot_col(mtx,&nz,&pivot_candidates,&(pivots[nz.row]),
                                sys->ptol,sys->pivot_zero);
    } else {
      /* no singular cols, and eliminate moved all left of diag crap out, so */
      pivot = mtx_get_pivot_col(mtx,&nz,mtx_ALL_COLS,&(pivots[nz.row]),
                                sys->ptol,sys->pivot_zero);
    }
#if RBADEBUG
      FPRINTF(gscr,"Pivot column found %d (org %d)\n",
        nz.col, mtx_col_to_org(mtx,nz.col));
#endif
    if( pivot < sys->pivot_zero ) { /* pivot is an epsilon */
      /*
       * Dependent row, drag to the end. The upper_mtx is a slave
       * of mtx, and will be dragged automatically.
       */
      mtx_drag(mtx, nz.row, last_row);
      number_drag(pivots, nz.row, last_row);
      --last_row;
      defect = 1;
    } else {
      /* Independent row: nz contains selected pivot */
      mtx_swap_cols(mtx, nz.row, nz.col);   /* this Fixes U as well */
      /* Move pivot to diagonal */
      mtx_drag(mtx, nz.row, sys->rng.low ); /* this Fix U as well */
      number_drag(pivots, nz.row, sys->rng.low);
      if( pivot < sys->smallest_pivot )
    sys->smallest_pivot = pivot;
      ++(nz.row);
    }
  }
#if RBADEBUG
/* get rid of this  asap*/
  {
    FILE *fp;
    int32 cc;

    fp = fopen("/usr1/ballan/tmp/lu/kw2p","w+");
    FPRINTF(fp,"kw2 final pivot sequence:\n");
    for (cc= sys->rng.low; cc <= last_row; cc++) {
      FPRINTF(fp,"orgrow,orgcol = %d,%d,%24.18g\n",
        mtx_row_to_org(mtx,cc),
        mtx_col_to_org(mtx,cc), pivots[cc]);
    }
    fclose(fp);
  }

#endif

  zero_diagonal_elements(sys->factors,sys->rng.low,sys->rng.high);
  sys->rank = last_row - sys->rng.low + 1;

#if RBADEBUG
  {
   FILE *fp;
   fp = fopen("/usr1/ballan/tmp/lu/kw2m","w+");
   mtx_write_region_human_rows(fp,mtx,mtx_ENTIRE_MATRIX);
   fclose(fp);
   fp = fopen("/usr1/ballan/tmp/lu/kw2um","w+");
   mtx_write_region_human_rows(fp,upper_mtx,mtx_ENTIRE_MATRIX);
   fclose(fp);
  }
#endif
#if RBADEBUG
  fclose(gscr);
#endif
}

/*
 * This function is the same as rankijz_factor except
 * that it uses 2 matrices; one to store L and one for U.
 * As such it uses eliminate_row2 rather than eliminate_row.
 *
 * We are now using the slave matrix feature of the code, and
 * thus we all drag/swaps automatically affect the slave (upper_mtx)
 * as well.
 */
static void rankijz2_factor(linsolqr_system_t sys)
{
  real64 biggest;
  mtx_coord_t nz, best;
  mtx_region_t candidates;
  real64 *tmp;
  real64 pivot, *pivots;
  mtx_matrix_t mtx, upper_mtx;
  /**
   ** tmp and pivots are cur indexed, tmp by cols, pivots by rows.
   ** rather than allocating tmp every time, which can get expensive,
   ** we allocate it once elsewhere and zero it all here.
   ** Because we don't know the state of it on entry (due to exceptions,
   ** etc), the rezeroing is unavoidable.
   ** Eliminate_row2 is responsible for returning the working region
   ** zeroed.
   **/
  tmp = sys->ludata->tmp;
  mtx_zero_real64(tmp,sys->capacity);
  pivots=sys->ludata->pivlist;
  mtx = sys->factors;
  upper_mtx = sys->inverse;

  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);
      /* make sure we aren't missing someone really nice in the column */
      if( fabs(biggest) >= sys->pivot_zero ) {
    if( pivot < sys->pivot_zero || pivot < sys->ptol*fabs(biggest) ) {
      mtx_swap_rows(mtx,nz.row,best.row); /* Fixes U as well */
      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_row2(mtx,upper_mtx,&(sys->rng),nz.row,tmp,pivots,sys->dtol);
    /* 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);     /* this Fixes U as well */
      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);
    }
  }
  zero_diagonal_elements(sys->factors,sys->rng.low,sys->rng.high);
  sys->rank = candidates.row.high - sys->rng.low + 1;
}

/*
 * See the comments attached to forward_substitute.
 * This code otherwise is the same. It uses a 2 bodied
 * matrix as well as making use of mtx_ALL_COLS and
 * mtx_ALL_ROWS whereever possible.
 * We make the following additional assumptions here:
 *   (if they do not hold, do NOT use this function)
 * - sys->inverse (U) has no incidence that is not multipliers
 *   (and the diagonal of 1s is NOT in sys->inverse.)
 *   As of 10/95, this is how ranki2 U is constructed.
 * - sys->factors (L) has no incidence on the upper triangle,
 *   including the diagonal, or outside the factored region.
 *    relaxation: incidence anywhere allowed if value = 0.0
 *                since 0 doesn't contribute to a dot product
 *                and the only thing we do with triangles is dot them.
 * - There may be singular rows and columns in the factorization,
 *   but any additions coming from these rows/columns during
 *   mtx_ALL_*O*S operations will not contribute to sums because the
 *   user zeroed the arr entries corresponding to these before
 *   calling this function.
 */
static
void forward_substitute2(linsolqr_system_t sys,
                real64 *arr,
                boolean transpose)
{
  mtx_coord_t nz;
  real64 sum, *pivlist;
  mtx_matrix_t mtx;
  int32 dotlim;
  boolean nonzero_found=FALSE;

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

      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,mtx_ALL_ROWS,TRUE);
    /* arr[org_col] = (arr[org_col] - sum)  / D_ONE */;
    arr[org_col] -= sum;
      }
    }
  } else {          /* arr is indexed by original row number */
    mtx=sys->factors;
    for( nz.row=sys->rng.low; nz.row < dotlim; ++(nz.row) ) {
      register int32 org_row;

      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,mtx_ALL_COLS,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];
      }
    }
  }
}

/*
 * See the comments attached to backward_substitute and
 * forward_substitute2.
 * When solving for the transpose, then we are actually
 * running of the lower triangle, hence we use sys->factors.
 * Otherwise we use sys->inverse which stores U.
 */
static
void backward_substitute2(linsolqr_system_t sys,
                 real64 *arr,
                 boolean transpose)
{
  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 */

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

      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,mtx_ALL_ROWS,TRUE);
    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 */
    mtx = sys->inverse;
    for( nz.row = sys->rng.low+sys->rank-1; nz.row >= dotlim ; --(nz.row) ) {
      register int32 org_row;

      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,mtx_ALL_COLS,TRUE);
    arr[org_row] -= sum;
      }
    }
  }
}

static void calc_dependent_rows_ranki2(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));
    /* must zero the whole thing to use the backward_substitute2 right */
    mtx_org_row_vec(mtx,nz.row,lc,&colrange);
    if( nz.row < rowrange.high || nz.row > rowrange.low )
      backward_substitute2(sys,lc,TRUE);
    forward_substitute2(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_ranki2(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_substitute2(sys,lc,FALSE);
    forward_substitute2(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);
    }
  }
}

static int ranki2_solve(linsolqr_system_t sys, struct rhs_list *rl)
{
  /* zero any unsolved for vars first so they don't contaminate
     mtx_ALL_*O*S dot products.
  */
  zero_unpivoted_vars(sys,rl->varvalue,rl->transpose);
  backward_substitute2(sys,rl->varvalue,rl->transpose);
  forward_substitute2(sys,rl->varvalue,rl->transpose);
  return 0;
}

/* include rankiba2 implementation. Note that this is VERY VERY bad style
and needs to be fixed immediately by splitting up the linsolqr file.
This is our fastest (far and away) ranki implementation -- drag free.
*/
#if 1
#include "rankiba2.c"
#endif

#ifdef BUILD_KIRK_CODE
/*
 * this forward declaration is necessary to keept the
 * system happy until we can organize this file
 */

static
int kirk1_factor(linsolqr_system_t sys,
          mtx_region_t *A11,
          int kirk_method);

#endif /* BUILD_KIRK_CODE */

/*
 * This is the entry point for all the ranki2_* schemes
 * which make use of a 2 bodied matrix.
 * it does not respond to ranki_kw and ranki_jz as if they
 * were calls for ranki_kw2 and ranki_jz2. with the new
 * factorization classes , there's really no call to be
 * sneaking around the header functions.
 * ranki_ka goes here as well. check sys->fmethod to see whos
 * expected to be in use.
 */
int ranki2_entry(linsolqr_system_t sys, mtx_region_t *region)
{
  struct rhs_list *rl;
  double time;

  CHECK_SYSTEM(sys);
  if (sys->factored) return 0;
  switch(sys->fmethod) {
  case ranki_ba2:
  case ranki_kw2:
  case ranki_jz2:    /* add new methods here */
#ifdef BUILD_KIRK_CODE
  case ranki_ka:
#endif /* BUILD_KIRK_CODE */
    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->inverse = mtx_create_slave(sys->factors);
  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_ba2:
    rankiba2_factor(sys);
    break;
  case ranki_kw2:
    rankikw2_factor(sys);
    break;
  case ranki_jz2:
    rankijz2_factor(sys);
    break;

#ifdef BUILD_KIRK_CODE
  case ranki_ka:
    kirk1_factor(sys,region,2);
    break;
#endif /* BUILD_KIRK_CODE */

  default:
    return 1;
  }
  sys->factored = TRUE;

#define KAA_DEBUG 1
#if KAA_DEBUG
  if (g_linsolqr_timing) {
    int anz;
    int fnz;
    anz = mtx_nonzeros_in_region(sys->coef,region);
    fnz = mtx_nonzeros_in_region(sys->factors,region) +
      mtx_nonzeros_in_region(sys->inverse,0);
    time = tm_cpu_time() - time;
    FPRINTF(stderr,"A-NNZ: %d Factor time: %f Fill %g\n",
      anz,time,( anz>0 ? (double)fnz/(double)anz : 0));
  }
#endif /* KAA_DEBUG */
#undef KAA_DEBUG
  return 0;
}

/***************************************************************************\
  End of RANKI implementation functions.
\***************************************************************************/


#ifdef BUILD_KIRK_CODE
/*
 ***************************************************************************
 * Start of kirk_* routines
 *
 * kirk1_factor:
 *  This routine is  based on rankikw2_factor, except that it takes
 *  an additional region, so as to do some *global* pivot restriction.
 *  This is for the case of a single border for the entire matrix.
 *
 ***************************************************************************
 */

/*
 * This is just a dummy structure so as to get the size
 * correct. It is *exactly the same size as mtx_linklist,
 * until we can sort out the software engineering issues.
 */
struct dlinklist {
  int index;
  struct dlinklist *prev;
};

static
int kirk1_factor(linsolqr_system_t sys,
          mtx_region_t *A11,
          int kirk_method)
{
  mtx_coord_t nz;
  int32 last_row;
  mtx_range_t pivot_candidates;
  real64 *tmp,*tmp_array_origin;
  real64 pivot, *pivots;
  int32 length;
  mtx_matrix_t mtx, upper_mtx;
  real64 maxa;

  int32 *inserted = NULL;       /* stuff for the link list */
  mem_store_t eltbuffer = NULL;

#ifdef NOP_DEBUG
  mtx_number_ops = 0;
#endif /* NOP_DEBUG */

  length = sys->rng.high + 1;
  tmp_array_origin = length > 0 ?
    (real64 *)asccalloc( length,sizeof(real64) ) : NULL;
  tmp = tmp_array_origin;
  pivots=sys->ludata->pivlist;
  mtx = sys->factors;
  upper_mtx = sys->inverse;
  inserted = (int32 *)asccalloc(length,sizeof(int32));
  eltbuffer = mem_create_store(2,256,sizeof(struct dlinklist),
                   2,256);

  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);
    maxa = mtx_row_max(mtx,&nz,&pivot_candidates,NULL);
    if ((pivot > sys->pivot_zero) && (pivot >= sys->ptol*maxa) &&
    !col_is_a_spike(sys,nz.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.
     **/
    if (kirk_method==1)
      eliminate_row2(mtx,upper_mtx,&(sys->rng),nz.row,tmp,pivots,sys->dtol);
    else{
      mtx_eliminate_row2(mtx,upper_mtx,&(sys->rng),
             nz.row,tmp,pivots,inserted,eltbuffer);
      mem_clear_store(eltbuffer);
    }
    /* 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. The upper_mtx is a slave
       * of mtx, and will be dragged automatically.
       */
      mtx_drag(mtx,nz.row,last_row);
      number_drag(pivots,nz.row,last_row);
      --last_row;
#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 */
      mtx_swap_cols(mtx,nz.row,nz.col);     /* this Fixes U as well */