sparse.c

Go to the documentation of this file.

#include <../../nrnconf.h>
 
/**************************************************************************
**
** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
**
**               Meschach Library
** 
** This Meschach Library is provided "as is" without any express 
** or implied warranty of any kind with respect to this software. 
** In particular the authors shall not be liable for any direct, 
** indirect, special, incidental or consequential damages arising 
** in any way from use of the software.
** 
** Everyone is granted permission to copy, modify and redistribute this
** Meschach Library, provided:
**  1.  All copies contain this copyright notice.
**  2.  All modified copies shall carry a notice stating who
**      made the last modification and the date of such modification.
**  3.  No charge is made for this software or works derived from it.  
**      This clause shall not be construed as constraining other software
**      distributed on the same medium as this software, nor is a
**      distribution fee considered a charge.
**
***************************************************************************/
 
/*
  Sparse matrix package
  See also: sparse.h, matrix.h
  */
 
#include    <stdio.h>
#include    <math.h>
#include        <stdlib.h>
#include    "sparse.h"
 
 
static char rcsid[] = "sparse.c,v 1.1 1997/12/04 17:55:48 hines Exp";
 
#define MINROWLEN   10
 
 
 
/* sp_get_val -- returns the (i,j) entry of the sparse matrix A */
double  sp_get_val(A,i,j)
SPMAT   *A;
int i, j;
{
   SPROW    *r;
   int  idx;
   
   if ( A == SMNULL )
     error(E_NULL,"sp_get_val");
   if ( i < 0 || i >= A->m || j < 0 || j >= A->n )
     error(E_SIZES,"sp_get_val");
   
   r = A->row+i;
   idx = sprow_idx(r,j);
   if ( idx < 0 )
     return 0.0;
   /* else */
   return r->elt[idx].val;
}
 
/* sp_set_val -- sets the (i,j) entry of the sparse matrix A */
double  sp_set_val(A,i,j,val)
SPMAT   *A;
int i, j;
double  val;
{
   SPROW    *r;
   int  idx, idx2, new_len;
   
   if ( A == SMNULL )
     error(E_NULL,"sp_set_val");
   if ( i < 0 || i >= A->m || j < 0 || j >= A->n )
     error(E_SIZES,"sp_set_val");
   
   r = A->row+i;
   idx = sprow_idx(r,j);
   /* printf("sp_set_val: idx = %d\n",idx); */
   if ( idx >= 0 )
   {    r->elt[idx].val = val;  return val; }
   /* else */ if ( idx < -1 )
   {
      /* Note: this destroys the column & diag access paths */
      A->flag_col = A->flag_diag = FALSE;
      /* shift & insert new value */
      idx = -(idx+2);   /* this is the intended insertion index */
      if ( r->len >= r->maxlen )
      {
     r->len = r->maxlen;
     new_len = max(2*r->maxlen+1,5);
     if (mem_info_is_on()) {
        mem_bytes(TYPE_SPMAT,A->row[i].maxlen*sizeof(row_elt),
                new_len*sizeof(row_elt));
     }
 
     r->elt = RENEW(r->elt,new_len,row_elt);
     if ( ! r->elt )    /* can't allocate */
       error(E_MEM,"sp_set_val");
     r->maxlen = 2*r->maxlen+1;
      }
      for ( idx2 = r->len-1; idx2 >= idx; idx2-- )
    MEM_COPY((char *)(&(r->elt[idx2])),
         (char *)(&(r->elt[idx2+1])),sizeof(row_elt));
      /************************************************************
    if ( idx < r->len )
    MEM_COPY((char *)(&(r->elt[idx])),(char *)(&(r->elt[idx+1])),
    (r->len-idx)*sizeof(row_elt));
    ************************************************************/
      r->len++;
      r->elt[idx].col = j;
      return r->elt[idx].val = val;
   }
   /* else -- idx == -1, error in index/matrix! */
   return 0.0;
}
 
/* sp_mv_mlt -- sparse matrix/dense vector multiply
   -- result is in out, which is returned unless out==NULL on entry
   --  if out==NULL on entry then the result vector is created */
VEC *sp_mv_mlt(A,x,out)
SPMAT   *A;
VEC *x, *out;
{
   int  i, j_idx, m, n, max_idx;
   Real sum, *x_ve;
   SPROW    *r;
   row_elt  *elts;
   
   if ( ! A || ! x )
     error(E_NULL,"sp_mv_mlt");
   if ( x->dim != A->n )
     error(E_SIZES,"sp_mv_mlt");
   if ( ! out || out->dim < A->m )
     out = v_resize(out,A->m);
   if ( out == x )
     error(E_INSITU,"sp_mv_mlt");
   m = A->m;    n = A->n;
   x_ve = x->ve;
   
   for ( i = 0; i < m; i++ )
   {
      sum = 0.0;
      r = &(A->row[i]);
      max_idx = r->len;
      elts    = r->elt;
      for ( j_idx = 0; j_idx < max_idx; j_idx++, elts++ )
    sum += elts->val*x_ve[elts->col];
      out->ve[i] = sum;
   }
   return out;
}
 
/* sp_vm_mlt -- sparse matrix/dense vector multiply from left
   -- result is in out, which is returned unless out==NULL on entry
   -- if out==NULL on entry then result vector is created & returned */
VEC *sp_vm_mlt(A,x,out)
SPMAT   *A;
VEC *x, *out;
{
   int  i, j_idx, m, n, max_idx;
   Real tmp, *x_ve, *out_ve;
   SPROW    *r;
   row_elt  *elts;
   
   if ( ! A || ! x )
     error(E_NULL,"sp_vm_mlt");
   if ( x->dim != A->m )
     error(E_SIZES,"sp_vm_mlt");
   if ( ! out || out->dim < A->n )
     out = v_resize(out,A->n);
   if ( out == x )
     error(E_INSITU,"sp_vm_mlt");
   
   m = A->m;    n = A->n;
   v_zero(out);
   x_ve = x->ve;    out_ve = out->ve;
   
   for ( i = 0; i < m; i++ )
   {
      r = A->row+i;
      max_idx = r->len;
      elts    = r->elt;
      tmp = x_ve[i];
      for ( j_idx = 0; j_idx < max_idx; j_idx++, elts++ )
    out_ve[elts->col] += elts->val*tmp;
   }
   
   return out;
}
 
 
/* sp_get -- get sparse matrix
   -- len is number of elements available for each row without
   allocating further memory */
SPMAT   *sp_get(m,n,maxlen)
int m, n, maxlen;
{
   SPMAT    *A;
   SPROW    *rows;
   int  i;
   
   if ( m < 0 || n < 0 )
     error(E_NEG,"sp_get");
 
   maxlen = max(maxlen,1);
   
   A = NEW(SPMAT);
   if ( ! A )       /* can't allocate */
     error(E_MEM,"sp_get");
   else if (mem_info_is_on()) {
      mem_bytes(TYPE_SPMAT,0,sizeof(SPMAT));
      mem_numvar(TYPE_SPMAT,1);
   }
   /* fprintf(stderr,"Have SPMAT structure\n"); */
   
   A->row = rows = NEW_A(m,SPROW);
   if ( ! A->row )      /* can't allocate */
     error(E_MEM,"sp_get");
   else if (mem_info_is_on()) {
      mem_bytes(TYPE_SPMAT,0,m*sizeof(SPROW));
   }
   /* fprintf(stderr,"Have row structure array\n"); */
   
   A->start_row = NEW_A(n,int);
   A->start_idx = NEW_A(n,int);
   if ( ! A->start_row || ! A->start_idx )  /* can't allocate */
     error(E_MEM,"sp_get");
   else if (mem_info_is_on()) {
      mem_bytes(TYPE_SPMAT,0,2*n*sizeof(int));
   }
   for ( i = 0; i < n; i++ )
     A->start_row[i] = A->start_idx[i] = -1;
   /* fprintf(stderr,"Have start_row array\n"); */
   
   A->m = A->max_m = m;
   A->n = A->max_n = n;
   
   for ( i = 0; i < m; i++, rows++ )
   {
      rows->elt = NEW_A(maxlen,row_elt);
      if ( ! rows->elt )
    error(E_MEM,"sp_get");
      else if (mem_info_is_on()) {
     mem_bytes(TYPE_SPMAT,0,maxlen*sizeof(row_elt));
      }
      /* fprintf(stderr,"Have row %d element array\n",i); */
      rows->len = 0;
      rows->maxlen = maxlen;
      rows->diag = -1;
   }
   
   return A;
}
 
 
/* sp_free -- frees up the memory for a sparse matrix */
int sp_free(A)
SPMAT   *A;
{
   SPROW    *r;
   int  i;
   
   if ( ! A )
     return -1;
   if ( A->start_row != (int *)NULL ) {
      if (mem_info_is_on()) {
     mem_bytes(TYPE_SPMAT,A->max_n*sizeof(int),0);
      }
      free((char *)(A->start_row));
   }
   if ( A->start_idx != (int *)NULL ) {
      if (mem_info_is_on()) {
     mem_bytes(TYPE_SPMAT,A->max_n*sizeof(int),0);
      }
      
      free((char *)(A->start_idx));
   }
   if ( ! A->row )
   {
      if (mem_info_is_on()) {
     mem_bytes(TYPE_SPMAT,sizeof(SPMAT),0);
     mem_numvar(TYPE_SPMAT,-1);
      }
      
      free((char *)A);
      return 0;
   }
   for ( i = 0; i < A->m; i++ )
   {
      r = &(A->row[i]);
      if ( r->elt != (row_elt *)NULL ) {
     if (mem_info_is_on()) {
        mem_bytes(TYPE_SPMAT,A->row[i].maxlen*sizeof(row_elt),0);
     }
     free((char *)(r->elt));
      }
   }
   
   if (mem_info_is_on()) {
      if (A->row) 
    mem_bytes(TYPE_SPMAT,A->max_m*sizeof(SPROW),0);
      mem_bytes(TYPE_SPMAT,sizeof(SPMAT),0);
      mem_numvar(TYPE_SPMAT,-1);
   }
   
   free((char *)(A->row));
   free((char *)A);
 
   return 0;
}
 
 
/* sp_copy -- constructs a copy of a given matrix
   -- note that the max_len fields (etc) are no larger in the copy
   than necessary
   -- result is returned */
SPMAT   *sp_copy(A)
SPMAT   *A;
{
   SPMAT    *out;
   SPROW    *row1, *row2;
   int  i;
   
   if ( A == SMNULL )
     error(E_NULL,"sp_copy");
   if ( ! (out=NEW(SPMAT)) )
     error(E_MEM,"sp_copy");
   else if (mem_info_is_on()) {
      mem_bytes(TYPE_SPMAT,0,sizeof(SPMAT));
      mem_numvar(TYPE_SPMAT,1);
   }
   out->m = out->max_m = A->m;  out->n = out->max_n = A->n;
   
   /* set up rows */
   if ( ! (out->row=NEW_A(A->m,SPROW)) )
     error(E_MEM,"sp_copy");
   else if (mem_info_is_on()) {
      mem_bytes(TYPE_SPMAT,0,A->m*sizeof(SPROW));
   }
   for ( i = 0; i < A->m; i++ )
   {
      row1 = &(A->row[i]);
      row2 = &(out->row[i]);
      if ( ! (row2->elt=NEW_A(max(row1->len,3),row_elt)) )
    error(E_MEM,"sp_copy");
      else if (mem_info_is_on()) {
     mem_bytes(TYPE_SPMAT,0,max(row1->len,3)*sizeof(row_elt));
      }
      row2->len = row1->len;
      row2->maxlen = max(row1->len,3);
      row2->diag = row1->diag;
      MEM_COPY((char *)(row1->elt),(char *)(row2->elt),
           row1->len*sizeof(row_elt));
   }
   
   /* set up start arrays -- for column access */
   if ( ! (out->start_idx=NEW_A(A->n,int)) ||
       ! (out->start_row=NEW_A(A->n,int)) )
     error(E_MEM,"sp_copy");
   else if (mem_info_is_on()) {
      mem_bytes(TYPE_SPMAT,0,2*A->n*sizeof(int));
   }
   MEM_COPY((char *)(A->start_idx),(char *)(out->start_idx),
        A->n*sizeof(int));
   MEM_COPY((char *)(A->start_row),(char *)(out->start_row),
        A->n*sizeof(int));
   
   return out;
}
 
/* sp_col_access -- set column access path; i.e. nxt_row, nxt_idx fields
   -- returns A */
SPMAT   *sp_col_access(A)
SPMAT   *A;
{
   int  i, j, j_idx, len, m, n;
   SPROW    *row;
   row_elt  *r_elt;
   int  *start_row, *start_idx;
   
   if ( A == SMNULL )
     error(E_NULL,"sp_col_access");
   
   m = A->m;    n = A->n;
   
   /* initialise start_row and start_idx */
   start_row = A->start_row;    start_idx = A->start_idx;
   for ( j = 0; j < n; j++ )
   {    *start_row++ = -1;  *start_idx++ = -1;  }
   
   start_row = A->start_row;    start_idx = A->start_idx;
   
   /* now work UP the rows, setting nxt_row, nxt_idx fields */
   for ( i = m-1; i >= 0; i-- )
   {
      row = &(A->row[i]);
      r_elt = row->elt;
      len   = row->len;
      for ( j_idx = 0; j_idx < len; j_idx++, r_elt++ )
      {
     j = r_elt->col;
     r_elt->nxt_row = start_row[j];
     r_elt->nxt_idx = start_idx[j];
     start_row[j] = i;
     start_idx[j] = j_idx;
      }
   }
   
   A->flag_col = TRUE;
   return A;
}
 
/* sp_diag_access -- set diagonal access path(s) */
SPMAT   *sp_diag_access(A)
SPMAT   *A;
{
   int  i, m;
   SPROW    *row;
   
   if ( A == SMNULL )
     error(E_NULL,"sp_diag_access");
   
   m = A->m;
   
   row = A->row;
   for ( i = 0; i < m; i++, row++ )
     row->diag = sprow_idx(row,i);
   
   A->flag_diag = TRUE;
   
   return A;
}
 
/* sp_m2dense -- convert a sparse matrix to a dense one */
MAT *sp_m2dense(A,out)
SPMAT   *A;
MAT *out;
{
   int  i, j_idx;
   SPROW    *row;
   row_elt  *elt;
   
   if ( ! A )
     error(E_NULL,"sp_m2dense");
   if ( ! out || out->m < A->m || out->n < A->n )
     out = m_get(A->m,A->n);
   
   m_zero(out);
   for ( i = 0; i < A->m; i++ )
   {
      row = &(A->row[i]);
      elt = row->elt;
      for ( j_idx = 0; j_idx < row->len; j_idx++, elt++ )
    out->me[i][elt->col] = elt->val;
   }
   
   return out;
}
 
 
/*  C = A+B, can be in situ */
SPMAT *sp_add(A,B,C)
SPMAT *A, *B, *C;
{
   int i, in_situ;
   SPROW *rc;
   static SPROW *tmp;
 
   if ( ! A || ! B )
     error(E_NULL,"sp_add");
   if ( A->m != B->m || A->n != B->n )
     error(E_SIZES,"sp_add");
   if (C == A || C == B)
     in_situ = TRUE;
   else in_situ = FALSE;
 
   if ( ! C )
     C = sp_get(A->m,A->n,5);
   else {
      if ( C->m != A->m || C->n != A->n  )
    error(E_SIZES,"sp_add");
      if (!in_situ) sp_zero(C);
   }
 
   if (tmp == (SPROW *)NULL && in_situ) {
      tmp = sprow_get(MINROWLEN);
      MEM_STAT_REG(tmp,TYPE_SPROW);
   }
 
   if (in_situ)
     for (i=0; i < A->m; i++) {
    rc = &(C->row[i]);
    sprow_add(&(A->row[i]),&(B->row[i]),0,tmp,TYPE_SPROW);
    sprow_resize(rc,tmp->len,TYPE_SPMAT);
    MEM_COPY(tmp->elt,rc->elt,tmp->len*sizeof(row_elt));
    rc->len = tmp->len;
     }
   else
     for (i=0; i < A->m; i++) {
    sprow_add(&(A->row[i]),&(B->row[i]),0,&(C->row[i]),TYPE_SPMAT);
     }
 
   C->flag_col = C->flag_diag = FALSE;
 
   return C;
}
 
/*  C = A-B, cannot be in situ */
SPMAT *sp_sub(A,B,C)
SPMAT *A, *B, *C;
{
   int i, in_situ;
   SPROW *rc;
   static SPROW *tmp;
   
   if ( ! A || ! B )
     error(E_NULL,"sp_sub");
   if ( A->m != B->m || A->n != B->n )
     error(E_SIZES,"sp_sub");
   if (C == A || C == B)
     in_situ = TRUE;
   else in_situ = FALSE;
 
   if ( ! C )
     C = sp_get(A->m,A->n,5);
   else {
      if ( C->m != A->m || C->n != A->n  )
    error(E_SIZES,"sp_sub");
      if (!in_situ) sp_zero(C);
   }
 
   if (tmp == (SPROW *)NULL && in_situ) {
      tmp = sprow_get(MINROWLEN);
      MEM_STAT_REG(tmp,TYPE_SPROW);
   }
 
   if (in_situ)
     for (i=0; i < A->m; i++) {
    rc = &(C->row[i]);
    sprow_sub(&(A->row[i]),&(B->row[i]),0,tmp,TYPE_SPROW);
    sprow_resize(rc,tmp->len,TYPE_SPMAT);
    MEM_COPY(tmp->elt,rc->elt,tmp->len*sizeof(row_elt));
    rc->len = tmp->len;
     }
   else
     for (i=0; i < A->m; i++) {
    sprow_sub(&(A->row[i]),&(B->row[i]),0,&(C->row[i]),TYPE_SPMAT);
     }
 
   C->flag_col = C->flag_diag = FALSE;
   
   return C;
}
 
/*  C = A+alpha*B, cannot be in situ */
SPMAT *sp_mltadd(A,B,alpha,C)
SPMAT *A, *B, *C;
double alpha;
{
   int i, in_situ;
   SPROW *rc;
   static SPROW *tmp;
 
   if ( ! A || ! B )
     error(E_NULL,"sp_mltadd");
   if ( A->m != B->m || A->n != B->n )
     error(E_SIZES,"sp_mltadd");
   if (C == A || C == B)
     in_situ = TRUE;
   else in_situ = FALSE;
 
   if ( ! C )
     C = sp_get(A->m,A->n,5);
   else {
      if ( C->m != A->m || C->n != A->n  )
    error(E_SIZES,"sp_mltadd");
      if (!in_situ) sp_zero(C);
   }
 
   if (tmp == (SPROW *)NULL && in_situ) {
      tmp = sprow_get(MINROWLEN);
      MEM_STAT_REG(tmp,TYPE_SPROW);
   }
 
   if (in_situ)
     for (i=0; i < A->m; i++) {
    rc = &(C->row[i]);
    sprow_mltadd(&(A->row[i]),&(B->row[i]),alpha,0,tmp,TYPE_SPROW);
    sprow_resize(rc,tmp->len,TYPE_SPMAT);
    MEM_COPY(tmp->elt,rc->elt,tmp->len*sizeof(row_elt));
    rc->len = tmp->len;
     }
   else
     for (i=0; i < A->m; i++) {
    sprow_mltadd(&(A->row[i]),&(B->row[i]),alpha,0,
             &(C->row[i]),TYPE_SPMAT);
     }
   
   C->flag_col = C->flag_diag = FALSE;
   
   return C;
}
 
 
 
/*  B = alpha*A, can be in situ */
SPMAT *sp_smlt(A,alpha,B)
SPMAT *A, *B;
double alpha;
{
   int i;
 
   if ( ! A )
     error(E_NULL,"sp_smlt");
   if ( ! B )
     B = sp_get(A->m,A->n,5);
   else
     if ( A->m != B->m || A->n != B->n )
       error(E_SIZES,"sp_smlt");
 
   for (i=0; i < A->m; i++) {
      sprow_smlt(&(A->row[i]),alpha,0,&(B->row[i]),TYPE_SPMAT);
   }
   return B;
}
 
 
 
/* sp_zero -- zero all the (represented) elements of a sparse matrix */
SPMAT   *sp_zero(A)
SPMAT   *A;
{
   int  i, idx, len;
   row_elt  *elt;
   
   if ( ! A )
     error(E_NULL,"sp_zero");
   
   for ( i = 0; i < A->m; i++ )
   {
      elt = A->row[i].elt;
      len = A->row[i].len;
      for ( idx = 0; idx < len; idx++ )
    (*elt++).val = 0.0;
   }
   
   return A;
}
 
/* sp_copy2 -- copy sparse matrix (type 2) 
   -- keeps structure of the OUT matrix */
SPMAT   *sp_copy2(A,OUT)
SPMAT   *A, *OUT;
{
   int  i /* , idx, len1, len2 */;
   SPROW    *r1, *r2;
   static SPROW *scratch = (SPROW *)NULL;
   /* row_elt   *e1, *e2; */
   
   if ( ! A )
     error(E_NULL,"sp_copy2");
   if ( ! OUT )
     OUT = sp_get(A->m,A->n,10);
   if ( ! scratch ) {
      scratch = sprow_xpd(scratch,MINROWLEN,TYPE_SPROW);
      MEM_STAT_REG(scratch,TYPE_SPROW);
   }
 
   if ( OUT->m < A->m )
   {
      if (mem_info_is_on()) {
     mem_bytes(TYPE_SPMAT,A->max_m*sizeof(SPROW),
              A->m*sizeof(SPROW));
      }
 
      OUT->row = RENEW(OUT->row,A->m,SPROW);
      if ( ! OUT->row )
    error(E_MEM,"sp_copy2");
      
      for ( i = OUT->m; i < A->m; i++ )
      {
     OUT->row[i].elt = NEW_A(MINROWLEN,row_elt);
     if ( ! OUT->row[i].elt )
       error(E_MEM,"sp_copy2");
     else if (mem_info_is_on()) {
        mem_bytes(TYPE_SPMAT,0,MINROWLEN*sizeof(row_elt));
     }
     
     OUT->row[i].maxlen = MINROWLEN;
     OUT->row[i].len = 0;
      }
      OUT->m = A->m;
   }
   
   OUT->flag_col = OUT->flag_diag = FALSE;
   /* sp_zero(OUT); */
 
   for ( i = 0; i < A->m; i++ )
   {
      r1 = &(A->row[i]);    r2 = &(OUT->row[i]);
      sprow_copy(r1,r2,scratch,TYPE_SPROW);
      if ( r2->maxlen < scratch->len )
    sprow_xpd(r2,scratch->len,TYPE_SPMAT);
      MEM_COPY((char *)(scratch->elt),(char *)(r2->elt),
           scratch->len*sizeof(row_elt));
      r2->len = scratch->len;
      /*******************************************************
    e1 = r1->elt;       e2 = r2->elt;
    len1 = r1->len;     len2 = r2->len;
    for ( idx = 0; idx < len2; idx++, e2++ )
    e2->val = 0.0;
    for ( idx = 0; idx < len1; idx++, e1++ )
    sprow_set_val(r2,e1->col,e1->val);
    *******************************************************/
   }
 
   sp_col_access(OUT);
   return OUT;
}
 
/* sp_resize -- resize a sparse matrix
   -- don't destroying any contents if possible
   -- returns resized matrix */
SPMAT   *sp_resize(A,m,n)
SPMAT   *A;
int m, n;
{
   int  i, len;
   SPROW    *r;
   
   if (m < 0 || n < 0)
     error(E_NEG,"sp_resize");
 
   if ( ! A )
     return sp_get(m,n,10);
 
   if (m == A->m && n == A->n)
     return A;
 
   if ( m <= A->max_m )
   {
      for ( i = A->m; i < m; i++ )
    A->row[i].len = 0;
      A->m = m;
   }
   else
   {
      if (mem_info_is_on()) {
     mem_bytes(TYPE_SPMAT,A->max_m*sizeof(SPROW),
             m*sizeof(SPROW));
      }
 
      A->row = RENEW(A->row,(unsigned)m,SPROW);
      if ( ! A->row )
    error(E_MEM,"sp_resize");
      for ( i = A->m; i < m; i++ )
      {
     if ( ! (A->row[i].elt = NEW_A(MINROWLEN,row_elt)) )
       error(E_MEM,"sp_resize");
     else if (mem_info_is_on()) {
        mem_bytes(TYPE_SPMAT,0,MINROWLEN*sizeof(row_elt));
     }
     A->row[i].len = 0; A->row[i].maxlen = MINROWLEN;
      }
      A->m = A->max_m = m;
   }
 
   /* update number of rows */
   A->n = n;
 
   /* do we need to increase the size of start_idx[] and start_row[] ? */
   if ( n > A->max_n )
   {    /* only have to update the start_idx & start_row arrays */
      if (mem_info_is_on())
      {
      mem_bytes(TYPE_SPMAT,2*A->max_n*sizeof(int),
            2*n*sizeof(int));
      }
 
      A->start_row = RENEW(A->start_row,(unsigned)n,int);
      A->start_idx = RENEW(A->start_idx,(unsigned)n,int);
      if ( ! A->start_row || ! A->start_idx )
    error(E_MEM,"sp_resize");
      A->max_n = n; /* ...and update max_n */
 
      return A;
   }
 
   if ( n <= A->n )
       /* make sure that all rows are truncated just before column n */
       for ( i = 0; i < A->m; i++ )
       {
       r = &(A->row[i]);
       len = sprow_idx(r,n);
       if ( len < 0 )
           len = -(len+2);
       if ( len < 0 )
           error(E_MEM,"sp_resize");
       r->len = len;
       }
   
   return A;
}
 
 
/* sp_compact -- removes zeros and near-zeros from a sparse matrix */
SPMAT   *sp_compact(A,tol)
SPMAT   *A;
double  tol;
{
   int  i, idx1, idx2;
   SPROW    *r;
   row_elt  *elt1, *elt2;
   
   if (  ! A )
     error(E_NULL,"sp_compact");
   if ( tol < 0.0 )
     error(E_RANGE,"sp_compact");
   
   A->flag_col = A->flag_diag = FALSE;
   
   for ( i = 0; i < A->m; i++ )
   {
      r = &(A->row[i]);
      elt1 = elt2 = r->elt;
      idx1 = idx2 = 0;
      while ( idx1 < r->len )
      {
     /* printf("# sp_compact: idx1 = %d, idx2 = %d\n",idx1,idx2); */
     if ( fabs(elt1->val) <= tol )
     {  idx1++; elt1++; continue;   }
     if ( elt1 != elt2 )
       MEM_COPY(elt1,elt2,sizeof(row_elt));
     idx1++;    elt1++;
     idx2++;    elt2++;
      }
      r->len = idx2;
   }
   
   return A;
}
 
/* varying number of arguments */
 
#ifdef ANSI_C
 
/* To allocate memory to many arguments. 
   The function should be called:
   sp_get_vars(m,n,deg,&x,&y,&z,...,NULL);
   where 
     int m,n,deg;
     SPMAT *x, *y, *z,...;
     The last argument should be NULL ! 
     m x n is the dimension of matrices x,y,z,...
     returned value is equal to the number of allocated variables
*/
 
int sp_get_vars(int m,int n,int deg,...) 
{
   va_list ap;
   int i=0;
   SPMAT **par;
   
   va_start(ap, deg);
   while ((par = va_arg(ap,SPMAT **))) {   /* NULL ends the list*/
      *par = sp_get(m,n,deg);
      i++;
   } 
 
   va_end(ap);
   return i;
}
 
 
/* To resize memory for many arguments. 
   The function should be called:
   sp_resize_vars(m,n,&x,&y,&z,...,NULL);
   where 
     int m,n;
     SPMAT *x, *y, *z,...;
     The last argument should be NULL ! 
     m X n is the resized dimension of matrices x,y,z,...
     returned value is equal to the number of allocated variables.
     If one of x,y,z,.. arguments is NULL then memory is allocated to this 
     argument. 
*/
  
int sp_resize_vars(int m,int n,...) 
{
   va_list ap;
   int i=0;
   SPMAT **par;
   
   va_start(ap, n);
   while ((par = va_arg(ap,SPMAT **))) {   /* NULL ends the list*/
      *par = sp_resize(*par,m,n);
      i++;
   } 
 
   va_end(ap);
   return i;
}
 
/* To deallocate memory for many arguments. 
   The function should be called:
   sp_free_vars(&x,&y,&z,...,NULL);
   where 
     SPMAT *x, *y, *z,...;
     The last argument should be NULL ! 
     There must be at least one not NULL argument.
     returned value is equal to the number of allocated variables.
     Returned value of x,y,z,.. is VNULL.
*/
 
int sp_free_vars(SPMAT **va,...)
{
   va_list ap;
   int i=1;
   SPMAT **par;
   
   sp_free(*va);
   *va = (SPMAT *) NULL;
   va_start(ap, va);
   while ((par = va_arg(ap,SPMAT **))) {   /* NULL ends the list*/
      sp_free(*par); 
      *par = (SPMAT *)NULL;
      i++;
   } 
 
   va_end(ap);
   return i;
}
 
 
#elif VARARGS
 
/* To allocate memory to many arguments. 
   The function should be called:
   sp_get_vars(m,n,deg,&x,&y,&z,...,NULL);
   where 
     int m,n,deg;
     SPMAT *x, *y, *z,...;
     The last argument should be NULL ! 
     m x n is the dimension of matrices x,y,z,...
     returned value is equal to the number of allocated variables
*/
 
int sp_get_vars(va_alist) va_dcl
{
   va_list ap;
   int i=0, m, n, deg;
   SPMAT **par;
   
   va_start(ap);
   m = va_arg(ap,int);
   n = va_arg(ap,int);
   deg = va_arg(ap,int);
   while ((par = va_arg(ap,SPMAT **))) {   /* NULL ends the list*/
      *par = sp_get(m,n,deg);
      i++;
   } 
 
   va_end(ap);
   return i;
}
 
 
/* To resize memory for many arguments. 
   The function should be called:
   sp_resize_vars(m,n,&x,&y,&z,...,NULL);
   where 
     int m,n;
     SPMAT *x, *y, *z,...;
     The last argument should be NULL ! 
     m X n is the resized dimension of matrices x,y,z,...
     returned value is equal to the number of allocated variables.
     If one of x,y,z,.. arguments is NULL then memory is allocated to this 
     argument. 
*/
 
int sp_resize_vars(va_alist) va_dcl
{
   va_list ap;
   int i=0, m, n;
   SPMAT **par;
   
   va_start(ap);
   m = va_arg(ap,int);
   n = va_arg(ap,int);
   while ((par = va_arg(ap,SPMAT **))) {   /* NULL ends the list*/
      *par = sp_resize(*par,m,n);
      i++;
   } 
 
   va_end(ap);
   return i;
}
 
 
 
/* To deallocate memory for many arguments. 
   The function should be called:
   sp_free_vars(&x,&y,&z,...,NULL);
   where 
     SPMAT *x, *y, *z,...;
     The last argument should be NULL ! 
     There must be at least one not NULL argument.
     returned value is equal to the number of allocated variables.
     Returned value of x,y,z,.. is VNULL.
*/
 
int sp_free_vars(va_alist) va_dcl
{
   va_list ap;
   int i=0;
   SPMAT **par;
   
   va_start(ap);
   while ((par = va_arg(ap,SPMAT **))) {   /* NULL ends the list*/
      sp_free(*par); 
      *par = (SPMAT *)NULL;
      i++;
   } 
 
   va_end(ap);
   return i;
}
 
 
 
#endif