SparseMatrix.h

/*  _________________________________________________________________________
 *
 *  UTILIB: A utility library for developing portable C++ codes.
 *  Copyright (c) 2001, Sandia National Laboratories.
 *  This software is distributed under the GNU Lesser General Public License.
 *  For more information, see the README file in the top UTILIB directory.
 *  _________________________________________________________________________
 */

//
// SparseMatrix.h
//
#ifndef __SparseMatrix_h
#define __SparseMatrix_h

//--- INCLUDES ---//

#include "_math.h"
#include "BitArray.h"
#include "CharString.h"
#include "DoubleVector.h"
#include "IntVector.h"
#include "Ereal.h"
 
//--- CONSTANTS ---//

//--- TYPES ---//

template <class T>
class UTILIB_API CMSparseMatrix; // Forward decl

template <class T>
class UTILIB_API RMSparseMatrix; // Forward decl

//
// Note:  is this numerically stable?
//
template <class T>
inline void product(NumArray<T>& res, const RMSparseMatrix<T>& mat, 
                                                const NumArray<T>& vec)
{
res.resize(mat.get_nrows());

int nrows=mat.get_nrows();
for (int i=0; i<nrows; i++) {
  res[i] = 0.0;
  const T* val = &(mat.matval[mat.matbeg[i]]);
  const int*    col = &(mat.matind[mat.matbeg[i]]);
  for (int k=0; k<mat.matcnt[i]; k++, val++, col++) {
    res[i] += vec[*col] * (*val);
    }
  }
}


//---------------------------------------------------------------------------

template <class T>
class UTILIB_API SparseMatrix
{
public:

  int get_ncols( void ) const { return ncols; }
  int get_nrows( void ) const { return nrows; }
  int get_nnzero( void ) const {return nnzeros;}
  IntVector& get_matbeg( void ) { return matbeg; }
  IntVector& get_matcnt( void ) { return matcnt; }
  IntVector& get_matind( void ) { return matind; }
  NumArray<T>& get_matval( void ) { return matval; }

  virtual T& operator()(const int row, const int col) = 0;
  virtual const T& operator()(const int row, const int col) const = 0;

  virtual int pretty_print(ostream& ostr) const;
  virtual int write(ostream& ostr) const = 0;
  virtual int read(istream& istr) = 0;
  int write(PackBuffer& ostr) const;
  int read(UnPackBuffer& istr);

  virtual void delete_row(const int i) = 0;
  virtual void delete_col(const int i) = 0;
  virtual void delete_element(const int row, const int col) = 0;

protected:

  SparseMatrix( void ) :
        ncols(0), nrows(0), nnzeros(0) {}
  virtual ~SparseMatrix() {}

  virtual void initialize(int nrow, int ncol, int nnzero) = 0;
  void setup(int ncol, int nrow, int nnzero, int flag);
  int          ncols;
  int          nrows;
  int          nnzeros;
  IntVector    matbeg;
  IntVector    matcnt;
  IntVector    matind;
  NumArray<T> matval;

};


template <class T>
inline ostream& operator<<(ostream& output, const SparseMatrix<T>& array)
{ array.write(output); return(output); }

template <class T>
inline istream& operator>>(istream& input, SparseMatrix<T>& array)
{ array.read(input); return(input); }


template <class T>
inline UTILIB_API PackBuffer& operator<<(PackBuffer& output, const SparseMatrix<T
>& array)
{ array.write(output); return(output); }

template <class T>
inline UTILIB_API UnPackBuffer& operator>>(UnPackBuffer& input, SparseMatrix<T>&
array)
{ array.read(input); return(input); }

//---------------------------------------------------------------------------

template <class T>
class UTILIB_API CMSparseMatrix : public SparseMatrix<T>
{
public:

  friend RMSparseMatrix<T>;

  CMSparseMatrix( void ) {}
  CMSparseMatrix( int nrow, int ncol, int nnzero )
                {setup(nrow,ncol,nnzero,1);}
  CMSparseMatrix( CMSparseMatrix<Ereal<T> >& mat);
  ~CMSparseMatrix() {}

  void initialize(int nrow, int ncol, int nnzero)
                {
                setup(nrow,ncol,nnzero,1);
                ncols=ncol;
                nrows=nrow;
                nnzeros=nnzero;
                }
  
  T& operator()(const int row, const int col);
  const T& operator()(const int row, const int col) const;

  int adjoinColumn(int count, int* rowPosition, T* value);

  int write(ostream& ostr) const;
  int read(istream& istr);
  void convert(const RMSparseMatrix<T>& rowmajor);

  void delete_row(const int i);
  void delete_col(const int i);
  void delete_element(const int row, const int col);
};


//---------------------------------------------------------------------------

template <class T>
class UTILIB_API RMSparseMatrix : public SparseMatrix<T>
{
public:

  friend CMSparseMatrix<T>;

  friend UTILIB_API void product<T>(NumArray<T>& res, 
                                const RMSparseMatrix<T>& mat, 
                                const NumArray<T>& vec);

  RMSparseMatrix( void ) {}
  RMSparseMatrix( int nrow, int ncol, int nnzero )
                {setup(nrow,ncol,nnzero,0);}
  virtual ~RMSparseMatrix() {}

  void initialize(int nrow, int ncol, int nnzero)
                {
                setup(nrow,ncol,nnzero,0);
                ncols=ncol;
                nrows=nrow;
                nnzeros=nnzero;
                }
  void resize(const int rowndx, const int rowlen);

  T& operator()(const int row, const int col);
  const T& operator()(const int row, const int col) const;

  int adjoinRow(int count, int* colPosition, T* value);

  int write(ostream& ostr) const;
  int read(istream& istr);
  void convert(const CMSparseMatrix<T>& colmajor);

  void delete_row(const int i);
  void delete_col(const int i);
  void delete_element(const int row, const int col);
};




//
// This just sets up the memory, but does not modify the ncols, nrows or
// nnzeros values
//
template <class T>
void SparseMatrix<T>::setup(int nrow, int ncol, int nnzero, int flag)
{
matind.resize(nnzero);
matval.resize(nnzero);
if (flag) {
   matbeg.resize(ncol);
   matcnt.resize(ncol);
   }
else {
   matbeg.resize(nrow);
   matcnt.resize(nrow);
   }
}


template <class T>
CMSparseMatrix<T>::CMSparseMatrix( CMSparseMatrix<Ereal<T> >& mat)
{
initialize(mat.get_nrows(), mat.get_ncols(), mat.get_nnzero());
matind << mat.get_matind();
matval << mat.get_matval();
matbeg << mat.get_matbeg();
matcnt << mat.get_matcnt();
}


template <class T>
T& CMSparseMatrix<T>::operator()(const int row, const int col)
{
if ((row < 0) || (row >= nrows) || (col < 0) || (col >= ncols))
   ErrAbort(errmsg("CMSparseMatrix<T>::operator : iterator out of range. %dx%d not in %dx%d",row,col,nrows,ncols));

int pt = matbeg[col];
for (int i=0; i<matcnt[col]; i++)
  if (matind[pt+i] == row)
     return matval[pt+i];

return matval[0];
}


template <class T>
const T& CMSparseMatrix<T>::operator()(const int row, const int col) const
{
if ((row < 0) || (row >= nrows) || (col < 0) || (col >= ncols))
   ErrAbort(errmsg("CMSparseMatrix<T>::operator : iterator out of range. %dx%d not in %dx%d",row,col,nrows,ncols));

int pt = matbeg[col];
for (int i=0; i<matcnt[col]; i++)
  if (matind[pt+i] == row)
     return matval[pt+i];

return matval[0];
}


template <class T>
T& RMSparseMatrix<T>::operator()(const int row, const int col)
{
if ((row < 0) || (row >= nrows) || (col < 0) || (col >= ncols))
   ErrAbort(errmsg("RMSparseMatrix<T>::operator : iterator out of range. %dx%d not in %dx%d",row,col,nrows,ncols));

int pt = matbeg[row];
for (int i=0; i<matcnt[row]; i++)
  if (matind[pt+i] == col)
     return matval[pt+i];

return matval[0];
}


template <class T>
const T& RMSparseMatrix<T>::operator()(const int row, const int col) const
{
if ((row < 0) || (row >= nrows) || (col < 0) || (col >= ncols))
   ErrAbort(errmsg("RMSparseMatrix<T>::operator : iterator out of range. %dx%d not in %dx%d",row,col,nrows,ncols));

int pt = matbeg[row];
for (int i=0; i<matcnt[row]; i++)
  if (matind[pt+i] == col)
     return matval[pt+i];

return matval[0];
}


template <class T>
int SparseMatrix<T>::pretty_print(ostream& ostr) const
{
for (int i=0; i<nrows; i++) {
  for (int j=0; j<ncols; j++)
    ostr << operator()(i,j) << " ";
  ostr << endl;
  }
return OK;
}


template <class T>
int CMSparseMatrix<T>::write(ostream& ostr) const 
{
ostr << nrows << " " << ncols << " " << nnzeros << ":" << endl;
int ndx=0;
for (int i=0; i<ncols; i++)
  for (int j=0; j<matcnt[i]; j++) {
    ostr << " (" << matind[ndx] << "," << i << ") " << matval[ndx] << endl;
    ndx++;
    }

return OK;
}


template <class T>
int CMSparseMatrix<T>::read(istream& istr)
{
ErrAbort("CMSparseMatrix<T>::read - not implemented!");
return ERR;
}


template <class T>
int SparseMatrix<T>::write(PackBuffer& ostr) const 
{
ostr << nrows << ncols << nnzeros;
ostr << matbeg << matcnt << matind << matval;

return OK;
}


template <class T>
int SparseMatrix<T>::read(UnPackBuffer& istr)
{
int ntmp, ctmp, ztmp;
istr >> ntmp >> ctmp >> ztmp;
initialize(ntmp,ctmp,ztmp);
istr >> matbeg >> matcnt >> matind >> matval;
return OK;
}


template <class T>
int RMSparseMatrix<T>::write(ostream& ostr) const 
{
ostr << nrows << " " << ncols << " " << nnzeros << ":" << endl;
int ndx=0;
for (int i=0; i<nrows; i++)
  for (int j=0; j<matcnt[i]; j++) {
    ostr << " (" << i << "," << matind[ndx] << ") " << matval[ndx] << endl;
    ndx++;
    }

return OK;
}


template <class T>
int RMSparseMatrix<T>::read(istream& istr)
{
ErrAbort("RMSparseMatrix<T>::read - not implemented!");
return ERR;
}



template <class T>
void RMSparseMatrix<T>::delete_row(const int i)
{
//
// Delete all of the elements in the row
//
while (matcnt[i] > 0)
  delete_element(i,matind[matbeg[i]]);
//
// Delete the row itself
//
for (int j=i; j<(nrows-1); j++) {
  matbeg[j] = matbeg[j+1];
  matcnt[j] = matcnt[j+1];
  }
nrows--;
}


template <class T>
void RMSparseMatrix<T>::delete_col(const int i)
{
//
// Delete all of the elements in the column
//
for (int j=0; j<nrows; j++) {
  //
  // Look for column i in row j
  //
  int k=0;
  while (k < matcnt[j]) {
    if (matind[matbeg[j]+k] == i)
       delete_element(j,i);
    else {
       if (matind[matbeg[j]+k] > i)
          matind[matbeg[j]+k]--;
       k++;
       }
    }
  }
//
// Delete the column itself (the column indeces are decremented above!)
//
ncols--;
}


template <class T>
void RMSparseMatrix<T>::delete_element(const int rowndx, const int colndx)
{
int k=0;
for (; k<matcnt[rowndx]; k++)
  if (matind[matbeg[rowndx]+k] == colndx) {
     break;
     }
if (k == matcnt[rowndx]) return;

for (int i=matbeg[rowndx]+k; i<(nnzeros-1); i++) {
  matind[i] = matind[i+1];
  matval[i] = matval[i+i];
  }
matcnt[rowndx]--;
for (int j=rowndx+1; j<nrows; j++)
  matbeg[j]--;
nnzeros--;
}


template <class T>
void RMSparseMatrix<T>::resize(const int rowndx, const int rowlen)
{
if (rowlen == matcnt[rowndx])
   return;
else if (rowlen > matcnt[rowndx]) {
   // TODO - resize the array up!
   }
else {
   int diff = rowlen-matcnt[rowndx];

   if (diff < 0) {
      for (int i=matbeg[rowndx]+rowlen; i<(nnzeros+diff); i++) {
        matind[i] = matind[i-diff];
        matval[i] = matval[i-diff];
        }
      }
   matcnt[rowndx] += diff;
   for (int j=rowndx+1; j<nrows; j++)
     matbeg[j] += diff;
   nnzeros += diff;
   }
}


template <class T>
void CMSparseMatrix<T>::convert(const RMSparseMatrix<T>& rowmajor)
{
matval.resize(rowmajor.nnzeros);
matind.resize(rowmajor.nnzeros);

ncols = rowmajor.ncols;
nrows = rowmajor.nrows;
nnzeros = rowmajor.nnzeros;
matcnt.resize(ncols);
matcnt=0;
matbeg.resize(ncols);

IntVector rowndx(nrows);
rowndx=0;
int prev_cndx=0;

for (int i=0; i<nnzeros; i++) {
  //
  // Set dummy values for cndx and rndx
  //
  int cndx=ncols;
  int rndx=-1;
  //
  // Find the point in the matrix with the next highest column index
  //
  for (int j=0; j<nrows; j++) {
    if ((rowmajor.matcnt[j] > 0) && (rowmajor.matcnt[j] > rowndx[j])) {
       if (rowmajor.matind[rowmajor.matbeg[j]+rowndx[j]] < cndx) {
          cndx = rowmajor.matind[rowmajor.matbeg[j]+rowndx[j]];
          rndx = j;
          }
       }
    if (cndx == prev_cndx)      // Can stop early because we're in the same col
       break;
    }
  //
  // Add this new point to the CM matrix
  //
  if (matcnt[cndx] == 0)
     matbeg[cndx] = i;
  matcnt[cndx]++;
  matind[i] = rndx;
  matval[i] = rowmajor.matval[rowmajor.matbeg[rndx]+rowndx[rndx]];

  rowndx[rndx]++;
  prev_cndx=cndx;
  }
}


template <class T>
void RMSparseMatrix<T>::convert(const CMSparseMatrix<T>& colmajor)
{
matval.resize(colmajor.nnzeros);
matind.resize(colmajor.nnzeros);

ncols = colmajor.ncols;
nrows = colmajor.nrows;
nnzeros = colmajor.nnzeros;
matcnt.resize(nrows);
matcnt=0;
matbeg.resize(nrows);

IntVector colndx(ncols);
colndx=0;
int prev_rndx=0;

for (int i=0; i<nnzeros; i++) {
  //
  // Find the point in the matrix with the next highest row index
  //
  int cndx=-1;
  int rndx=nrows;
  for (int j=0; j<ncols; j++) {
    if ((colmajor.matcnt[j] > 0) && (colmajor.matcnt[j] > colndx[j])) {
       if (colmajor.matind[colmajor.matbeg[j]+colndx[j]] < rndx) {
          cndx = j;
          rndx = colmajor.matind[colmajor.matbeg[j]+colndx[j]];
          }
       }
    if (prev_rndx == rndx)      // Can stop because we are on the same row
       break;
    }
  
  //
  // Add this new point to the RM matrix
  //
  if (matcnt[rndx] == 0)
     matbeg[rndx] = i;
  matcnt[rndx]++;
  matind[i] = cndx;
  matval[i] = colmajor.matval[colmajor.matbeg[cndx]+colndx[cndx]];

  colndx[cndx]++;
  prev_rndx=rndx;
  }
}



template <class T>
int CMSparseMatrix<T>::adjoinColumn(int count, int* rowPosition, T* value)
{
if (ncols == (int) matcnt.size()) {
   matcnt.resize(ncols+32);
   matbeg.resize(ncols+32);
   }

if (nnzeros + count > (int) matval.size()) {
   matval.resize(nnzeros+max(count,64));
   matind.resize(nnzeros+max(count,64));
   }

matbeg[ncols] = nnzeros;
matcnt[ncols++] = count;
for (int i=0; i<count; i++) {
  matind[nnzeros]   = rowPosition[i];
  // If we're building from scratch, number of rows is the highest row
  // number seen so far (+1 since 0 based)
  if (rowPosition[i] >= nrows)
    nrows = rowPosition[i] + 1;
  matval[nnzeros++] = value[i];
  }

return nnzeros;
}


template <class T>
void CMSparseMatrix<T>::delete_col(const int i)
{
//
// Delete all of the elements in the col 
//
while (matcnt[i])
  delete_element(matind[matbeg[i]],i);
//
// Delete the col itself
//
for (int j=i; j<(ncols-1); j++) {
  matbeg[j] = matbeg[j+1];
  matcnt[j] = matcnt[j+1];
  }
ncols--;
}


template <class T>
void CMSparseMatrix<T>::delete_row(const int i)
{
//
// Delete all of the elements in the column
//
for (int j=0; j<ncols; j++) {
  //
  // Look for row i in row j
  //
  int k=0;
  while (k < matcnt[j]) {
    if (matind[matbeg[j]+k] == i)
       delete_element(i,j);
    else {
       if (matind[matbeg[j]+k] > i)
          matind[matbeg[j]+k]--;
       k++;
       }
    }
  }
//
// Delete the row itself (the row indeces are decremented above!)
//
nrows--;
}


template <class T>
void CMSparseMatrix<T>::delete_element(const int rowndx, const int colndx)
{
int k=0;
for (; k<matcnt[colndx]; k++)
  if (matind[matbeg[colndx]+k] == rowndx) {
     break;
     }
if (k == matcnt[colndx]) return;

for (int i=matbeg[colndx]+k; i<(nnzeros-1); i++) {
  matind[i] = matind[i+1];
  matval[i] = matval[i+1];
  }
matcnt[colndx]--;
for (int j=colndx+1; j<ncols; j++)
  matbeg[j]--;
nnzeros--;
}



template <class T>
int RMSparseMatrix<T>::adjoinRow(int count, int* colPosition, T* value)
{
if (nrows == (int) matcnt.size()) {
   matcnt.resize(nrows+32);
   matbeg.resize(nrows+32);
   }

if (nnzeros + count > (int) matval.size()) {
   matval.resize(nnzeros+max(count,64));
   matind.resize(nnzeros+max(count,64));
   }

matbeg[nrows] = nnzeros;
matcnt[nrows++] = count;
for (int i=0; i<count; i++) {
  matind[nnzeros]   = colPosition[i];
  matval[nnzeros++] = value[i];
  }

return nnzeros;
}



#endif