#ifndef __SPARSEMATRIX_H
#define __SPARSEMATRIX_H
/*!
  \file    sparsematrix.hh
  \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 <vector>
#include <valarray>
#include <set>
#include <algorithm>
#include <cassert>
#include <iostream>
#include "vector.hh"

using namespace std;

/// class for the non-zero pattern of the matrix
class SparseShape
{  
public:
  /// default constructor
  SparseShape ();
  
  /// constructor 
  SparseShape (int rows, int cols);
  
  /// copy constructor
  SparseShape (const SparseShape& ss);
  
  /// assignment
  SparseShape& operator= (const SparseShape& ss);

  /// destructor
  ~SparseShape ();
  
  /// number of rows
  int rows () const { return rows_; };

  /// number of columns
  int cols () const { return cols_; };
  
  /// number of non-zero elements
  int elems () const { return elems_; };
  
  /// number of (non-zero) entries in row
  int elemsInRow (int row) const 
  {
    assert (0 <= row && row < rows_);
    return num_[row]; 
  }
  
  /// clear current pattern and resize
  void resize (int rows, int cols);
  
  // add non-zero entry to pattern
  void addPosition (int row, int col);
  
  /// is this a non-zero entry?
  bool isNonZero (int row, int col) const;
  
  /// pretty output of the pattern
  void print (ostream& os) const;
  
  /// compressed output of the pattern
  void printShort (ostream& os) const;
  
private:
  
  int rows_;    ///< number of rows
  int cols_;    ///< number of columns
  int elems_;   ///< total number of elements
  
  vector<int> num_;         ///< number of elements per row
  vector<set<int> >pos_;    ///< position array (indices counted from 0 on)
  
  friend class SparseMatrix;
  
}; // class SparseShape



/// sparse matrix in compressed row storage
class SparseMatrix
{
public:
  
  /// default constructor
  SparseMatrix ();

  /// create matrix from SparseShape
  SparseMatrix (const SparseShape& ss);

  /// copy constructor
  SparseMatrix (const SparseMatrix& A);

  /// assignment
  SparseMatrix& operator= (const SparseMatrix& A);

  /// destructor
  ~SparseMatrix ();
  
  /// number of rows
  int rows () const { return rows_; }

  /// number of columns
  int cols () const { return cols_; }
  
  /// number of non-zero telements
  int elems () const { return elems_; }
  
  /// clear matrix and rebuild according to ss
  void resize (const SparseShape& ss);
  
  /// return all column indices of row as a set
  void getColIndices (int row, set<int>& indices) const;

  /// set all matrix elements to zero
  
  /// is this an entry off the non-zero pattern?
  bool isZero (int row, int col) const;
  
  /// const element access
  const double operator() (int row, int col) const;
  
  /// element access
  double& operator() (int row, int col);
  
  /// set all elements to one value
  SparseMatrix& operator= (double value);
  
  /// multiply a vector v by this matrix and store the result in w
  void apply (const Vector& v, Vector& w) const;
  
  /// multiply a vector v by the transpose of this matrix and store the result in w
  void applyTrans (const Vector& v, Vector& w) const;
  
  /// matrix-vector multiplication
  Vector operator* (const Vector& v) const
  {
    Vector w(rows_);
    apply(v, w);
    return w; 
  }
  
  /// band width of row
  int widthInRow (int row) const;
  
  /// band width of this matrix
  int maxWidth () const;
  
  /// pretty-print this matrix
  void print (ostream& os, bool matlabformat=false) const
  { print (os, 14, 0, rows_-1, 0, cols_-1, matlabformat); }
  
  /// pretty-print sub-block of this matrix
  void print (ostream& os, int precision, int rS=0, int rE=-1, int cS=0, int cE=-1,
	      bool matlabformat=false) const;
  
  /// print this matrix in compressed form
  void printShort (ostream& os, int precision=14) const;
  
private:
  
  int rows_;   ///< number or rows
  int cols_;   ///< number of columns
  int elems_;  ///< total number of elements
  
  vector<int> _start;      ///< index in pos_/value_ where the .-th row starts
  vector<int> pos_;        ///< column indices
  valarray<double> value_; ///< values of the entries
  
  /// get index p in pos_ / value_ for (row, col); found indicates if this is really a matrix entry
  void getP (int row, int col, int& p, bool& found) const;
  
}; // class SparseMatrix



#endif // __SPARSEMATRIX_H