ascend4/solver/mtx_basic.c

/*
 *  mtx: Ascend Sparse Matrix Package
 *  by Karl Michael Westerberg
 *  Created: 2/6/90
 *  Version: $Revision: 1.20 $
 *  Version control file: $RCSfile: mtx_basic.c,v $
 *  Date last modified: $Date: 2000/01/25 02:27:07 $
 *  Last modified by: $Author: ballan $
 *
 *  This file is part of the SLV solver.
 *
 *  Copyright (C) 1990 Karl Michael Westerberg
 *  Copyright (C) 1993 Joseph Zaher
 *  Copyright (C) 1994 Joseph Zaher, Benjamin Andrew Allan
 *  Copyright (C) 1995 Kirk Andre Abbott, Benjamin Andrew Allan
 *  Copyright (C) 1996 Benjamin Andrew Allan
 *
 *  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.
 */

#include <math.h>
#include "utilities/ascConfig.h"
#include "utilities/ascMalloc.h"
#include "utilities/mem.h"
#include "solver/mtx.h"
/* grab our private parts */
#define __MTX_C_SEEN__
#include "solver/mtx_use_only.h"

/* this returns an error count */
int super_check_matrix( mtx_matrix_t mtx)
/**
 ***  Really check the matrix.
 **/
{
  int32 ndx,errcnt=0;
  int32 rowc,colc;

  /* Test consistency of permutation arrays */
  for( ndx=ZERO ; ndx < mtx->capacity ; ++ndx ) {
    if( mtx->perm.row.org_to_cur[mtx->perm.row.cur_to_org[ndx]] != ndx ) {
      FPRINTF(g_mtxerr,"ERROR:  (mtx) super_check_matrix\n");
      FPRINTF(g_mtxerr,"        Permutation violation in row %d.\n", ndx);
        errcnt++;
    }
    if( mtx->perm.col.org_to_cur[mtx->perm.col.cur_to_org[ndx]] != ndx ) {
      FPRINTF(g_mtxerr,"ERROR:  (mtx) super_check_matrix\n");
      FPRINTF(g_mtxerr,"        Permutation violation in col %d.\n",ndx);
      errcnt++;
    }
  }

  if( mtx->order > mtx->capacity ) {
    FPRINTF(g_mtxerr,"ERROR:  (mtx) super_check_matrix\n");
    FPRINTF(g_mtxerr,"        Capacity %d is less than order %d\n",
      mtx->capacity,mtx->order);
    errcnt++;
  }

  /* Make sure rows and columns which don't exist are blank */
  for( ndx = mtx->order ; ndx < mtx->capacity ; ++ndx ) {
    int32 org_row,org_col;
    org_row = mtx->perm.row.cur_to_org[ndx];
    org_col = mtx->perm.col.cur_to_org[ndx];
    if( NOTNULL(mtx->hdr.row[org_row]) ) {
      FPRINTF(g_mtxerr,"ERROR:  (mtx) super_check_matrix\n");
      FPRINTF(g_mtxerr,"        Non-zeros found in non-existent row %d.\n",ndx);
      errcnt++;
    }
    if( NOTNULL(mtx->hdr.col[org_col]) ) {
      FPRINTF(g_mtxerr,"ERROR:  (mtx) super_check_matrix\n");
      FPRINTF(g_mtxerr,"        Non-zeros found in non-existent col %d.\n",ndx);
      errcnt++;
    }
  }

  /* Verify fishnet structure */
  for( ndx=rowc=colc=ZERO ; ndx < mtx->capacity ; ++ndx ) {
    struct element_t *elt;
    int32 org_row,org_col;
    org_row = mtx->perm.row.cur_to_org[ndx];
    org_col = mtx->perm.col.cur_to_org[ndx];
    for( elt = mtx->hdr.row[org_row] ; NOTNULL(elt) ; elt = elt->next.col ) {
      if( elt->row != mtx_NONE ) {
        ++rowc;
        if( elt->row != org_row ) {
          FPRINTF(g_mtxerr,"ERROR:  (mtx) super_check_matrix\n");
          FPRINTF(g_mtxerr,"        Element mismatch in row %d.\n", ndx);
          errcnt++;
        }
      } else {
        FPRINTF(g_mtxerr,"ERROR:  (mtx) super_check_matrix\n");
        FPRINTF(g_mtxerr,"        Disclaimed element in row %d.\n", ndx);
        errcnt++;
      }
    }
    for( elt = mtx->hdr.col[org_col] ; NOTNULL(elt) ; elt = elt->next.row ) {
      if( elt->col != mtx_NONE ) {
        ++colc;
        if( elt->col != org_col ) {
          FPRINTF(g_mtxerr,"ERROR:  (mtx) super_check_matrix\n");
          FPRINTF(g_mtxerr,"        Element mismatch in col %d.\n", ndx);
          errcnt++;
        }
      } else {
        FPRINTF(g_mtxerr,"ERROR:  (mtx) super_check_matrix\n");
        FPRINTF(g_mtxerr,"        Disclaimed element in col %d.\n", ndx);
        errcnt++;
      }
    }
  }
  if( rowc != colc ) {
    FPRINTF(g_mtxerr,"ERROR:  (mtx) super_check_matrix\n");
    FPRINTF(g_mtxerr,"        Non-zero discrepancy, %d by row, %d by col.\n",
      rowc, colc);
    errcnt++;
  }
  return errcnt;
}

boolean check_matrix( mtx_matrix_t mtx, char *file, int line)
/**
 ***  Checks the integrity flag of the matrix.
 **/
{

   if( ISNULL(mtx) ) {
      FPRINTF(g_mtxerr,"ERROR:  (mtx) check_matrix\n");
      FPRINTF(g_mtxerr,"        NULL matrix in %s line %d.\n",file,line);
      return 0;
   }

   switch( mtx->integrity ) {
      case OK:
         break;
      case DESTROYED:
         FPRINTF(g_mtxerr,"ERROR:  (mtx) check_matrix\n");
         FPRINTF(g_mtxerr,
           "        Matrix deceased found in %s line %d.\n",file,line);
         return 0;
      default:
         FPRINTF(g_mtxerr,"ERROR:  (mtx) check_matrix\n");
         FPRINTF(g_mtxerr,"        Matrix garbage found in %s line %d.\n",
           file,line);
         return 0;
   }

#if MTX_DEBUG
   super_check_matrix(mtx);
#endif

   if (ISSLAVE(mtx)) {
     return (check_matrix(mtx->master,file,line));
   }
   return 1;
}

/* this function checks the consistency of a sparse as best it can.
   returns FALSE if something wierd. */
boolean check_sparse(const mtx_sparse_t * const sp, char *file, int line)
{
  if ( ISNULL(sp) ||
       sp->len > sp->cap ||
       ( sp->cap > 0 && ( ISNULL(sp->idata) || ISNULL(sp->data) ) )
     ) {
    FPRINTF(g_mtxerr,"ERROR:  (mtx) check_sparse - Inconsistent or\n");
    FPRINTF(g_mtxerr,"        NULL sparse in %s line %d.\n",file,line);
    return 0;
  }
  return 1;
}

#define free_unless_null(ptr) if( NOTNULL(ptr) ) ascfree(ptr)
/*
static void free_unless_null(POINTER ptr)
{
   if( NOTNULL(ptr) )
      ascfree(ptr);
}
*/

void mtx_zero_int32(int32 *data, int len)
{
  int i;
  if (data==NULL) return;
  for (i=0; i<len; i++) data[i]=0;
}

void mtx_zero_real64(real64 *data, int len)
{
  int i;
  if (data==NULL) return;
  for (i=0; i<len; i++) data[i]=0.0;
}

void mtx_zero_ptr(void **data, int len)
{
  int i;
  if (data==NULL) return;
  for (i=0; i<len; i++) data[i]=NULL;
}

static mtx_matrix_t alloc_header(void)
/**
 ***  Allocates a matrix header and returns a pointer to it.
 **/
{
   mtx_matrix_t mtx;
   mtx = (mtx_matrix_t)ascmalloc( sizeof(struct mtx_header) );
   mtx->integrity = OK;
   return(mtx);
}

static void free_header(mtx_matrix_t mtx)
/**
 ***  Frees a matrix header.
 **/
{
   mtx->integrity = DESTROYED;
   ascfree(mtx);
}

#ifdef DELETE_THIS_UNUSED_FUNCTION
/* not IN use? */
static struct element_t **alloc_nz_header(int32 cap)
/**
 ***  Allocates a row or column header randomly initialized.
 **/
{
   return(cap > 0 ? (struct element_t **)
          ascmalloc(cap*sizeof(struct element_t *)) : NULL);
}
#endif /* DELETE_THIS_UNUSED_FUNCTION */


static struct element_t **calloc_nz_header(int32 cap)
/**
 ***  Allocates a zeroed row or column header.
 **/
{
   return(cap > 0 ? (struct element_t **)
          asccalloc(cap,sizeof(struct element_t *)) : NULL);
}

static void copy_nz_header(struct element_t **tarhdr,
             struct element_t **srchdr, int32 cap)
/**
 ***  Copies srchdr to tarhdr given the capacity of srchdr.
 **/
{
   mem_copy_cast(srchdr,tarhdr,cap*sizeof(struct element_t *));
}

#define free_nz_header(hdr) free_unless_null( (POINTER)(hdr) )
/**
 ***  Frees a row or column header
 **/


/********************************************************************/
/*
Any method of deallocating elements that doesn't use
delete_from_row/col or nuke_fishnet is DOOMED to die.
On ANY entry to delete_from_row or delete_from_col,
last_value_matrix should be already set to the mtx being deleted from.
What this really means is that before any of the following list of functions
is called, the last_value_matrix should be set:

disclaim element
delete_from_row/col
blast_row/col
clean_row/col

Note that the above functions call each other, with blast_/clean_ being
the top housekeeping/garbage collection pair. Users of delete, blast or
clean should set last_value_matrix at the scope where the matrix is known.

Excuses for this kluge:
1) mtx_value/set_value need it.
2) the mem_store scheme REQUIRES it.
3) it is faster than passing the mtx
up and down through the functions listed above, and we do a LOT of
these operations.
*/
static mtx_matrix_t last_value_matrix = NULL;

static void disclaim_element(struct element_t *element)
/**
 ***  Indicates that one of the row or column (it doesn't matter) in which
 ***  this element appears will no longer point to this element (it has
 ***  already been removed from the list).  Elements disclaimed once have
 ***  their row,col indices set to mtx_NONE.  Elements disclaimed twice are
 ***  discarded.
 ***  NOTE WELL: last_value_mtx must be set before this function is called.
 **/
{
  if( element->row == mtx_NONE ) {
    mem_free_element(last_value_matrix->ms,(void *)element);
  } else {
    element->row = element->col = mtx_NONE;
  }
}

static void delete_from_row(struct element_t **link)
/**
 ***  Deletes the element from the row it is in, given a pointer to
 ***  the row link which leads to this element.
 ***  On return the pointer pointed to by link points to the new next element.
 ***  NOTE WELL: last_value_mtx must be set before this function is called.
 **/
{
   struct element_t *p;
   p = *link;
   *link = p->next.col;
   disclaim_element(p);
   /* conservatively cause mtx_set_value to forget */
   last_value_matrix->last_value = NULL;
}

static void delete_from_col(struct element_t **link)
/**
 ***  Deletes the element from the col it is in, given a pointer to
 ***  the col link which leads to this element.
 ***  NOTE WELL: last_value_mtx must be set before this function is called.
 **/
{
   struct element_t *p;
   p = *link;
   *link = p->next.row;
   disclaim_element(p);
   /* conservatively cause mtx_set_value to forget */
   last_value_matrix->last_value = NULL;
}

static void nuke_fishnet(mtx_matrix_t mtx)
/**
 ***  Deletes every element even remotely associated with a matrix
 ***  and nulls the headers. Also nulls the headers of slaves
 ***  and cleans out their elements.
 **/
{
  int32 i;

  if ( ISNULL(mtx) || ISNULL(mtx->hdr.row) ||
       ISNULL(mtx->hdr.col) || ISSLAVE(mtx)) {
    FPRINTF(g_mtxerr,"WARNING:     (mtx) nuke_fishnet called\n");
    FPRINTF(g_mtxerr,"                   with bad matrix.\n");
    return;
  }
  if (mtx->capacity<1) return;
  mtx->last_value = NULL;
  mem_clear_store(mtx->ms);
  zero(mtx->hdr.row,mtx->capacity,struct element_t *);
  zero(mtx->hdr.col,mtx->capacity,struct element_t *);
  for (i = 0; i < mtx->nslaves; i++) {
    zero(mtx->slaves[i]->hdr.row,mtx->capacity,struct element_t *);
    zero(mtx->slaves[i]->hdr.col,mtx->capacity,struct element_t *);
  }
}

static void blast_row( mtx_matrix_t mtx, int32 org, mtx_range_t *rng)
/**
 ***  Destroys all elements in the given row with (current)
 ***  col index in the given col range: if rng == mtx_ALL_COLS,
 ***  entire row is destroyed.
 ***  NOTE WELL: last_value_mtx must be set before this function is called.
 **/
{
  struct element_t **link;

  link = &(mtx->hdr.row[org]);
  if( rng == mtx_ALL_COLS ) {
    while( NOTNULL(*link) )
      delete_from_row(link);
  } else {
    int32 *tocur = mtx->perm.col.org_to_cur;
    int32 col;
    while( NOTNULL(*link) ) {
      col = (*link)->col;
      if( col == mtx_NONE || in_range(rng,tocur[col]) )
        delete_from_row(link);
      else
        link = &((*link)->next.col);
    }
  }
}

static void blast_col( mtx_matrix_t mtx, int32 org, mtx_range_t *rng)
/**
 ***  Destroys all elements in the given col with (current)
 ***  row index in the given row range: if rng == mtx_ALL_ROWS,
 ***  entire col is destroyed.
 ***  NOTE WELL: last_value_mtx must be set before this function is called.
 **/
{
  struct element_t **link;

  link = &(mtx->hdr.col[org]);
  if( rng == mtx_ALL_ROWS ) {
    while( NOTNULL(*link) )
      delete_from_col(link);
  } else {
    int32 *tocur = mtx->perm.row.org_to_cur, row;

    while( NOTNULL(*link) ) {
      row = (*link)->row;
      if( row == mtx_NONE || in_range(rng,tocur[row]) )
        delete_from_col(link);
      else
        link = &((*link)->next.row);
    }
  }
}

static void clean_row(mtx_matrix_t mtx, int32 org)
/**
 ***  Disclaims all elements in the given row which have already been
 ***  disclaimed once.
 ***  NOTE WELL: last_value_mtx must be set before this function is called.
 **/
{
  struct element_t **link;

  link = &(mtx->hdr.row[org]);
  while( NOTNULL(*link) )
    if( (*link)->row == mtx_NONE )
      delete_from_row(link);
    else
      link = &((*link)->next.col);
}

static void clean_col(mtx_matrix_t mtx, int32 org)
/**
 ***  Disclaims all elements in the given col which have already been
 ***  disclaimed once.
 ***  NOTE WELL: last_value_mtx must be set before this function is called.
 **/
{
  struct element_t **link;

  link = &(mtx->hdr.col[org]);
  while( NOTNULL(*link) )
    if( (*link)->col == mtx_NONE )
      delete_from_col(link);
    else
      link = &((*link)->next.row);
}
/********************************************************************/

mtx_coord_t *mtx_coord(mtx_coord_t *coordp, int32 row, int32 col)
{
   coordp->row = row;
   coordp->col = col;
   return(coordp);
}

mtx_range_t *mtx_range(mtx_range_t *rangep, int32 low, int32 high)
{
   rangep->low = low;
   rangep->high = high;
   return(rangep);
}

mtx_region_t *mtx_region( mtx_region_t *regionp, int32 rowlow, int32 rowhigh,
            int32 collow, int32 colhigh)
{
   mtx_range(&(regionp->row),rowlow,rowhigh);
   mtx_range(&(regionp->col),collow,colhigh);
   return(regionp);
}

static int g_mtx_debug_redirect = 0;

void mtx_debug_redirect_freeze() {
  g_mtx_debug_redirect = 1;
}

static void mtx_redirectErrors(FILE *f) {
  if (!g_mtx_debug_redirect) {
    assert(f != NULL);
    g_mtxerr = f;
  }
}

mtx_matrix_t mtx_create(void)
{
  mtx_matrix_t mtx;

  mtx_redirectErrors(stderr); /* call mtx_debug_redirect_freeze() to bypass */

  mtx = alloc_header();
  mtx->order = ZERO;
  mtx->capacity = ZERO;
  mtx->hdr.row = mtx->hdr.col = NULL;
  mtx->perm.row.org_to_cur = mtx_alloc_perm(ZERO);
  mtx->perm.row.cur_to_org = mtx_alloc_perm(ZERO);
  mtx->perm.row.parity = EVEN;
  mtx->perm.col.org_to_cur = mtx_alloc_perm(ZERO);
  mtx->perm.col.cur_to_org = mtx_alloc_perm(ZERO);
  mtx->perm.col.parity = EVEN;
  mtx->perm.transpose = mtx_NONE;
  mtx->data =
    (struct structural_data_t *)ascmalloc(sizeof(struct structural_data_t));
  mtx->data->symbolic_rank = -1;
  mtx->data->nblocks = -1;
  mtx->data->block = NULL;
  mtx->master = NULL;
  mtx->nslaves = ZERO;
  mtx->slaves = NULL;
  mtx->last_value = NULL;
  mtx->ms = mem_create_store(LENMAGIC, WIDTHMAGIC/sizeof(struct element_t) - 1,
              sizeof(struct element_t),10,2048);
  return(mtx);
}

/**
 ** Slave headers share the master's:
 ** perm.anything, data->anything, ms->anything, order and capacity.
 ** Unique to the slave are:
 ** hdr.anything
 ** last_value.
 ** All slaves of a master appear as pointers in the masters slaves list
 ** which is nslaves long. slaves do not have slaves.
 **/
mtx_matrix_t mtx_create_slave(mtx_matrix_t master)
{
  mtx_matrix_t mtx, *mlist;
  int32 newcnt;
  mtx_redirectErrors(stderr); /* call mtx_debug_redirect_freeze() to bypass */
  if(!mtx_check_matrix(master) || ISSLAVE(master)) return NULL;

  mtx = alloc_header();
  if (ISNULL(mtx)) return mtx;
  mtx->order = master->order;
  mtx->capacity = master->capacity;
  mtx->nslaves = ZERO;
  mtx->hdr.row = calloc_nz_header(master->capacity);
  if (ISNULL(mtx->hdr.row)) {
    ascfree(mtx);
    FPRINTF(g_mtxerr,"mtx_create_slave: Insufficient memory.\n");
    return NULL;
  }
  mtx->hdr.col = calloc_nz_header(master->capacity);
  if (ISNULL(mtx->hdr.col)) {
    ascfree(mtx->hdr.row);
    ascfree(mtx);
    FPRINTF(g_mtxerr,"mtx_create_slave: Insufficient memory.\n");
    return NULL;
  }
  mtx->perm.row.org_to_cur = master->perm.row.org_to_cur;
  mtx->perm.row.cur_to_org = master->perm.row.cur_to_org;
  mtx->perm.row.parity = EVEN; /* dont look at this again*/
  mtx->perm.col.org_to_cur = master->perm.col.org_to_cur;
  mtx->perm.col.cur_to_org = master->perm.col.cur_to_org;
  mtx->perm.col.parity = EVEN; /* dont look at this again*/
  mtx->perm.transpose = mtx_NONE; /* dont look at this again*/
  mtx->slaves = NULL;
  mtx->data = master->data;
  mtx->master = master;
  mtx->last_value = NULL;
  mtx->ms = master->ms;
  newcnt = master->nslaves + 1;
  if (newcnt == 1) {
    mlist = (mtx_matrix_t *)ascmalloc(sizeof(mtx_matrix_t));
  } else {
    mlist =
      (mtx_matrix_t *)ascrealloc(master->slaves,newcnt*sizeof(mtx_matrix_t));
  }
  if (ISNULL(mlist)) {
    ascfree(mtx->hdr.row);
    ascfree(mtx->hdr.col);
    ascfree(mtx);
    FPRINTF(g_mtxerr,"mtx_create_slave: Insufficient memory.\n");
    return NULL;
  } else {
    master->slaves = mlist;
    master->slaves[master->nslaves] = mtx;
    master->nslaves = newcnt;
  }
  return(mtx);
}

static void destroy_slave(mtx_matrix_t master, mtx_matrix_t mtx)
{
  int32 sid,i;

  if(!mtx_check_matrix(mtx) || !mtx_check_matrix(master)) {
    FPRINTF(g_mtxerr,
      "destroy_slave(mtx.c) called with corrupt slave or master\n");
    return;
  }
  if (ISSLAVE(master) || !ISSLAVE(mtx)) {
    FPRINTF(g_mtxerr,
      "destroy_slave(mtx.c) called with mismatched slave/master pair\n");
    return;
  }
  /* get slave index */
  sid = -1;
  for (i=0; i<master->nslaves; i++) {
    if (master->slaves[i] == mtx) {
      sid = i;
    }
  }
  if (sid < 0) {
    FPRINTF(g_mtxerr,
      "destroy_slave(mtx.c) called with mismatched slave/master pair\n");
    return;
  }
  /* move the rest of the slaves up the list */
  for (i = sid; i < master->nslaves -1;) {
    master->slaves[i] = master->slaves[++i];
  }
  (master->nslaves)--;
  /* we will not realloc smaller. */
  if (master->nslaves == 0 && NOTNULL(master->slaves)) {
     ascfree(master->slaves);
     master->slaves = NULL;
  }
  mtx_clear_region(mtx,mtx_ENTIRE_MATRIX);
  free_nz_header(mtx->hdr.row);
  free_nz_header(mtx->hdr.col);
  free_header(mtx);
}

void mtx_destroy(mtx_matrix_t mtx)
{
  int32 i;
  if(!mtx_check_matrix(mtx)) return;
  if (ISSLAVE(mtx)) {
    destroy_slave(mtx->master,mtx);
    return;
  }

  for (i=0; i<mtx->nslaves; i++) {
    if (mtx_check_matrix(mtx->slaves[i])) {
      free_nz_header(mtx->slaves[i]->hdr.row);
      free_nz_header(mtx->slaves[i]->hdr.col);
      free_header(mtx->slaves[i]);
      mtx->slaves[i] = NULL;
    } else {
      FPRINTF(g_mtxerr,
        "mtx_destroy: Corrupt slave found while destroying master.\n");
      FPRINTF(g_mtxerr,
        "             Slave %d being abandoned.\n",i);
    }
  }
  mtx->nslaves = 0;

  mtx_clear_region(mtx,mtx_ENTIRE_MATRIX);
  /* mtx_check_matrix is called again if MTX_DEBUG*/
  /**
   ***  The fishnet structure is
   ***  destroyed, just leaving
   ***  the headers and maybe ms.
   **/
  mem_destroy_store(mtx->ms);

  free_nz_header(mtx->hdr.row);
  free_nz_header(mtx->hdr.col);

  mtx_free_perm(mtx->perm.row.org_to_cur);
  mtx_free_perm(mtx->perm.row.cur_to_org);
  mtx_free_perm(mtx->perm.col.org_to_cur);
  mtx_free_perm(mtx->perm.col.cur_to_org);

  if( NOTNULL(mtx->data->block) ) {
    ascfree(mtx->data->block);
  }
  ascfree(mtx->data);

  free_header(mtx);
}

mtx_sparse_t *mtx_create_sparse(int32 cap) {
  mtx_sparse_t *ret;
  ret = (mtx_sparse_t *)ascmalloc(sizeof(mtx_sparse_t));
  if (ISNULL(ret)) {
    FPRINTF(g_mtxerr,"ERROR: (mtx_create_sparse) Insufficient memory.\n");
    return ret;
  }
  ret->data = (real64 *)ascmalloc(cap*sizeof(real64));
  if (ISNULL(ret->data)) {
    FPRINTF(g_mtxerr,"ERROR: (mtx_create_sparse) Insufficient memory.\n");
    ascfree(ret);
    return NULL;
  }
  ret->idata = (int32 *)ascmalloc(cap*sizeof(int32));
  if (ISNULL(ret->idata)) {
    FPRINTF(g_mtxerr,"ERROR: (mtx_create_sparse) Insufficient memory.\n");
    ascfree(ret->data);
    ascfree(ret);
    return NULL;
  }
  ret->cap = cap;
  ret->len = 0;
  return ret;
}

void mtx_destroy_sparse(mtx_sparse_t *ret)
{
  if (ISNULL(ret)) return;
  if (NOTNULL(ret->idata)) ascfree(ret->idata);
  if (NOTNULL(ret->data)) ascfree(ret->data);
  ascfree(ret);
}

void mtx_destroy_blocklist(mtx_block_t *bl)
{
  if (ISNULL(bl)) return;
  if (NOTNULL(bl->block) && bl->nblocks>0) ascfree(bl->block);
  ascfree(bl);
}

mtx_matrix_t mtx_duplicate_region(mtx_matrix_t mtx, mtx_region_t *reg,
                                  real64 drop)
{
  mtx_matrix_t new;
  int32 org_row;
  struct element_t *elt;

  if (!mtx_check_matrix(mtx) ||
      (mtx->master != NULL && !mtx_check_matrix(mtx->master))) {
    return NULL;
  }
  drop = fabs(drop);
  if (ISSLAVE(mtx)) {
    new = mtx_create_slave(mtx->master);
  } else {
    new = mtx_create_slave(mtx);
  }
  if (new == NULL) {
    FPRINTF(g_mtxerr,"ERROR: (mtx_duplicate_region) Insufficient memory.\n");
    return NULL;
  }
  if (reg==mtx_ENTIRE_MATRIX){
    if (drop >0.0) { /* copy with drop tolerance */
      for( org_row=0 ; org_row < mtx->order ; ++org_row ) {
        for( elt = mtx->hdr.row[org_row] ;
             NOTNULL(elt) ; elt = elt->next.col ) {
          if( elt->row != mtx_NONE  && fabs(elt->value) > drop) {
            mtx_create_element_value(new,elt->row,elt->col,elt->value);
          }
        }
      }
    } else {
      for( org_row=0 ; org_row < mtx->order ; ++org_row ) {
        for( elt = mtx->hdr.row[org_row] ;
             NOTNULL(elt) ; elt = elt->next.col ) {
          if( elt->row != mtx_NONE ) {
            mtx_create_element_value(new,elt->row,elt->col,elt->value);
          }
        }
      }
    }
  } else {
    int32 *rowcur, *colcur, row,col;
    mtx_range_t *colrng,*rowrng;
    rowcur=mtx->perm.row.org_to_cur;
    colcur=mtx->perm.col.org_to_cur;
    rowrng=&(reg->row);
    colrng=&(reg->col);
    if (drop >0.0) { /* copy with drop tolerance */
      for( org_row=0 ; org_row < mtx->order ; ++org_row ) {
        row=rowcur[org_row];
        if (in_range(rowrng,row)) {
          for( elt = mtx->hdr.row[org_row] ;
               NOTNULL(elt) ; elt = elt->next.col ) {
            col=colcur[elt->col];
            if( in_range(colrng,col) && elt->row != mtx_NONE &&
                fabs(elt->value) > drop ) {
                                        /* don't copy^bad^elements */
              mtx_create_element_value(new,elt->row,elt->col,elt->value);
            }
          }
        }
      }
    } else {
      for( org_row=0 ; org_row < mtx->order ; ++org_row ) {
        row=rowcur[org_row];
        if (in_range(rowrng,row)) {
          for( elt = mtx->hdr.row[org_row] ;
               NOTNULL(elt) ; elt = elt->next.col ) {
            col=colcur[elt->col];
            if( in_range(colrng,col) && elt->row != mtx_NONE ) {
                                        /* don't copy^bad^elements */
              mtx_create_element_value(new,elt->row,elt->col,elt->value);
            }
          }
        }
      }
    }
  }

  return new;
}

/* WARNING: this function assumes sufficient memory is available.
 * it needs to be check for null returns in allocs and doesn't.
 */
mtx_matrix_t mtx_copy_options(mtx_matrix_t mtx, boolean blocks,
                              boolean incidence, mtx_region_t *reg,
                              real64 drop)
{
  mtx_matrix_t copy;
  int32 org_row;
  struct element_t *elt;
  int32 bnum;

  if (!mtx_check_matrix(mtx)) return NULL;
  drop = fabs(drop);

  /* create same size matrix */
  copy = alloc_header();
  copy->order = mtx->order;
  copy->capacity = mtx->capacity;
  copy->master = NULL;
  copy->nslaves = ZERO;
  copy->slaves = NULL;
  copy->last_value = NULL;

  /* copy headers and fishnet */
  copy->hdr.row = calloc_nz_header(copy->capacity);
  copy->hdr.col = calloc_nz_header(copy->capacity);
  if (incidence) {
    struct mem_statistics stat;
    int s_expool;
    mem_get_stats(&stat,mtx->ms);
    s_expool = MAX(2048,stat.elt_total/stat.str_wid);
    copy->ms = mem_create_store(LENMAGIC,
                 WIDTHMAGIC/sizeof(struct element_t) - 1,
                 sizeof(struct element_t),
                 10,s_expool-LENMAGIC);
    /* copy of a slave or master matrix will end up with an
     * initial mem_store that is the size cumulative of the
     * master and all its slaves since they share the ms.
     */
    if (reg==mtx_ENTIRE_MATRIX){
      if (drop >0.0) { /* copy with drop tolerance */
        for( org_row=0 ; org_row < mtx->order ; ++org_row ) {
          for( elt = mtx->hdr.row[org_row] ;
               NOTNULL(elt) ; elt = elt->next.col ) {
            if( elt->row != mtx_NONE  && fabs(elt->value) > drop) {
              mtx_create_element_value(copy,elt->row,elt->col,elt->value);
            }
          }
        }
      } else {
        for( org_row=0 ; org_row < mtx->order ; ++org_row ) {
          for( elt = mtx->hdr.row[org_row] ;
               NOTNULL(elt) ; elt = elt->next.col ) {
            if( elt->row != mtx_NONE ) {
              mtx_create_element_value(copy,elt->row,elt->col,elt->value);
            }
          }
        }
      }
    } else {
      int32 *rowcur, *colcur, row,col;
      mtx_range_t *colrng,*rowrng;
      rowcur=mtx->perm.row.org_to_cur;
      colcur=mtx->perm.col.org_to_cur;
      rowrng=&(reg->row);
      colrng=&(reg->col);
      if (drop >0.0) { /* copy with drop tolerance */
        for( org_row=0 ; org_row < mtx->order ; ++org_row ) {
          row=rowcur[org_row];
          if (in_range(rowrng,row)) {
            for( elt = mtx->hdr.row[org_row] ;
                 NOTNULL(elt) ; elt = elt->next.col ) {
              col=colcur[elt->col];
              if( in_range(colrng,col) && elt->row != mtx_NONE &&
                  fabs(elt->value) > drop ) {
                                          /* don't copy^bad^elements */
                mtx_create_element_value(copy,elt->row,elt->col,elt->value);
              }
            }
          }
        }
      } else {
        for( org_row=0 ; org_row < mtx->order ; ++org_row ) {
          row=rowcur[org_row];
          if (in_range(rowrng,row)) {
            for( elt = mtx->hdr.row[org_row] ;
                 NOTNULL(elt) ; elt = elt->next.col ) {
              col=colcur[elt->col];
              if( in_range(colrng,col) && elt->row != mtx_NONE ) {
                                          /* don't copy^bad^elements */
                mtx_create_element_value(copy,elt->row,elt->col,elt->value);
              }
            }
          }
        }
      }
    }
  } else {
    copy->ms = mem_create_store(LENMAGIC,
                 WIDTHMAGIC/sizeof(struct element_t) - 1,
                 sizeof(struct element_t),10,2048);
    /* copy of a slave or master matrix will end up with an
       initial mem_store that is the default size */
  }

  /* copy permutation information */
  copy->perm.row.org_to_cur = mtx_alloc_perm(copy->capacity);
  mtx_copy_perm(copy->perm.row.org_to_cur,
            mtx->perm.row.org_to_cur,copy->capacity);
  copy->perm.row.cur_to_org = mtx_alloc_perm(copy->capacity);
  mtx_copy_perm(copy->perm.row.cur_to_org,
            mtx->perm.row.cur_to_org,copy->capacity);
  copy->perm.row.parity = mtx_row_parity(mtx);

  copy->perm.col.org_to_cur = mtx_alloc_perm(copy->capacity);
  mtx_copy_perm(copy->perm.col.org_to_cur,
            mtx->perm.col.org_to_cur,copy->capacity);
  copy->perm.col.cur_to_org = mtx_alloc_perm(copy->capacity);
  mtx_copy_perm(copy->perm.col.cur_to_org,
            mtx->perm.col.cur_to_org,copy->capacity);
  copy->perm.col.parity = mtx_col_parity(mtx);

  /* copy (or not) block data */
  copy->data =
    (struct structural_data_t *)ascmalloc(sizeof(struct structural_data_t));
  if (blocks) {
    copy->data->symbolic_rank = mtx->data->symbolic_rank;
    copy->data->nblocks = mtx->data->nblocks;
    copy->data->block = mtx->data->nblocks > 0 ? (mtx_region_t *)
      ascmalloc( mtx->data->nblocks*sizeof(mtx_region_t) ) : NULL;
    for( bnum=0; bnum < mtx->data->nblocks; bnum++ ) {
      copy->data->block[bnum].row.low = mtx->data->block[bnum].row.low;
      copy->data->block[bnum].row.high = mtx->data->block[bnum].row.high;
      copy->data->block[bnum].col.low = mtx->data->block[bnum].col.low;
      copy->data->block[bnum].col.high = mtx->data->block[bnum].col.high;
    }
  } else {
    copy->data->symbolic_rank=0;
    copy->data->nblocks=0;
    copy->data->block =NULL;
  }

  return(copy);
}

int32 mtx_order( mtx_matrix_t mtx)
{
   if (!mtx_check_matrix(mtx)) return -1;
   return(mtx->order);
}

int32 mtx_capacity( mtx_matrix_t mtx)
{
   if (!mtx_check_matrix(mtx)) return -1;
   return(mtx->capacity);
}

/* this is only to be called from in mtx_set_order */
static void trim_incidence(mtx_matrix_t mtx, int32 order)
{
   int32 ndx;
   int32 *toorg;

   last_value_matrix = mtx;
   toorg = mtx->perm.col.cur_to_org;
   for( ndx = order ; ndx < mtx->order ; ++ndx )
     blast_col(mtx,toorg[ndx],mtx_ALL_ROWS);
   for( ndx = ZERO ; ndx < order ; ++ndx )
     clean_row(mtx,toorg[ndx]);

   toorg = mtx->perm.row.cur_to_org;
   for( ndx = order ; ndx < mtx->order ; ++ndx )
     blast_row(mtx,toorg[ndx],mtx_ALL_COLS);
   for( ndx = ZERO ; ndx < order ; ++ndx )
     clean_col(mtx,toorg[ndx]);
}

/* this is only to be called from in mtx_set_order */
static void enlarge_nzheaders(mtx_matrix_t mtx, int32 order)
{
  struct element_t **newhdr;

  /* realloc does not initialize, so calloc is the best we can do here */

  newhdr = calloc_nz_header(order);
  copy_nz_header(newhdr,mtx->hdr.row,mtx->capacity);
  free_nz_header(mtx->hdr.row);
  mtx->hdr.row = newhdr;

  newhdr = calloc_nz_header(order);
  copy_nz_header(newhdr,mtx->hdr.col,mtx->capacity);
  free_nz_header(mtx->hdr.col);
  mtx->hdr.col = newhdr;
}

void mtx_set_order( mtx_matrix_t mtx, int32 order)
/**
 ***  This function will preserve the fact that all of the arrays are
 ***  "correct" out to the capacity of the arrays, not just out to the order
 ***  of the matrix.  In other words, extensions to orders which still do
 ***  not exceed the capacity of the arrays are trivial.
 **/
{
  int32 i;
  if(!mtx_check_matrix(mtx)) return;
  if (ISSLAVE(mtx)) {
    mtx_set_order(mtx->master,order);
    return;
  }
  /* we now have a master matrix */
  if( order < mtx->order ) {   /* Truncate */
    trim_incidence(mtx,order); /* clean master */
    for (i = 0; i < mtx->nslaves; i++) {
      trim_incidence(mtx->slaves[i],order);
    }
  }
  if (mtx->perm.transpose == mtx_NONE) {
    mtx->perm.transpose = 0;
  }

  if( order > mtx->capacity ) {
    int32 *newperm;
    int32 ndx;

    enlarge_nzheaders(mtx,order);
    for (i = 0; i < mtx->nslaves; i++) {
      enlarge_nzheaders(mtx->slaves[i],order);
    }

/* realloc not in order here.  Happens only on the master. */
    newperm = mtx_alloc_perm(order);
    mtx_copy_perm(newperm,mtx->perm.row.org_to_cur,mtx->capacity);
    for( ndx=mtx->capacity ; ndx < order ; ++ndx )
       newperm[ndx] = ndx;
    mtx_free_perm(mtx->perm.row.org_to_cur);
    mtx->perm.row.org_to_cur = newperm;

    newperm = mtx_alloc_perm(order);
    mtx_copy_perm(newperm,mtx->perm.row.cur_to_org,mtx->capacity);
    for( ndx=mtx->capacity ; ndx < order ; ++ndx )
       newperm[ndx] = ndx;
    mtx_free_perm(mtx->perm.row.cur_to_org);
    mtx->perm.row.cur_to_org = newperm;

    newperm = mtx_alloc_perm(order);
    mtx_copy_perm(newperm,mtx->perm.col.org_to_cur,mtx->capacity);
    for( ndx=mtx->capacity ; ndx < order ; ++ndx )
       newperm[ndx] = ndx;
    mtx_free_perm(mtx->perm.col.org_to_cur);
    mtx->perm.col.org_to_cur = newperm;

    newperm = mtx_alloc_perm(order);
    mtx_copy_perm(newperm,mtx->perm.col.cur_to_org,mtx->capacity);
    for( ndx=mtx->capacity ; ndx < order ; ++ndx )
       newperm[ndx] = ndx;
    mtx_free_perm(mtx->perm.col.cur_to_org);
    mtx->perm.col.cur_to_org = newperm;

    mtx->capacity = order;
    for (i = 0; i < mtx->nslaves; i++) {
      mtx->slaves[i]->capacity = order;
      /* point slaves at master perm again in case anything has changed */
      mtx->slaves[i]->perm.row.org_to_cur = mtx->perm.row.org_to_cur;
      mtx->slaves[i]->perm.row.cur_to_org = mtx->perm.row.cur_to_org;
      mtx->slaves[i]->perm.col.org_to_cur = mtx->perm.col.org_to_cur;
      mtx->slaves[i]->perm.col.cur_to_org = mtx->perm.col.cur_to_org;
    }
  }
  mtx->order = order;
  for (i = 0; i < mtx->nslaves; i++) {
    mtx->slaves[i]->order = order;
  }
}

void mtx_clear_coord(mtx_matrix_t mtx, int32 row, int32 col)
{
   struct element_t **rlink, **clink;
   int32 org_row, org_col;

#if MTX_DEBUG
   if( !mtx_check_matrix(mtx) ) return; /*ben*/
#endif

   last_value_matrix = mtx;
   org_row = mtx->perm.row.cur_to_org[row];
   org_col = mtx->perm.col.cur_to_org[col];
   rlink = &(mtx->hdr.row[org_row]);
   for( ; NOTNULL(*rlink) && (*rlink)->col != org_col ; )
      rlink = &((*rlink)->next.col);
   if( ISNULL(*rlink) ) return;
   clink = &(mtx->hdr.col[org_col]);
   for( ; NOTNULL(*clink)  && (*clink)->row != org_row ; )
      clink = &((*clink)->next.row);
   delete_from_row(rlink);
   delete_from_col(clink);
}

void mtx_clear_row( mtx_matrix_t mtx, int32 row, mtx_range_t *rng)
{
   struct element_t **rlink, **clink;
   int32 org_row;

#if MTX_DEBUG
   if( !mtx_check_matrix(mtx) ) return; /*ben*/
#endif

   last_value_matrix = mtx;
   org_row = mtx->perm.row.cur_to_org[row];
   rlink = &(mtx->hdr.row[org_row]);
   if( rng == mtx_ALL_COLS ) {
      while( NOTNULL(*rlink)  ) {
         clink = &(mtx->hdr.col[(*rlink)->col]);
         for( ; NOTNULL(*clink) && (*clink)->row != org_row ; )
            clink = &((*clink)->next.row);
         delete_from_row(rlink);
         delete_from_col(clink);
      }
   } else if( rng->high >= rng->low ) {
      while( NOTNULL(*rlink) ) {
         if( in_range(rng,mtx->perm.col.org_to_cur[(*rlink)->col]) ) {
            clink = &(mtx->hdr.col[(*rlink)->col]);
            for( ; NOTNULL(*clink) && (*clink)->row != org_row ; )
               clink = &((*clink)->next.row);
            delete_from_row(rlink);
            delete_from_col(clink);
         } else
            rlink = &((*rlink)->next.col);
      }
   }
}

void mtx_clear_col( mtx_matrix_t mtx, int32 col, mtx_range_t *rng)
{
   struct element_t **clink, **rlink;
   int32 org_col;

#if MTX_DEBUG
   if( !mtx_check_matrix(mtx) ) return;
#endif

   last_value_matrix = mtx;
   org_col = mtx->perm.col.cur_to_org[col];
   clink = &(mtx->hdr.col[org_col]);
   if( rng == mtx_ALL_ROWS ) {
      while( NOTNULL(*clink) ) {
         rlink = &(mtx->hdr.row[(*clink)->row]);
         for( ; NOTNULL(*rlink) && (*rlink)->col != org_col ; )
            rlink = &((*rlink)->next.col);
         delete_from_col(clink);
         delete_from_row(rlink);
      }
   } else if( rng->high >= rng->low ) {
      while( NOTNULL(*clink) ) {
         if( in_range(rng,mtx->perm.row.org_to_cur[(*clink)->row]) ) {
            rlink = &(mtx->hdr.row[(*clink)->row]);
            for( ; NOTNULL(*rlink) && (*rlink)->col != org_col ; )
               rlink = &((*rlink)->next.col);
            delete_from_col(clink);
            delete_from_row(rlink);
         } else
            clink = &((*clink)->next.row);
      }
   }
}

void mtx_clear_region(mtx_matrix_t mtx, mtx_region_t *region)
{
#if MTX_DEBUG
  if(!mtx_check_matrix(mtx)) return;
#endif

  if( region == mtx_ENTIRE_MATRIX && !ISSLAVE(mtx)) {
    /* damn the torpedos, wipe that sucker and slaves fast */
    nuke_fishnet(mtx);
  } else {
    int32 ndx;
    int32 *toorg;
    mtx_region_t reg;

    if (region == mtx_ENTIRE_MATRIX) {
      reg.row.low = reg.col.low = 0;
      reg.row.high = reg.col.high = mtx->order-1;
    } else {
      reg = *region;
    }
    last_value_matrix = mtx;
    toorg = mtx->perm.row.cur_to_org;
    for( ndx = reg.row.low ; ndx <= reg.row.high ; ++ndx )
      blast_row(mtx,toorg[ndx],&(reg.col));
    toorg = mtx->perm.col.cur_to_org;
    for( ndx = reg.col.low ; ndx <= reg.col.high ; ++ndx )
      clean_col(mtx,toorg[ndx]);
  }
}


void mtx_reset_perm(mtx_matrix_t mtx)
{
  int32 ndx;

#if MTX_DEBUG
  if(!mtx_check_matrix(mtx)) return;
#endif
  if (ISSLAVE(mtx)) {
    mtx_reset_perm(mtx->master);
    return;
  }

  for( ndx=ZERO ; ndx < mtx->capacity ; ++ndx ) {
    mtx->perm.row.org_to_cur[ndx] = ndx;
    mtx->perm.row.cur_to_org[ndx] = ndx;
    mtx->perm.col.org_to_cur[ndx] = ndx;
    mtx->perm.col.cur_to_org[ndx] = ndx;
  }
  mtx->perm.row.parity = mtx->perm.col.parity = EVEN;

  mtx->data->symbolic_rank = -1;
  mtx->data->nblocks = -1;
  if( NOTNULL(mtx->data->block) ) {
    ascfree(mtx->data->block);
    mtx->data->block = NULL;
  }
}

void mtx_clear(mtx_matrix_t mtx)
{
  if (ISSLAVE(mtx)) {
    mtx_clear_region(mtx->master,mtx_ENTIRE_MATRIX);
  }
  mtx_clear_region(mtx,mtx_ENTIRE_MATRIX);
  mtx_reset_perm(mtx);
}

real64 mtx_value(mtx_matrix_t mtx, mtx_coord_t *coord)
{
   struct element_t *elt;
   register int32 org_row,org_col;

#if MTX_DEBUG
   if(!mtx_check_matrix(mtx)) return D_ZERO;
#endif
   org_row = mtx->perm.row.cur_to_org[coord->row];
   org_col = mtx->perm.col.cur_to_org[coord->col];
   elt = mtx_find_element(mtx,org_row,org_col);
   mtx->last_value = elt;
   return( ISNULL(elt) ? D_ZERO : elt->value );
}

void mtx_set_value( mtx_matrix_t mtx, mtx_coord_t *coord, real64 value)
{
  register int32 org_row,org_col;

#if MTX_DEBUG
  if(!mtx_check_matrix(mtx)) return;
#endif
  org_row = mtx->perm.row.cur_to_org[coord->row];
  org_col = mtx->perm.col.cur_to_org[coord->col];
  if ( NOTNULL(mtx->last_value) &&
      mtx->last_value->row==org_row &&
      mtx->last_value->col==org_col ) {
      mtx->last_value->value = value;
  } else {
    struct element_t *elt;
    if( ISNULL(elt = mtx_find_element(mtx,org_row,org_col)) ) {
      if (value != D_ZERO ) {
        elt = mtx_create_element_value(mtx,org_row,org_col,value);
      }
    } else {
      elt->value = value;
    }
  }
}

void mtx_fill_value(mtx_matrix_t mtx, mtx_coord_t *coord, real64 value)
{
   register int32 org_row,org_col;

#if MTX_DEBUG
   if(!mtx_check_matrix(mtx)) return;
#endif
   org_row = mtx->perm.row.cur_to_org[coord->row];
   org_col = mtx->perm.col.cur_to_org[coord->col];
   mtx_create_element_value(mtx,org_row,org_col,value);
}

void mtx_fill_org_value(mtx_matrix_t mtx, const mtx_coord_t *coord,
            real64 value)
{
#if MTX_DEBUG
   if(!mtx_check_matrix(mtx)) return;
#endif
   mtx_create_element_value(mtx,coord->row,coord->col,value);
}

/* sparse matrix assembly of potentially duplicate fills */
/* takes a matrix, assumed to have redundant and otherwise insane incidences
  created by 'misusing mtx_fill_value' and sums all like entries, eliminating
  the duplicates and the zeroes. Returns -# of elements removed, or 1 if bad.
  returns 1 if fails for some reason, 0 otherwise.
  Could stand to have the error messages it emits improved.
  Could stand to take a rowrange or a rowlist, a colrange or a collist,droptol.
  algorithm: O(3tau)
    establish a vector mv of columns needing cleaning markers
    set lowmark=-1, highmark= -2
    for rows
      if row empty continue
      establish a vector ev of elt pointers set null
      prevlink = &header 1st elt pointer
      for each elt in row
        if elt is not marked in ev
          mark ev with elt, update prevlink
        else
          add value to elt pointed at by ev
          mark col in mv as needing cleaning, updating lowmark,highmark
          delete elt
        endif
      endfor
      for each elt in row
        unmark ev with elt
        if elt is 0, delete and mark col in mv.
      endfor
    endfor
    for i = lowmark to highmark
      clean col
    endfor
   My god, the tricks we play on a linked list.
*/
int32 mtx_assemble(mtx_matrix_t mtx)
{
  char *real_mv=NULL, *mv;
  struct element_t **real_ev=NULL, **ev, **prevlink, *elt;
  int32 orgrow, orgcol, lowmark,highmark,dup;
  /* note elt and prevlink could be combined, but the code is
     unreadable utterly if you do so. */

  if (ISNULL(mtx)) return 1;
  if (mtx->order <1) return 0;
  real_mv = mtx_null_mark(mtx->order+1);
  mv = real_mv+1;
  if (ISNULL(real_mv))  return 1;
  real_ev = mtx_null_vector(mtx->order+1);
  ev = real_ev+1;
  if (ISNULL(real_ev)) {
    mtx_null_mark_release();
    return 1;
  }
  /* we have allocated arrays which include a -1 element to buy
     ourselves an awful convenient lot of safety. */
  lowmark=-1;
  highmark=-2;
  dup = 0;
  last_value_matrix = mtx;
  for (orgrow=0; orgrow < mtx->order; orgrow++) {
    elt = mtx->hdr.row[orgrow];
    prevlink = &(mtx->hdr.row[orgrow]);
    while (NOTNULL(elt)) {
      if (ISNULL(ev[elt->col])) {
        ev[elt->col] = elt;         /* mark first elt found for this col */
        prevlink= &(elt->next.col); /* collect pointer to where we go next */
        elt = elt->next.col;        /* go there */
      } else {
        /* elt is duplicate and must die */
        dup++;
        ev[elt->col]->value += elt->value; /* accumulate value */
        /* update lowmark, highmark . this is a debatable implementation. */
        /* for small mods on large matrices, this is sane */
        if (lowmark > -1) { /* not first mark */
          if (elt->col < lowmark && elt->col > -1) {
            lowmark = elt->col;
          }
          if (elt->col > highmark && elt->col < mtx->order) {
            highmark = elt->col;
          }
        } else { /* very first mark */
          if (elt->col > -1 && elt->col < mtx->order) {
            lowmark = highmark = elt->col;
          }
        }
        mv[elt->col] = 1; /* mark column as to be cleaned */
        delete_from_row(prevlink);
        elt = *prevlink;
      }
    }
    elt = mtx->hdr.row[orgrow];
    prevlink = &(mtx->hdr.row[orgrow]);
    while (NOTNULL(elt)) {  /* renull the accumulator and trash 0s */
      ev[elt->col] = NULL; /* regardless, reset accum */
      if (elt->value != D_ZERO) {
        prevlink= &(elt->next.col); /* collect pointer to where we go next */
        elt = elt->next.col;        /* go there */
      } else {
        /* this is still a debatable implementation. */
        if (lowmark > -1) { /* not first mark */
          if (elt->col < lowmark && elt->col > -1) {
            lowmark = elt->col;
          }
          if (elt->col > highmark && elt->col < mtx->order) {
            highmark = elt->col;
          }
        } else { /* very first mark */
          if (elt->col > -1 && elt->col < mtx->order) {
            lowmark = highmark = elt->col;
          }
        }
        mv[elt->col] = 1; /* mark column as to be cleaned */
        delete_from_row(prevlink);
        elt = *prevlink;
      }
    }
  }
  for (orgcol = lowmark; orgcol <= highmark; orgcol++) {
    if (mv[orgcol]) {
      clean_col(mtx,orgcol); /* scrap dups and 0s */
    }
  }
  mtx_null_mark_release();
  mtx_null_vector_release();
  return -dup;
}

void mtx_del_zr_in_row(mtx_matrix_t mtx, int32 row)
{
   register int32 org;
   struct element_t **rlink, **clink;

#if MTX_DEBUG
   if(!mtx_check_matrix(mtx)) return;
#endif
   org = mtx->perm.row.cur_to_org[row];
   rlink = &(mtx->hdr.row[org]);

   last_value_matrix = mtx;
   while( NOTNULL(*rlink) )
      if( (*rlink)->value == D_ZERO ) {
         clink = &(mtx->hdr.col[(*rlink)->col]);
         for( ; NOTNULL(*clink) && (*clink)->row != org ; )
            clink = &((*clink)->next.row);
         delete_from_row(rlink);
         delete_from_col(clink);
      } else
         rlink = &((*rlink)->next.col);
}

void mtx_del_zr_in_col(mtx_matrix_t mtx, int32 col)
{
   register int32 org;
   struct element_t **clink, **rlink;

#if MTX_DEBUG
   if(!mtx_check_matrix(mtx)) return;
#endif
   org = mtx->perm.col.cur_to_org[col];
   clink = &(mtx->hdr.col[org]);

   last_value_matrix = mtx;
   while( NOTNULL(*clink) )
      if( (*clink)->value == D_ZERO ) {
         rlink = &(mtx->hdr.row[(*clink)->row]);
         for( ; NOTNULL(*rlink) && (*rlink)->col != org ; )
            rlink = &((*rlink)->next.col);
         delete_from_col(clink);
         delete_from_row(rlink);
      } else
         clink = &((*clink)->next.row);
}

void mtx_del_zr_in_rowrange(mtx_matrix_t mtx, mtx_range_t *rng)
{
  register int32 org,row,rowhi, *toorg;
  struct element_t **rlink, **clink;

#if MTX_DEBUG
  if(!mtx_check_matrix(mtx)) return;
#endif
  rowhi=rng->high;
  toorg= mtx->perm.row.cur_to_org;
  last_value_matrix = mtx;
  for (row=rng->low; row <=rowhi; row++) {
    org = toorg[row];
    rlink = &(mtx->hdr.row[org]);

    while( NOTNULL(*rlink) ) {
      if( (*rlink)->value == D_ZERO ) {
        clink = &(mtx->hdr.col[(*rlink)->col]);
        for( ; NOTNULL(*clink) && (*clink)->row != org ; )
          clink = &((*clink)->next.row);
        delete_from_row(rlink);
        delete_from_col(clink);
      } else {
        rlink = &((*rlink)->next.col);
      }
    }

  }
}

void mtx_del_zr_in_colrange(mtx_matrix_t mtx, mtx_range_t *rng)
{
  register int32 org,col,colhi, *toorg;
  struct element_t **clink, **rlink;

#if MTX_DEBUG
  if(!mtx_check_matrix(mtx)) return;
#endif
  colhi=rng->high;
  toorg= mtx->perm.col.cur_to_org;
  last_value_matrix = mtx;
  for (col=rng->low; col <=colhi; col++) {
    org = toorg[col];
    clink = &(mtx->hdr.col[org]);

    while( NOTNULL(*clink) ) {
      if( (*clink)->value == D_ZERO ) {
        rlink = &(mtx->hdr.row[(*clink)->row]);
        for( ; NOTNULL(*rlink) && (*rlink)->col != org ; )
          rlink = &((*rlink)->next.col);
        delete_from_col(clink);
        delete_from_row(rlink);
      } else {
        clink = &((*clink)->next.row);
      }
    }

  }
}

void mtx_steal_org_row_vec(mtx_matrix_t mtx, int32 row,
                           real64 *vec, mtx_range_t *rng)
{
  struct element_t **rlink, **clink;
  int32 org_row;

#if MTX_DEBUG
  if( !mtx_check_matrix(mtx) ) return;
#endif

  last_value_matrix = mtx;
  org_row = mtx->perm.row.cur_to_org[row];
  rlink = &(mtx->hdr.row[org_row]);
  if( rng == mtx_ALL_COLS ) {
    while( NOTNULL(*rlink)  ) {
      vec[(*rlink)->col] = (*rlink)->value;
      clink = &(mtx->hdr.col[(*rlink)->col]);
      for( ; NOTNULL(*clink) && (*clink)->row != org_row ; )
        clink = &((*clink)->next.row);
      delete_from_row(rlink);
      delete_from_col(clink);
    }
  } else if( rng->high >= rng->low ) {
    int32 *tocur;
    tocur = mtx->perm.col.org_to_cur;
    while( NOTNULL(*rlink) ) {
      if( in_range(rng,tocur[(*rlink)->col]) ) {
        vec[(*rlink)->col] = (*rlink)->value;
        clink = &(mtx->hdr.col[(*rlink)->col]);
        for( ; NOTNULL(*clink) && (*clink)->row != org_row ; )
          clink = &((*clink)->next.row);
        delete_from_row(rlink);
        delete_from_col(clink);
      } else {
         rlink = &((*rlink)->next.col);
      }
    }
  }
}

void mtx_steal_org_col_vec(mtx_matrix_t mtx, int32 col,
                           real64 *vec, mtx_range_t *rng)
{
  struct element_t **clink, **rlink;
  int32 org_col;

#if MTX_DEBUG
  if( !mtx_check_matrix(mtx) ) return;
#endif

  last_value_matrix = mtx;
  org_col = mtx->perm.col.cur_to_org[col];
  clink = &(mtx->hdr.col[org_col]);
  if( rng == mtx_ALL_ROWS ) {
    while( NOTNULL(*clink) ) {
      vec[(*clink)->row] = (*clink)->value;
      rlink = &(mtx->hdr.row[(*clink)->row]);
      for( ; NOTNULL(*rlink) && (*rlink)->col != org_col ; )
        rlink = &((*rlink)->next.col);
      delete_from_col(clink);
      delete_from_row(rlink);
    }
  } else if( rng->high >= rng->low ) {
    int32 *tocur;
    tocur = mtx->perm.row.org_to_cur;
    while( NOTNULL(*clink) ) {
      if( in_range(rng,tocur[(*clink)->row]) ) {
        vec[(*clink)->row] = (*clink)->value;
        rlink = &(mtx->hdr.row[(*clink)->row]);
        for( ; NOTNULL(*rlink) && (*rlink)->col != org_col ; )
          rlink = &((*rlink)->next.col);
        delete_from_col(clink);
        delete_from_row(rlink);
      } else {
        clink = &((*clink)->next.row);
      }
    }
  }
}

void mtx_steal_cur_row_vec(mtx_matrix_t mtx, int32 row,
                           real64 *vec, mtx_range_t *rng)
{
  struct element_t **rlink, **clink;
  int32 org_row, *tocur;

#if MTX_DEBUG
  if( !mtx_check_matrix(mtx) ) return;
#endif

  tocur = mtx->perm.col.org_to_cur;
  last_value_matrix = mtx;
  org_row = mtx->perm.row.cur_to_org[row];
  rlink = &(mtx->hdr.row[org_row]);
  if( rng == mtx_ALL_COLS ) {
    while( NOTNULL(*rlink)  ) {
      vec[tocur[(*rlink)->col]] = (*rlink)->value;
      clink = &(mtx->hdr.col[(*rlink)->col]);
      for( ; NOTNULL(*clink) && (*clink)->row != org_row ; )
        clink = &((*clink)->next.row);
      delete_from_row(rlink);
      delete_from_col(clink);
    }
  } else if( rng->high >= rng->low ) {
    while( NOTNULL(*rlink) ) {
      if( in_range(rng,tocur[(*rlink)->col]) ) {
        vec[tocur[(*rlink)->col]] = (*rlink)->value;
        clink = &(mtx->hdr.col[(*rlink)->col]);
        for( ; NOTNULL(*clink) && (*clink)->row != org_row ; )
          clink = &((*clink)->next.row);
        delete_from_row(rlink);
        delete_from_col(clink);
      } else {
         rlink = &((*rlink)->next.col);
      }
    }
  }
}

void mtx_steal_cur_col_vec(mtx_matrix_t mtx, int32 col,
                           real64 *vec, mtx_range_t *rng)
{
  struct element_t **clink, **rlink;
  int32 org_col, *tocur;

#if MTX_DEBUG
  if( !mtx_check_matrix(mtx) ) return;
#endif

  tocur = mtx->perm.row.org_to_cur;
  last_value_matrix = mtx;
  org_col = mtx->perm.col.cur_to_org[col];
  clink = &(mtx->hdr.col[org_col]);
  if( rng == mtx_ALL_ROWS ) {
    while( NOTNULL(*clink) ) {
      vec[tocur[(*clink)->row]] = (*clink)->value;
      rlink = &(mtx->hdr.row[(*clink)->row]);
      for( ; NOTNULL(*rlink) && (*rlink)->col != org_col ; )
        rlink = &((*rlink)->next.col);
      delete_from_col(clink);
      delete_from_row(rlink);
    }
  } else if( rng->high >= rng->low ) {
    while( NOTNULL(*clink) ) {
      if( in_range(rng,tocur[(*clink)->row]) ) {
        vec[tocur[(*clink)->row]] = (*clink)->value;
        rlink = &(mtx->hdr.row[(*clink)->row]);
        for( ; NOTNULL(*rlink) && (*rlink)->col != org_col ; )
          rlink = &((*rlink)->next.col);
        delete_from_col(clink);
        delete_from_row(rlink);
      } else {
        clink = &((*clink)->next.row);
      }
    }
  }
}

#ifdef DELETE_THIS_UNUSED_FUNCTION
/* Little function to enlarge the capacity of a sparse to the len given.
 * New memory allocated, if any, is not initialized in any way.
 * The pointer given is not itself changed, only its data.
 * This function in no way shrinks the data in a sparse.
 * (exception: buggy realloc implementations may shrink it, ack!)
 * Data/idata values up to the old capacity (ret->cap) are preserved.
 */
static int enlarge_sparse(mtx_sparse_t *ret, int32 len)
{
  int32 *inew;
  real64 *dnew;
  if (ISNULL(ret)) {
    FPRINTF(g_mtxerr,"ERROR: (mtx.c)    enlarge_sparse called with NULL.\n");
    return 1;
  }
  if (len <= ret->cap || len <1) return 0; /* already big enough */

  if (ret->idata == NULL) {
    inew = (int32 *)ascmalloc(sizeof(int32)*len);
  } else {
    inew = (int32 *)ascrealloc(ret->idata,sizeof(int32)*len);
  }
  if (ISNULL(inew)) {
    FPRINTF(g_mtxerr,"ERROR: (mtx.c)    enlarge_sparse\n");
    FPRINTF(g_mtxerr,"                  Insufficient memory.\n");
    return 1;
  }
  ret->idata = inew; /* dnew can still fail without losing inew memory. */

  if (ret->data == NULL) {
    dnew = (real64 *)ascmalloc(sizeof(real64)*len);
  } else {
    dnew = (real64 *)ascrealloc(ret->data,sizeof(real64)*len);
  }
  if (ISNULL(dnew)) {
    FPRINTF(g_mtxerr,"ERROR: (mtx.c)    enlarge_sparse\n");
    FPRINTF(g_mtxerr,"                  Insufficient memory.\n");
    return 1;
  }
  ret->data = dnew; /* we succeeded */
  ret->cap = len;
  return 0;
}
#endif /* DELETE_THIS_UNUSED_FUNCTION */


/* going to try to make steal also handle sparse creation ...*/
/* don't you dare, whoever you are! */
boolean mtx_steal_org_row_sparse(mtx_matrix_t mtx, int32 row,
                                 mtx_sparse_t *sp, mtx_range_t *rng)
{
  struct element_t **rlink, **clink;
  int32 org_row;
  int32 len,k;

#if MTX_DEBUG
  if( !mtx_check_matrix(mtx) ) return TRUE;
#endif
  if (sp == mtx_CREATE_SPARSE) {
    FPRINTF(g_mtxerr,"ERROR: (mtx.c)  mtx_steal_org_row_sparse called with\n");
    FPRINTF(g_mtxerr,"                mtx_CREATE_SPARSE. Not supported.\n");
    return TRUE;
  }
  if (rng == mtx_ALL_COLS) {
    len = mtx->order;
  } else {
    len = rng->high - rng->low +1;
  }
  if (sp->cap < len) {
    sp->len = 0;
    FPRINTF(g_mtxerr,"ERROR: (mtx.c)  mtx_steal_org_row_sparse called with\n");
    FPRINTF(g_mtxerr,"                sparse of insufficient capacity.\n");
    return TRUE;
  }

  last_value_matrix = mtx;
  org_row = mtx->perm.row.cur_to_org[row];
  rlink = &(mtx->hdr.row[org_row]);
  k = 0;

  if( rng == mtx_ALL_COLS ) {
    while( NOTNULL(*rlink)  ) {
      sp->idata[k] = (*rlink)->col;
      sp->data[k] = (*rlink)->value;
      k++;
      clink = &(mtx->hdr.col[(*rlink)->col]);
      for( ; NOTNULL(*clink) && (*clink)->row != org_row ; )
        clink = &((*clink)->next.row);
      delete_from_row(rlink);
      delete_from_col(clink);
    }
  } else if( rng->high >= rng->low ) {
    int32 *tocur;
    tocur = mtx->perm.col.org_to_cur;
    while( NOTNULL(*rlink) ) {
      if( in_range(rng,tocur[(*rlink)->col]) ) {
        sp->idata[k] = (*rlink)->col;
        sp->data[k] = (*rlink)->value;
        k++;
        clink = &(mtx->hdr.col[(*rlink)->col]);
        for( ; NOTNULL(*clink) && (*clink)->row != org_row ; )
          clink = &((*clink)->next.row);
        delete_from_row(rlink);
        delete_from_col(clink);
      } else {
         rlink = &((*rlink)->next.col);
      }
    }
  }
  sp->len = k;
  return FALSE;
}

boolean mtx_steal_org_col_sparse(mtx_matrix_t mtx, int32 col,
                                 mtx_sparse_t *sp, mtx_range_t *rng)
{
  struct element_t **clink, **rlink;
  int32 org_col;
  int32 len,k;

#if MTX_DEBUG
  if( !mtx_check_matrix(mtx) ) return TRUE;
#endif
  if (sp == mtx_CREATE_SPARSE) {
    FPRINTF(g_mtxerr,"ERROR: (mtx.c)  mtx_steal_org_col_sparse called with\n");
    FPRINTF(g_mtxerr,"                mtx_CREATE_SPARSE. Not supported.\n");
    return TRUE;
  }
  if (rng == mtx_ALL_ROWS) {
    len = mtx->order;
  } else {
    len = rng->high - rng->low +1;
  }
  if (sp->cap < len) {
    sp->len = 0;
    FPRINTF(g_mtxerr,"ERROR: (mtx.c)  mtx_steal_org_col_sparse called with\n");
    FPRINTF(g_mtxerr,"                sparse of insufficient capacity.\n");
    return TRUE;
  }

  last_value_matrix = mtx;
  org_col = mtx->perm.col.cur_to_org[col];
  clink = &(mtx->hdr.col[org_col]);
  k = 0;

  if( rng == mtx_ALL_ROWS ) {
    while( NOTNULL(*clink) ) {
      sp->idata[k] = (*clink)->row;
      sp->data[k] = (*clink)->value;
      k++;
      rlink = &(mtx->hdr.row[(*clink)->row]);
      for( ; NOTNULL(*rlink) && (*rlink)->col != org_col ; )
        rlink = &((*rlink)->next.col);
      delete_from_col(clink);
      delete_from_row(rlink);
    }
  } else if( rng->high >= rng->low ) {
    int32 *tocur;
    tocur = mtx->perm.row.org_to_cur;
    while( NOTNULL(*clink) ) {
      if( in_range(rng,tocur[(*clink)->row]) ) {
        sp->idata[k] = (*clink)->row;
        sp->data[k] = (*clink)->value;
        k++;
        rlink = &(mtx->hdr.row[(*clink)->row]);
        for( ; NOTNULL(*rlink) && (*rlink)->col != org_col ; )
          rlink = &((*rlink)->next.col);
        delete_from_col(clink);
        delete_from_row(rlink);
      } else {
        clink = &((*clink)->next.row);
      }
    }
  }
  sp->len = k;
  return FALSE;
}

boolean mtx_steal_cur_row_sparse(mtx_matrix_t mtx, int32 row,
                                 mtx_sparse_t *sp, mtx_range_t *rng)
{
  struct element_t **rlink, **clink;
  int32 org_row, *tocur;
  int32 len,k;

#if MTX_DEBUG
  if( !mtx_check_matrix(mtx) ) return TRUE;
#endif
  if (sp == mtx_CREATE_SPARSE) {
    FPRINTF(g_mtxerr,"ERROR: (mtx.c)  mtx_steal_cur_row_sparse called with\n");
    FPRINTF(g_mtxerr,"                mtx_CREATE_SPARSE. Not supported.\n");
    return TRUE;
  }
  if (rng == mtx_ALL_COLS) {
    len = mtx->order;
  } else {
    len = rng->high - rng->low +1;
  }
  if (sp->cap < len) {
    sp->len = 0;
    FPRINTF(g_mtxerr,"ERROR: (mtx.c)  mtx_steal_cur_row_sparse called with\n");
    FPRINTF(g_mtxerr,"                sparse of insufficient capacity.\n");
    return TRUE;
  }

  tocur = mtx->perm.col.org_to_cur;
  last_value_matrix = mtx;
  org_row = mtx->perm.row.cur_to_org[row];
  rlink = &(mtx->hdr.row[org_row]);
  k = 0;

  if( rng == mtx_ALL_COLS ) {
    while( NOTNULL(*rlink)  ) {
      sp->idata[k] = tocur[(*rlink)->col];
      sp->data[k] = (*rlink)->value;
      k++;
      clink = &(mtx->hdr.col[(*rlink)->col]);
      for( ; NOTNULL(*clink) && (*clink)->row != org_row ; )
        clink = &((*clink)->next.row);
      delete_from_row(rlink);
      delete_from_col(clink);
    }
  } else if( rng->high >= rng->low ) {
    while( NOTNULL(*rlink) ) {
      if( in_range(rng,tocur[(*rlink)->col]) ) {
        sp->idata[k] = tocur[(*rlink)->col];
        sp->data[k] = (*rlink)->value;
        k++;
        clink = &(mtx->hdr.col[(*rlink)->col]);
        for( ; NOTNULL(*clink) && (*clink)->row != org_row ; )
          clink = &((*clink)->next.row);
        delete_from_row(rlink);
        delete_from_col(clink);
      } else {
         rlink = &((*rlink)->next.col);
      }
    }
  }
  sp->len = k;
  return FALSE;
}

boolean mtx_steal_cur_col_sparse(mtx_matrix_t mtx, int32 col,
                                 mtx_sparse_t *sp, mtx_range_t *rng)
{
  struct element_t **clink, **rlink;
  int32 org_col, *tocur;
  int32 len,k;

#if MTX_DEBUG
  if( !mtx_check_matrix(mtx) ) return TRUE;
#endif
  if (sp == mtx_CREATE_SPARSE) {
    FPRINTF(g_mtxerr,"ERROR: (mtx.c)  mtx_steal_cur_col_sparse called with\n");
    FPRINTF(g_mtxerr,"                mtx_CREATE_SPARSE. Not supported.\n");
    return TRUE;
  }
  if (rng == mtx_ALL_ROWS) {
    len = mtx->order;
  } else {
    len = rng->high - rng->low +1;
  }
  if (sp->cap < len) {
    sp->len = 0;
    FPRINTF(g_mtxerr,"ERROR: (mtx.c)  mtx_steal_cur_col_sparse called with\n");
    FPRINTF(g_mtxerr,"                sparse of insufficient capacity.\n");
    return TRUE;
  }

  tocur = mtx->perm.row.org_to_cur;
  last_value_matrix = mtx;
  org_col = mtx->perm.col.cur_to_org[col];
  clink = &(mtx->hdr.col[org_col]);
  k = 0;

  if( rng == mtx_ALL_ROWS ) {
    while( NOTNULL(*clink) ) {
      sp->idata[k] = tocur[(*clink)->row];
      sp->data[k] = (*clink)->value;
      k++;
      rlink = &(mtx->hdr.row[(*clink)->row]);
      for( ; NOTNULL(*rlink) && (*rlink)->col != org_col ; )
        rlink = &((*rlink)->next.col);
      delete_from_col(clink);
      delete_from_row(rlink);
    }
  } else if( rng->high >= rng->low ) {
    while( NOTNULL(*clink) ) {
      if( in_range(rng,tocur[(*clink)->row]) ) {
        sp->idata[k] = tocur[(*clink)->row];
        sp->data[k] = (*clink)->value;
        k++;
        rlink = &(mtx->hdr.row[(*clink)->row]);
        for( ; NOTNULL(*rlink) && (*rlink)->col != org_col ; )
          rlink = &((*rlink)->next.col);
        delete_from_col(clink);
        delete_from_row(rlink);
      } else {
        clink = &((*clink)->next.row);
      }
    }
  }
  sp->len = k;
  return FALSE;
}

void mtx_fill_org_row_vec(mtx_matrix_t mtx, int32 row,
                          real64 *vec, mtx_range_t *rng)
{
  int32 org_row,highcol, *toorg;
  register int32 org_col;
  register real64 value;

#if MTX_DEBUG
  if(!mtx_check_matrix(mtx)) return;
#endif

  org_row = mtx->perm.row.cur_to_org[row];
  if( rng == mtx_ALL_COLS ) {
    highcol=mtx->order;
    for( org_col=0 ; org_col <highcol ; org_col++ ) {
      if ((value=vec[org_col]) != D_ZERO)
        mtx_create_element_value(mtx,org_row,org_col,value);
    }
  } else {
    register int32 cur_col;

    toorg = mtx->perm.col.cur_to_org;
    highcol= rng->high;
    for ( cur_col=rng->low; cur_col<=highcol; cur_col++) {
      if ((value=vec[(org_col=toorg[cur_col])]) != D_ZERO)
        mtx_create_element_value(mtx,org_row,org_col,value);
    }
  }
}

void mtx_fill_org_col_vec(mtx_matrix_t mtx, int32 col,
                          real64 *vec, mtx_range_t *rng)
{
  int32 org_col,highrow, *toorg;
  register int32 org_row;
  register real64 value;

#if MTX_DEBUG
  if(!mtx_check_matrix(mtx)) return;
#endif

  org_col = mtx->perm.col.cur_to_org[col];
  if( rng == mtx_ALL_ROWS ) {
    highrow=mtx->order;
    for( org_row=0 ; org_row <highrow ; org_row++ ) {
      if ((value=vec[org_row]) != D_ZERO)
        mtx_create_element_value(mtx,org_row,org_col,value);
    }
  } else {
    register int32 cur_row;

    toorg = mtx->perm.row.cur_to_org;
    highrow= rng->high;
    for ( cur_row=rng->low; cur_row<=highrow; cur_row++) {
      if ((value=vec[(org_row=toorg[cur_row])]) != D_ZERO)
        mtx_create_element_value(mtx,org_row,org_col,value);
    }
  }
}

void mtx_fill_cur_row_vec(mtx_matrix_t mtx, int32 row,
                          real64 *vec, mtx_range_t *rng)
{
  int32 org_row,highcol,lowcol, *toorg;
  register int32 cur_col;
  register real64 value;

#if MTX_DEBUG
  if(!mtx_check_matrix(mtx)) return;
#endif

  org_row = mtx->perm.row.cur_to_org[row];
  toorg=mtx->perm.col.cur_to_org;
  if( rng == mtx_ALL_COLS ) {
    highcol = mtx->order-1;
    lowcol = 0;
  } else {
    highcol= rng->high;
    lowcol=rng->low;
  }
  for( cur_col=lowcol ; cur_col <= highcol ; cur_col++ ) {
    if ((value=vec[cur_col]) != D_ZERO)
      mtx_create_element_value(mtx,org_row,toorg[cur_col],value);
  }
}

void mtx_fill_cur_col_vec(mtx_matrix_t mtx, int32 col,
                          real64 *vec, mtx_range_t *rng)
{
  int32 org_col,highrow,lowrow, *toorg;
  register int32 cur_row;
  register real64 value;

#if MTX_DEBUG
  if(!mtx_check_matrix(mtx)) return;
#endif

  org_col = mtx->perm.col.cur_to_org[col];
  toorg=mtx->perm.row.cur_to_org;
  if( rng == mtx_ALL_ROWS ) {
    highrow=mtx->order-1;
    lowrow=0;
  } else {
    highrow= rng->high;
    lowrow=rng->low;
  }
  for( cur_row=lowrow ; cur_row <= highrow ; cur_row++ ) {
    if ((value=vec[cur_row]) != D_ZERO)
      mtx_create_element_value(mtx,toorg[cur_row],org_col,value);
  }
}

void mtx_dropfill_cur_col_vec(mtx_matrix_t mtx, int32 col,
                              real64 *vec, mtx_range_t *rng,
                              real64 tol)
{
  int32 org_col,highrow,lowrow, *toorg;
  register int32 cur_row;
  register real64 value;

  if (tol==0.0) {
    mtx_fill_cur_col_vec(mtx,col,vec,rng);
    return;
  }

#if MTX_DEBUG
  if(!mtx_check_matrix(mtx)) return;
#endif

  org_col = mtx->perm.col.cur_to_org[col];
  toorg=mtx->perm.row.cur_to_org;
  if( rng == mtx_ALL_ROWS ) {
    highrow=mtx->order-1;
    lowrow=0;
  } else {
    highrow= rng->high;
    lowrow=rng->low;
  }
  for( cur_row=lowrow ; cur_row <= highrow ; cur_row++ ) {
    if (fabs(value=vec[cur_row]) > tol)
      mtx_create_element_value(mtx,toorg[cur_row],org_col,value);
  }
}

void mtx_dropfill_cur_row_vec(mtx_matrix_t mtx, int32 row,
                              real64 *vec, mtx_range_t *rng,
                              real64 tol)
{
  int32 org_row,highcol,lowcol, *toorg;
  register int32 cur_col;
  register real64 value;

  if (tol==0.0) {
    mtx_fill_cur_row_vec(mtx,row,vec,rng);
    return;
  }

#if MTX_DEBUG
  if(!mtx_check_matrix(mtx)) return;
#endif

  org_row = mtx->perm.row.cur_to_org[row];
  toorg=mtx->perm.col.cur_to_org;
  if( rng == mtx_ALL_COLS ) {
    highcol = mtx->order-1;
    lowcol = 0;
  } else {
    highcol= rng->high;
    lowcol=rng->low;
  }
  for( cur_col=lowcol ; cur_col <= highcol ; cur_col++ ) {
    if (fabs(value=vec[cur_col]) > tol)
      mtx_create_element_value(mtx,org_row,toorg[cur_col],value);
  }
}

void mtx_fill_org_row_sparse(mtx_matrix_t mtx, int32 row,
                             const mtx_sparse_t *sp)
{
  int32 orgrow,i;
#if MTX_DEBUG
  if(!mtx_check_matrix(mtx) || !mtx_check_sparse(sp)) return;
#endif
  orgrow = mtx->perm.row.cur_to_org[row];
  for (i=0; i < sp->len; i++) {
    if (sp->data[i] != D_ZERO
#if MTX_DEBUG
        && sp->idata[i] >=0 && sp->idata[i] < mtx->order
#endif
       ) {
      mtx_create_element_value(mtx,orgrow,sp->idata[i],sp->data[i]);
    }
  }
}

void mtx_fill_org_col_sparse(mtx_matrix_t mtx, int32 col,
                             const mtx_sparse_t *sp)
{
  int32 orgcol,i;
#if MTX_DEBUG
  if(!mtx_check_matrix(mtx) || !mtx_check_sparse(sp)) return;
#endif
  orgcol = mtx->perm.col.cur_to_org[col];
  for (i=0; i < sp->len; i++) {
    if (sp->data[i] != D_ZERO
#if MTX_DEBUG
        && sp->idata[i] >=0 && sp->idata[i] < mtx->order
#endif
       ) {
      mtx_create_element_value(mtx,sp->idata[i],orgcol,sp->data[i]);
    }
  }
}

void mtx_fill_cur_row_sparse(mtx_matrix_t mtx, int32 row,
                             const mtx_sparse_t *sp)
{
  int32 orgrow,i;
#if MTX_DEBUG
  if(!mtx_check_matrix(mtx) || !mtx_check_sparse(sp)) return;
#endif
  orgrow = mtx->perm.row.cur_to_org[row];
  for (i=0; i < sp->len; i++) {
    if (
#if MTX_DEBUG
        sp->idata[i] >=0 && sp->idata[i] < mtx->order &&
#endif
        sp->data[i] != D_ZERO) {
      mtx_create_element_value(mtx,orgrow,
                           mtx->perm.col.cur_to_org[sp->idata[i]],
                           sp->data[i]);
    }
  }
}

void mtx_fill_cur_col_sparse(mtx_matrix_t mtx, int32 col,
                             const mtx_sparse_t *sp)
{
  int32 orgcol,i;
#if MTX_DEBUG
  if(!mtx_check_matrix(mtx) || !mtx_check_sparse(sp)) return;
#endif
  orgcol = mtx->perm.col.cur_to_org[col];
  for (i=0; i < sp->len; i++) {
    if (
#if MTX_DEBUG
        sp->idata[i] >=0 && sp->idata[i] < mtx->order &&
#endif
        sp->data[i] != D_ZERO) {
      mtx_create_element_value(mtx,
                           mtx->perm.col.cur_to_org[sp->idata[i]],
                           orgcol,
                           sp->data[i]);
    }
  }
}

void mtx_mult_row(mtx_matrix_t mtx, int32 row, real64 factor, mtx_range_t *rng)
{
   struct element_t *elt;
   int32 *tocur;

   if( factor == D_ZERO ) {
      mtx_clear_row(mtx,row,rng);
      return;
   }

   if (factor == D_ONE) return;
#if MTX_DEBUG
   if(!mtx_check_matrix(mtx)) return;
#endif

   tocur = mtx->perm.col.org_to_cur;
   elt = mtx->hdr.row[mtx->perm.row.cur_to_org[row]];
   if( rng == mtx_ALL_COLS )
      for( ; NOTNULL(elt); elt = elt->next.col ) elt->value *= factor;
   else if( rng->high >= rng->low )
      for( ; NOTNULL(elt); elt = elt->next.col )
         if( in_range(rng,tocur[elt->col]) ) elt->value *= factor;
}

void mtx_mult_col(mtx_matrix_t mtx, int32 col, real64 factor, mtx_range_t *rng)
{
   struct element_t *elt;
   int32 *tocur;

   if( factor == D_ZERO ) {
      mtx_clear_col(mtx,col,rng);
      return;
   }
   if (factor == D_ONE) return;

#if MTX_DEBUG
   if(!mtx_check_matrix(mtx)) return;
#endif

   tocur = mtx->perm.row.org_to_cur;
   elt = mtx->hdr.col[mtx->perm.col.cur_to_org[col]];
   if( rng == mtx_ALL_ROWS )
      for( ; NOTNULL(elt); elt = elt->next.row ) elt->value *= factor;
   else if( rng->high >= rng->low )
      for( ; NOTNULL(elt); elt = elt->next.row )
         if( in_range(rng,tocur[elt->row]) ) elt->value *= factor;
}

void mtx_mult_row_zero(mtx_matrix_t mtx, int32 row, mtx_range_t *rng)
{
   struct element_t *elt;
   int32 *tocur;

#if MTX_DEBUG
   if(!mtx_check_matrix(mtx)) return;
#endif

   tocur = mtx->perm.col.org_to_cur;
   elt = mtx->hdr.row[mtx->perm.row.cur_to_org[row]];
   if( rng == mtx_ALL_COLS )
      for( ; NOTNULL(elt); elt = elt->next.col ) elt->value = D_ZERO;
   else if( rng->high >= rng->low )
      for( ; NOTNULL(elt); elt = elt->next.col )
         if( in_range(rng,tocur[elt->col]) ) elt->value = D_ZERO;
}

void mtx_mult_col_zero(mtx_matrix_t mtx, int32 col, mtx_range_t *rng)
{
   struct element_t *elt;
   int32 *tocur;

#if MTX_DEBUG
   if(!mtx_check_matrix(mtx)) return;
#endif

   tocur = mtx->perm.row.org_to_cur;
   elt = mtx->hdr.col[mtx->perm.col.cur_to_org[col]];
   if( rng == mtx_ALL_ROWS )
      for( ; NOTNULL(elt); elt = elt->next.row ) elt->value = D_ZERO;
   else if( rng->high >= rng->low )
      for( ; NOTNULL(elt); elt = elt->next.row )
         if( in_range(rng,tocur[elt->row]) ) elt->value = D_ZERO;
}

void mtx_add_row(mtx_matrix_t mtx, int32 s_cur, int32 t_cur, real64 factor,
        mtx_range_t *rng)
{
   register int32 org_col;
   int32 t_org,*tocur;
   struct element_t **arr,*elt;

#if MTX_DEBUG
   if(!mtx_check_matrix(mtx)) return;
#endif

   t_org = mtx->perm.row.cur_to_org[t_cur];
   if( rng == mtx_ALL_COLS ) {
      arr = mtx_expand_row(mtx,t_org);   /* Expand the target row */
      elt = mtx->hdr.row[mtx->perm.row.cur_to_org[s_cur]];
      for( ; NOTNULL(elt); elt = elt->next.col ) {
         if( ISNULL(arr[(org_col=elt->col)]) ) {
            arr[org_col] =
             mtx_create_element_value(mtx,t_org,org_col,(factor * elt->value));
         } else {
            arr[org_col]->value += factor * elt->value;
         }
      }
      mtx_renull_using_row(mtx,t_org,arr);
   } else if( rng->high >= rng->low ) {
      register int32 cur_col;

      tocur = mtx->perm.col.org_to_cur;
      arr = mtx_expand_row(mtx,t_org);   /* Expand the target row */
      elt = mtx->hdr.row[mtx->perm.row.cur_to_org[s_cur]];
      for( ; NOTNULL(elt); elt = elt->next.col ) {
         cur_col=tocur[(org_col=elt->col)];
         if( in_range(rng,cur_col) )
            if( NOTNULL(arr[org_col]) ) {
               arr[org_col]->value += factor * elt->value;
            } else {
              arr[org_col] =
                mtx_create_element_value(mtx,t_org,org_col,
                                         (factor * elt->value));
               /* hit rate usually lower */
            }
      }
      mtx_renull_using_row(mtx,t_org,arr);
   }
   mtx_null_vector_release();
}

void mtx_add_col(mtx_matrix_t mtx, int32 s_cur, int32 t_cur, real64 factor,
         mtx_range_t *rng)
{
   int32 t_org,*tocur;
   struct element_t **arr,*elt;

#if MTX_DEBUG
   if(!mtx_check_matrix(mtx)) return;
#endif

   if( rng == mtx_ALL_ROWS ) {
      t_org = mtx->perm.col.cur_to_org[t_cur];
      arr = mtx_expand_col(mtx,t_org);   /* Expand the target col */
      elt = mtx->hdr.col[mtx->perm.col.cur_to_org[s_cur]];
      for( ; NOTNULL(elt); elt = elt->next.row )
         if( ISNULL(arr[elt->row]) ) {
            arr[elt->row] =
               mtx_create_element_value(mtx,elt->row,t_org,
                                        (factor * elt->value));
         } else arr[elt->row]->value += factor * elt->value;
      mtx_renull_using_col(mtx,t_org,arr);
   } else if( rng->high >= rng->low ) {
      tocur = mtx->perm.row.org_to_cur;
      t_org = mtx->perm.col.cur_to_org[t_cur];
      arr = mtx_expand_col(mtx,t_org);   /* Expand the target col */
      elt = mtx->hdr.col[mtx->perm.col.cur_to_org[s_cur]];
      for( ; NOTNULL(elt); elt = elt->next.row )
         if( in_range(rng,tocur[elt->row]) )
            if( NOTNULL(arr[elt->row]) ) {
               arr[elt->row]->value += factor * elt->value;
            } else {
               arr[elt->row] =
                 mtx_create_element_value(mtx,elt->row,t_org,
                                          (factor * elt->value));
            }
      mtx_renull_using_col(mtx,t_org,arr);
   }
   mtx_null_vector_release();
}

static struct element_t **mtx_expand_row_series( mtx_matrix_t mtx, int32 org)
/**
 ***  Expands the given row into an array of pointers, indexed on original
 ***  col number.  The array is obtained from mtx_null_row_vector().
 ***  Be sure to call mtx_null_row_vector_release() when done with the vector and
 ***  you have rezeroed it.
 **/
{
   struct element_t **arr;
   struct element_t *elt;

   arr = mtx_null_row_vector(mtx->order);
   for( elt=mtx->hdr.row[org]; NOTNULL(elt); elt=elt->next.col )
      arr[elt->col] = elt;
   return(arr);
}

static struct element_t **mtx_expand_col_series(mtx_matrix_t mtx, int32 org)
/**
 ***  Expands the given col into an array of pointers, indexed on original
 ***  row number.  The array is obtained from mtx_null_col_vector().
 *** Be sure to call mtx_null_col_vector_release() when done with the vector and
 ***  you have rezeroed it.
 **/
{
   struct element_t **arr;
   struct element_t *elt;

   arr = mtx_null_col_vector(mtx->order);
   for( elt = mtx->hdr.col[org] ; NOTNULL(elt) ; elt = elt->next.row )
      arr[elt->row] = elt;
   return(arr);
}

struct add_series_data {
  mtx_matrix_t mtx;        /* matrix we're operating on */
  struct element_t **arr;  /* target row/col expansion array */
  int32 *tocur;      /* col/row permutation vector */
  int32 t_org;       /* org row/col which is expanded */
};

static struct add_series_data
  rsdata={NULL,NULL,NULL,mtx_NONE},
  csdata={NULL,NULL,NULL,mtx_NONE};

static void add_series_data_release(void)
{
  /* if apparently sane, and have array release the arr */
  if (NOTNULL(rsdata.mtx) && NOTNULL(rsdata.tocur)
      && rsdata.t_org >= 0 && NOTNULL(rsdata.arr)) {
    mtx_null_row_vector_release();
  }
  if (NOTNULL(csdata.mtx) && NOTNULL(csdata.tocur)
      && csdata.t_org >= 0 && NOTNULL(csdata.arr)) {
    mtx_null_col_vector_release();
  }
  /* reinit data */

  rsdata.mtx = NULL;
  rsdata.arr = NULL;
  rsdata.tocur = NULL;
  rsdata.t_org = mtx_NONE;

  csdata.mtx = NULL;
  csdata.arr = NULL;
  csdata.tocur = NULL;
  csdata.t_org = mtx_NONE;
}

void mtx_add_row_series_init(mtx_matrix_t mtx,int32 t_cur,boolean use)
{

#if MTX_DEBUG
  if(!mtx_check_matrix(mtx)) return;
#endif

  if ( !(rsdata.mtx) ) {
    /* this is a grabbing call due to memory reuse */
    if ( mtx && (t_cur >= 0) && (t_cur < mtx->order) ) {
      rsdata.mtx   = mtx;
      rsdata.tocur = mtx->perm.col.org_to_cur;
      rsdata.t_org = mtx->perm.row.cur_to_org[t_cur];
      rsdata.arr   = mtx_expand_row_series(mtx,rsdata.t_org);
    } else {
      FPRINTF(g_mtxerr,"ERROR:     mtx_add_row_series_init called\n");
      FPRINTF(g_mtxerr,"           for grab with invalid column or mtx.\n");
      FPRINTF(g_mtxerr,"           col number %d.\n",t_cur);
    }
    return;
  } else {
    /* this is supposed to be a releasing call */
    if (t_cur != mtx_NONE) {
      FPRINTF(g_mtxerr,"ERROR:     mtx_add_row_series_init called for\n");
      FPRINTF(g_mtxerr,"           release without mtx_NONE.\n");
      return;
    }
    if (mtx != rsdata.mtx) {
      FPRINTF(g_mtxerr,"ERROR:     mtx_add_row_series_init called for\n");
      FPRINTF(g_mtxerr,"           release with ungrabbed matrix.\n");
      return;
    }
    if (use) {
      mtx_renull_using_row(rsdata.mtx, rsdata.t_org, rsdata.arr);
    } else {
      mtx_renull_all(rsdata.mtx, rsdata.arr);
    }
    mtx_null_row_vector_release();
    rsdata.mtx   = NULL;
    rsdata.arr   = NULL;
    rsdata.tocur = NULL;
    rsdata.t_org = mtx_NONE;
  }
}

void mtx_add_row_series(int32 s_cur,real64 factor,mtx_range_t *rng)
{
  register int32 org_col;
  int32 t_org,*tocur;
  struct element_t **arr,*elt;
  mtx_matrix_t mtx;

  if ( !(rsdata.mtx) ) {
    FPRINTF(g_mtxerr,"ERROR:     mtx_add_row_series called for\n");
    FPRINTF(g_mtxerr,"           without grabbing target row first.\n");
    return;
  }
  mtx   = rsdata.mtx;
  arr   = rsdata.arr;
  tocur = rsdata.tocur;
  t_org = rsdata.t_org;

#if MTX_DEBUG
  if(!mtx_check_matrix(mtx)) return;
#endif

  elt = mtx->hdr.row[mtx->perm.row.cur_to_org[s_cur]];
  if( rng == mtx_ALL_COLS ) {
    for( ; NOTNULL(elt); elt = elt->next.col ) {
      if( ISNULL(arr[(org_col=elt->col)]) ) {
        arr[org_col] =
          mtx_create_element_value(mtx,t_org,org_col,(factor * elt->value));
      } else {
        arr[org_col]->value += factor * elt->value;
      }
    }
    return;
  }
  /* fast_in_range is a 10% winner on the alpha, and should be even more
     on the other platforms. */
  if( rng->high >= rng->low ) {
    register int32 cur_col, lo,hi;
    lo=rng->low; hi= rng->high;
    for( ; NOTNULL(elt); elt = elt->next.col ) {
      cur_col=tocur[(org_col=elt->col)];
      if( fast_in_range(lo,hi,cur_col) ) {
        if( NOTNULL(arr[org_col]) ) {
          arr[org_col]->value += factor * elt->value;
        } else {
          arr[org_col] =
            mtx_create_element_value(mtx,t_org,org_col,(factor*elt->value));
        }
      }
    }
  }
}

void mtx_add_col_series_init(mtx_matrix_t mtx,int32 t_cur,boolean use)
{

#if MTX_DEBUG
  if(!mtx_check_matrix(mtx)) return;
#endif

  if ( !(csdata.mtx) ) {
    /* this is a grabbing call */
    if ( mtx && (t_cur >= 0) && (t_cur < mtx->order) ) {
      csdata.mtx   = mtx;
      csdata.tocur = mtx->perm.row.org_to_cur;
      csdata.t_org = mtx->perm.col.cur_to_org[t_cur];
      csdata.arr   = mtx_expand_col_series(mtx,csdata.t_org);
    } else {
      FPRINTF(g_mtxerr,"ERROR:     mtx_add_col_series_init called\n");
      FPRINTF(g_mtxerr,"           for grab with invalid row or mtx.\n");
      FPRINTF(g_mtxerr,"           row number %d.\n",t_cur);
    }
    return;
  } else {
    /* this is supposed to be a releasing call */
    if (t_cur != mtx_NONE) {
      FPRINTF(g_mtxerr,"ERROR:     mtx_add_col_series_init called for\n");
      FPRINTF(g_mtxerr,"           release without mtx_NONE.\n");
      return;
    }
    if (mtx != csdata.mtx) {
      FPRINTF(g_mtxerr,"ERROR:     mtx_add_col_series_init called for\n");
      FPRINTF(g_mtxerr,"           release with ungrabbed matrix.\n");
      return;
    }
    if (use) {
      mtx_renull_using_col(csdata.mtx, csdata.t_org, csdata.arr);
    } else {
      mtx_renull_all(csdata.mtx, csdata.arr);
    }
    mtx_null_col_vector_release();
    csdata.mtx   = NULL;
    csdata.arr   = NULL;
    csdata.tocur = NULL;
    csdata.t_org = mtx_NONE;
  }
}

void mtx_add_col_series(int32 s_cur,real64 factor,mtx_range_t *rng)
{
  register int32 org_row;
  int32 t_org,*tocur;
  struct element_t **arr,*elt;
  mtx_matrix_t mtx;

  if ( !(csdata.mtx) ) {
    FPRINTF(g_mtxerr,"ERROR:     mtx_add_col_series called for\n");
    FPRINTF(g_mtxerr,"           without grabbing target col first.\n");
    return;
  }
  mtx   = csdata.mtx;
  arr   = csdata.arr;
  tocur = csdata.tocur;
  t_org = csdata.t_org;

#if MTX_DEBUG
  if(!mtx_check_matrix(mtx)) return;
#endif

  elt = mtx->hdr.col[mtx->perm.col.cur_to_org[s_cur]];
  if( rng == mtx_ALL_ROWS ) {
    for( ; NOTNULL(elt); elt = elt->next.row ) {
      if( ISNULL(arr[(org_row=elt->row)]) ) {
        arr[org_row] =
          mtx_create_element_value(mtx,org_row,t_org,(factor * elt->value));
      } else {
        arr[org_row]->value += factor * elt->value;
      }
    }
    return;
  }
  if( rng->high >= rng->low ) {
    register int32 cur_row;

    for( ; NOTNULL(elt); elt = elt->next.row ) {
      cur_row=tocur[(org_row=elt->row)];
      if( in_range(rng,cur_row) ) {
        if( NOTNULL(arr[org_row]) ) {
          arr[org_row]->value += factor * elt->value;
        } else {
          arr[org_row] =
            mtx_create_element_value(mtx,org_row,t_org,(factor*elt->value));
        }
      }
    }
  }
}

void mtx_old_add_row_sparse(mtx_matrix_t mtx,
                       int32 t_cur,   /* cur index of target row */
                       real64 *s_cur, /* dense source row, curindexed */
                       real64 factor, /* coefficient */
                       mtx_range_t *rng,
                       int32 *ilist)  /* list to add over */
{
  int32 t_org,*toorg,cindex,orgcindex,hilim;
  struct element_t **arr;
  real64 value;

#if MTX_DEBUG
  if(!mtx_check_matrix(mtx)) return;
#endif

  if (factor==D_ZERO) return; /* adding 0 rather silly */

  if( rng == mtx_ALL_COLS ) {
    t_org = mtx->perm.row.cur_to_org[t_cur]; /* org of target row */
    arr = mtx_expand_row(mtx,t_org);             /* Expand the target row */
    toorg = mtx->perm.col.cur_to_org;        /* current col perm */
    hilim=mtx->order;

    if (!ilist) { /* FULL ROW CASE no ilist */
      for (cindex=0; cindex< hilim; cindex++) {
        if ((value=s_cur[cindex])!=D_ZERO) {
          if ( ISNULL(arr[(orgcindex=toorg[cindex])]) ) { /* add element */
            arr[orgcindex] =
              mtx_create_element_value(mtx,t_org,orgcindex,(factor * value));
          } else {  /* increment element */
            arr[orgcindex]->value += factor * value;
          }
        }
      }
    } else { /* SPARSE ROW CASE with ilist */
      int32 i;
      i=0;
      while ((cindex=ilist[i])>=0) {
        value=s_cur[cindex];
        if ( ISNULL(arr[(orgcindex=toorg[cindex])]) ) { /* add element */
          arr[orgcindex] =
            mtx_create_element_value(mtx,t_org,orgcindex,(factor * value));
        } else {  /* increment element */
          arr[orgcindex]->value += factor * value;
        }
        i++;
      }
    }
    mtx_renull_using_row(mtx,t_org,arr);
  } else if( rng->high >= rng->low ) { /* DENSE RANGE CASE */
    t_org = mtx->perm.row.cur_to_org[t_cur]; /* org of target row */
    arr = mtx_expand_row(mtx,t_org);             /* Expand the target row */
    toorg = mtx->perm.col.cur_to_org;        /* current col perm */
    hilim = rng->high;

    for (cindex=rng->low; cindex<= hilim; cindex++) {
      if ((value=s_cur[cindex])!=D_ZERO) {
        if ( ISNULL(arr[(orgcindex=toorg[cindex])]) ) { /* add element */
          arr[orgcindex] =
            mtx_create_element_value(mtx,t_org,orgcindex,(factor * value));
        } else {  /* increment element */
          arr[orgcindex]->value += factor * value;
        }
      }
    }
    mtx_renull_using_row(mtx,t_org,arr);
  }
  mtx_null_vector_release();
}

void mtx_old_add_col_sparse(mtx_matrix_t mtx,
                       int32 t_cur,   /* cur index of target col */
                       real64 *s_cur, /* dense source col, curindexed */
                       real64 factor, /* coefficient */
                       mtx_range_t *rng,    /* range to add over or */
                       int32 *ilist)  /* list to add over */
{
  int32 t_org,*toorg,rindex,orgrindex,hilim;
  struct element_t **arr;
  real64 value;

#if MTX_DEBUG
  if(!mtx_check_matrix(mtx)) return;
#endif

  if (factor==D_ZERO) return;

  if( rng == mtx_ALL_ROWS ) {
    t_org = mtx->perm.col.cur_to_org[t_cur]; /* org of target col */
    arr = mtx_expand_col(mtx,t_org);             /* Expand the target col */
    toorg = mtx->perm.row.cur_to_org;        /* current row perm */
    hilim=mtx->order;

    if (!ilist) { /* FULL COL CASE no ilist */
      for (rindex=0; rindex< hilim; rindex++) {
        if ((value=s_cur[rindex])!=D_ZERO) {
          if ( ISNULL(arr[(orgrindex=toorg[rindex])]) ) { /* add element */
            arr[orgrindex] =
              mtx_create_element_value(mtx,orgrindex,t_org,(factor * value));
          } else { /* increment element */
            arr[orgrindex]->value += factor * value;
          }
        }
      }
    } else { /* SPARSE COL CASE with ilist */
      int32 i;
      i=0;
      while ((rindex=ilist[i])>=0) {
        value=s_cur[rindex];
        if ( ISNULL(arr[(orgrindex=toorg[rindex])])  ) { /* add element */
          arr[orgrindex] =
            mtx_create_element_value(mtx,orgrindex,t_org,(factor * value));
        } else { /* increment element */
          arr[orgrindex]->value += factor * value;
        }
        i++;
      }
    }
    mtx_renull_using_col(mtx,t_org,arr);

  } else if( rng->high >= rng->low ) { /* DENSE RANGE CASE */
    t_org = mtx->perm.col.cur_to_org[t_cur]; /* org of target col */
    arr = mtx_expand_col(mtx,t_org);             /* Expand the target col */
    toorg = mtx->perm.row.cur_to_org;        /* current row perm */
    hilim = rng->high;

    for (rindex=rng->low; rindex<= hilim; rindex++) {
      if ((value=s_cur[rindex])!=D_ZERO) {
        if ( ISNULL(arr[(orgrindex=toorg[rindex])]) ) { /* add element */
          arr[orgrindex] =
            mtx_create_element_value(mtx,orgrindex,t_org,(factor * value));
        } else { /* increment element */
          arr[orgrindex]->value += factor * value;
        }
      }
    }
    mtx_renull_using_col(mtx,t_org,arr);
  }
  mtx_null_vector_release();
}

size_t mtx_size(mtx_matrix_t mtx) {
  size_t size=0;
  if (ISNULL(mtx)) return size;
  size += (1+mtx->nslaves)*sizeof(struct mtx_header);
  if (ISSLAVE(mtx)) {
    return 2*sizeof(struct element_t *)*(size_t)mtx->capacity +size;
  }
  size += mem_sizeof_store(mtx->ms);
  /* headers */
  size += (1 + mtx->nslaves) * (size_t)mtx->capacity *
    (size_t)2 * sizeof(struct element_t *);
  /* block data */
  if (mtx->data->nblocks >0)
    size += sizeof(mtx_region_t)*(size_t)(mtx->data->nblocks-1);
  /* permutations */
  size += (size_t)4*(sizeof(int32)*(size_t)mtx->capacity);
  return size;
}

size_t mtx_chattel_size(mtx_matrix_t mtx) {
  size_t size=0;
  int32 i;
  if (ISNULL(mtx) || ISSLAVE(mtx)) return size;
  /* headers */
  size += (mtx->nslaves)*sizeof(struct mtx_header);
  /* incidence */
  for (i=0; i <mtx->nslaves; i++) {
    size += (size_t)mtx_nonzeros_in_region(mtx->slaves[i],mtx_ENTIRE_MATRIX);
  }
  /*nz headers */
  size += mtx->nslaves * (size_t)mtx->capacity *
    (