/*!
  \file    sparsematrix.cc
  \brief   sparse matrix in compressed row storage
  \date    30 April 2010
  \author  Clemens Pechstein

  Institute of Computational Mathematics
  Johannes Kepler University Linz, Austria

  Tutorials -- Numerical methods for elliptic PDEs
  summer semester 2010
*/
#include "sparsematrix.hh"

using namespace std;

////////////////////////////////////////////////////
//      S p a r s e   S h a p e     routines

SparseShape :: SparseShape ()
  : rows_(0), cols_(0), elems_(0), num_(), pos_()
{
}

  
SparseShape :: SparseShape (int rows, int cols) :
  rows_(rows), cols_(cols), elems_(0),
  num_(rows_), pos_(rows_)
{
  for (int i=0; i<rows_; ++i) { num_[i] = 0; pos_[i].clear (); }
}

  
SparseShape :: SparseShape (const SparseShape& ss)
{
  rows_ = ss.rows_;
  cols_ = ss.cols_;
  elems_ = ss.elems_;
  num_ = ss.num_;
  pos_ = ss.pos_;
}


SparseShape& SparseShape :: operator= (const SparseShape& ss)
{
  if (&ss != this)
    {
      rows_ = ss.rows_;
      cols_ = ss.cols_;
      elems_ = ss.elems_;
      num_ = ss.num_;
      pos_ = ss.pos_;
    }
  return (*this);
}


SparseShape :: ~SparseShape ()
{
}


void SparseShape :: resize (int rows, int cols)
{
  rows_ = rows;
  cols_ = cols;
  elems_ = 0;
  num_.resize(rows_);
  pos_.resize(rows_);
  for (int i=0; i<rows_; ++i) { num_[i] = 0; pos_[i].clear (); }
}

void SparseShape :: addPosition (int row, int col)
{
  assert (0 <= row && row < rows_);
  pos_[row].insert (col);
  if (num_[row] < (int)pos_[row].size ()) ++elems_;
  num_[row] = pos_[row].size ();
}


bool SparseShape :: isNonZero (int row, int col) const
{
  assert (0 <= row && row < rows_);
  return (pos_[row].count (col) == 1);
}


void SparseShape :: print (ostream& os) const
{
  for (int i=0; i<rows_; i++)
    {
      string str(cols_, ' ');
      set<int>::const_iterator it = pos_[i].begin ();
      while (it != pos_[i].end ())
	{
	  str[*it] = '*';
	  ++it;
	}
      os << '[' << str << ']' << endl;
    }
}

void SparseShape :: printShort (ostream& os) const
{
  for (int i=0; i<rows_; i++)
    {
      os << "row " << i+1 << " : [ ";
      set<int>::const_iterator it = pos_[i].begin ();
      while (it != pos_[i].end ())
	{
	  os << *it << ' ';
	  ++it;
	}
      os << ']' << endl;
    }
}

/////////////////////////////////////////////
//    S p a r s e M a t r i x     routines

void SparseMatrix :: getP (int row, int col, int& p, bool& found) const
{
  assert (0 <= row && row < rows_);
  assert (0 <= col && col < cols_);
  
  p = _start[row];
  
  while (p < elems_ && p < _start[row+1] && pos_[p] < col) ++p;
  while (p >= 0 && p >= _start[row] && pos_[p] > col) --p;
  
  found = (0 <= p && p < elems_ && _start[row] <= p && p < _start[row+1] && pos_[p] == col);
}

SparseMatrix :: SparseMatrix ()
  : rows_(0), cols_(0), elems_(0), _start(), pos_(), value_()
{
  ;
}
  
SparseMatrix :: SparseMatrix (const SparseShape& ss)
  : rows_(ss.rows_), cols_(ss.cols_), elems_(ss.elems_), _start(), pos_(), value_()
{
  resize (ss); 
}
  
SparseMatrix :: SparseMatrix (const SparseMatrix& A)
  : rows_(A.rows_), cols_(A.cols_), elems_(A.elems_),
    _start(A._start), pos_(A.pos_), value_(A.value_)
{ 
  ; 
}
  
SparseMatrix& SparseMatrix :: operator= (const SparseMatrix& A)
{
  if (this != &A)
    {
      rows_ = A.rows_;
      cols_ = A.cols_;
      elems_ = A.elems_;
      _start = A._start;
      pos_ = A.pos_;
      value_ = A.value_;
    }
  return (*this);
}

SparseMatrix :: ~SparseMatrix ()
{ 
  ; 
}

void SparseMatrix :: resize (const SparseShape& ss)
{
  rows_ = ss.rows_;
  cols_ = ss.cols_;
  elems_ = ss.elems_;
  
  _start.resize(rows_+1);
  pos_.reserve(elems_);
  value_.resize(elems_);
  
  _start[0] = 0;
  for (int i=0; i<rows_; ++i)
    {
      std::set<int>::const_iterator it;
      it = ss.pos_[i].begin ();
      while (it != ss.pos_[i].end ())
	{
	  pos_.push_back ((int)(*it));
	  ++it;
	}
      _start[i+1] = _start[i] + ss.num_[i];
    }
  assert ((int)pos_.size () == elems_);
}

void SparseMatrix :: getColIndices (int row, ::set<int>& indices) const
{
  assert (0 <= row && row < rows_);
  indices.clear ();
  int p = _start[row];
  while (p < _start[row+1])
    {
      indices.insert (pos_[p]);
      ++p;
    }
}

bool SparseMatrix :: isZero (int row, int col) const
{
  int p;
  bool found;
  getP (row, col, p, found);
  return ! found;
}

const double SparseMatrix :: operator() (int row, int col) const
{
  int p;
  bool found;
  getP (row, col, p, found);
  if (found) return value_[p]; else return 0;
}

double& SparseMatrix :: operator() (int row, int col)
{
  int p;
  bool found;
  getP (row, col, p, found);
  
  if (found)
    return value_[p];
  else
    assert (false && "FAILED to return reference to zero element in SparseMatrix !");
}

SparseMatrix& SparseMatrix :: operator= (double value)
{
  for (int i=0; i<elems_; ++i) value_[i] = value;
  return *this;
}

void SparseMatrix :: apply (const Vector& v, Vector& w) const
{
  assert ((int)v.size () == cols_);
  
  int i, p;
  w.resize (rows_, 0); // resize and initialize by zero
  
  p = 0;
  for (i=0; i<rows_; ++i)
    {
      while (p < _start[i+1])
	{
	  w[i] += value_[p] * v[pos_[p]];
	  ++p;
	}
    }
}

void SparseMatrix :: applyTrans (const Vector& v, Vector& w) const
{
  assert ((int)v.size () == rows_);
  
  int i, p;
  w.resize (cols_, 0); // resize and initialize by zero
  
  p = 0;
  for (i=0; i<rows_; ++i)
    {
      while (p < _start[i+1])
	{
	  w[pos_[p]] += value_[p] * v[i];
	  ++p;
	}
    }
}

int SparseMatrix :: widthInRow (int row) const
{
  assert (0 <= row && row < rows_);
  int i = _start[row];
  int min = pos_[i];
  int max = pos_[i];
  ++i;
  while (i < _start[row+1])
    {
      if (pos_[i] < min) { min = pos_[i]; }
      if (pos_[i] > max) { max = pos_[i]; }
      ++i;
    }
    
  return (max-min);
}

int SparseMatrix :: maxWidth () const
{
  int mw = 0;
  if (rows_ >= 0)
    {
      mw = widthInRow (0);
      int w;
      for (int i=1; i<rows_; i++)
	{
	  w = widthInRow (i);
	  if (w > mw) mw = w;
	}
    }
  return mw;
}

void SparseMatrix :: print (ostream& os, int precision, int rS, int rE, int cS, int cE, bool matlabformat) const
{
  if (rE < 0) rE = rows_-1;
  if (cE < 0) cE = cols_-1;
  assert (0 <= rS && rE < rows_);
  assert (0 <= cS && cE < cols_);

  int i, j;
  int lj;
    
  os.precision (precision);
    
  if (matlabformat) os << "[";

  for (i=rS; i<=rE; i++)
    {
      os << "[ ";
	
      j = 0;
      lj = _start[i];
      while (lj < _start[i+1])
	{
	  while (j < pos_[lj])
	    {
	      if (cS <= j && j <= cE)
		if (matlabformat)
		  {
		    os << "0"; if (j != cE) os << ", "; 
		  }
		else
		  {
		    os << "0 ";
		  }
		
	      j++;
	    }
	    
	  if (cS <= j && j <= cE)
	    {
	      os << value_[lj];
	      if (matlabformat) { if (j != cE) os << ", "; else os << ' '; }
	      else os << ' ';
	    }
	  lj++;
	  j++;
	}
      while (j <= cE)
	{
	  os << "0";
	  if (matlabformat) { if (j != cE) os << ", "; else os << ' '; }
	  else os << ' ';
	    
	  j++;
	}
	
      if (matlabformat) 
	{
	  if (i != rE) os << "], ";
	  else os << "]";
	}
      else
	os << "]" << endl;;
    }
    
  if (matlabformat) os << "]" << endl;
}

void SparseMatrix :: printShort (ostream& os, int precision) const
{
  int i, j;
    
  cout<<"SparseMatrix["<<elems_<<"]"<<endl;
  for (i=0; i<rows_; i++)
    {
      os << "row " << i << ": [";
      j = _start[i];
      while (j<_start[i+1])
	{
	  os << '(' << pos_[j] << ": ";
	  os.precision (precision);
	  os << value_[j] << ") ";
	  j++;
	}
      os << ']' << endl;
    }
}

// EOF