/**************************************************************************
**       Title: crs block matrix class
**        Date: 8.5.2006
**   Copyright: GPL Bernard Haasdonk
** Description: crs block matrix class
**************************************************************************/

#ifndef CRSBLOCK_MATRIX_HH
#define CRSBLOCK_MATRIX_HH

#include <iostream>

using namespace std;


template <class BlockType>
class CRSBlockMatrix
{
public:  
  CRSBlockMatrix(const int nrows, const int maxblocks)
        {
          nonZero_ = maxblocks;
          blocks_ = new BlockType[nrows*maxblocks];
          cols_ = new int[nrows*maxblocks];
          nrows_= nrows;
          for (int i=0;i<nrows*maxblocks;i++)
          {
            cols_[i] = -1;
            blocks_[i].setAll(0.0);
          }
          
        }
  
  ~CRSBlockMatrix()
        {
          delete[] cols_;
          delete[] blocks_;
        }
    
  int nrows() const
        {
          return nrows_;
        }
  
  void setBlock(const int row,const int col, const BlockType block)
        {
          int ci = 0;
//          cout << "entered setBlock \n";
          
          while (cols_[row*nonZero_+ci]!=col && 
                 cols_[row*nonZero_+ci]!=-1)
          {
            ci++;
            
            if (ci==nonZero_)
            {
              std::cerr << "error: too many block in row.\n";
              abort();
            }
          }
          
//          cout << "block before setting: \n" << 
//              blocks_[row*nonZero_+ci];
          blocks_[row*nonZero_+ci] = block;
          if (cols_[row*nonZero_+ci]==-1)
              cols_[row*nonZero_+ci]= col;            
//          cout << "block after setting: \n" << 
//              blocks_[row*nonZero_+ci];        
        } 
  
  int numEntries() const
        {
          return nrows_* nonZero_ * blocks_[0].nrows()* blocks_[0].ncols();
        }
  
  void mult(const double *v, double *w) const
        {
          int rows_per_block = blocks_[0].nrows();
          int cols_per_block = blocks_[0].ncols();
          double* wbegin = w;
          double* wend = w + rows_per_block;
          double temp[rows_per_block];
          
          for (double* wp = w; wp!= w+rows_per_block*nrows_;wp++)
              *wp = 0.0;
          
          for (int i=0;i<nrows_;i++)
          {
            for (int ci=0; ci<nonZero_; ci++)
                if (cols_[i*nonZero_+ci]!=-1)
                {
                  blocks_[i*nonZero_+ci].mult(
                      v+cols_[i*nonZero_+ci]*cols_per_block,
                      temp);
                  for (int j=0;j<rows_per_block;j++)
                      w[i*rows_per_block+j] += temp[j];                  
                }
          } 
        }
  
  template <class BT>
  friend std::ostream& operator<<(std::ostream &os,
                                   const CRSBlockMatrix<BT> &mat);  

private:
  
  int nrows_, nonZero_;
  int * cols_;
  BlockType* blocks_;
};

template <class BlockType>
std::ostream& operator <<(std::ostream &os,
                          const CRSBlockMatrix<BlockType> &mat)
{
  os << "number of block rows: "<< mat.nrows() << " \n";
  os << "number of nonzero blocks per row: "<< mat.nonZero_ << " \n";
  
  os << "[\n";
  for (int i=0;i<mat.nrows();i++)
  {
    int ci=0;
//    while (mat.cols_[ci]!=-1)
    for (ci=0;ci<mat.nonZero_;ci++)
    {
      os << "Block-Matrix row = " << i << 
          ", column(" << ci << ") = " << mat.cols_[i*mat.nonZero_+ci] << 
          ": \n";
      os << mat.blocks_[i*mat.nonZero_+ci];
//      ci++;
    }
  }
  os << "]\n";
  return os;
}


#endif