#ifndef __MAT_H
#define __MAT_H
/*!
  \file    mat.hh
  \brief   static matrix
  \date    15 April 2010
  \author  Clemens Pechstein

  Institute of Computational Mathematics
  Johannes Kepler University Linz, Austria

  Tutorials -- Numerical methods for elliptic PDEs
  summer semester 2010
*/

#include <iostream>
#include <cassert>
#include "vec.hh"

template <int Rows, int Cols>
class Mat
{
public:
  /// default constructor
  Mat ()
  {
    for (int k=0; k<Rows*Cols; ++k) data_[k] = 0;
  }
  
  /// destructor
  ~Mat ()
  {}
  
  /// copy constructor
  Mat (const Mat& m)
  {
    for (int k=0; k<Rows*Cols; ++k) data_[k] = m.data_[k];
  }
  
  /// assignment
  Mat& operator= (const Mat& m)
  {
    for (int k=0; k<Rows*Cols; ++k) data_[k] = m.data_[k];
    return *this;
  }
  
  /// assign all entries to val
  Mat& operator= (const double val)
  {
    for (int k=0; k<Rows*Cols; ++k) data_[k] = val;
    return *this;
  }
  
  /// element access
  const double operator() (int i, int j) const
  {
    assert (0 <= i && i < Rows && 0 <= j && j < Cols && "range check");
    return data_[Cols*i + j];
  }
  
  /// element access by reference
  double& operator() (int i, int j)
  {
    assert (0 <= i && i < Rows && 0 <= j && j < Cols && "range check");
    return data_[Cols*i + j];
  }
  
  /// matrix vector multiplication
  void apply (const Vec<Cols>& v, Vec<Rows>& res) const
  {
    assert (&v != &res && "check to avoid side effect");
    for (int i=0; i<Rows; ++i)
      {
	res[i] = 0;
	for (int j=0; j<Cols; ++j)
	  {
	    res[i] += data_[Cols*i + j] * v[j];
	  }
      }
  }

  /// matrix vector multiplication
  Vec<Rows> operator* (const Vec<Cols>& v) const
  {
    Vec<Rows> res;
    apply (v, res);
    return res;
  }
  
  /// transposed matrix vector multiplication
  void applyTrans (const Vec<Rows>& v, Vec<Cols>& res) const
  {
    assert (&v != &res && "check to avoid side effect");
    res = 0;
    for (int i=0; i<Rows; ++i)
      {
	for (int j=0; j<Cols; ++j)
	  {
	    res[j] += data_[Cols*i + j] * v[i];
	  }
      }
  }
  
  /// increment
  Mat& operator+= (const Mat& m)
  {
    for (int k=0; k<Rows*Cols; ++k) data_[k] += m.data_[k];
    return *this;
  }
  
  /// decrement
  Mat& operator-= (const Mat& m)
  {
    for (int k=0; k<Rows*Cols; ++k) data_[k] += m.data_[k];
    return *this;
  }
  
  /// factor
  Mat& operator*= (const double factor)
  {
    for (int k=0; k<Rows*Cols; ++k) data_[k] *= factor;
    return *this;
  }
  
  /// output
  void print (std::ostream& os) const
  {
    os << std::endl;
    for (int i=0; i<Rows; ++i)
      {
	os << "\t[ " << data_[i*Cols + 0];
	for (int j=1; j<Cols; ++j) os << '\t' << data_[i*Cols + j];
	os << " ]" << std::endl;
      }
  }

  /// transpose the matrix (works only for Rows==Cols)
  void transpose ()
  {
    assert (Rows == Cols && "transposition in place works only if Rows==Cols");
    for (int i=0; i<Rows; ++i)
      for (int j=0; j<i; ++j)
	{
	  swap (data_[i*Cols + j], data_[j*Cols + i]);
	}
  }
  
private:
  
  double data_[Rows*Cols];
  
}; // class Mat



/// matrix addition
template <int Rows, int Cols>
Mat<Rows, Cols> operator+ (const Mat<Rows, Cols>& m1, const Mat<Rows, Cols>& m2)
{
  Mat<Rows, Cols> res(m1); // copy
  res += m2;
  return res;
}

 
/// matrix subtraction
template <int Rows, int Cols>
Mat<Rows, Cols> operator- (const Mat<Rows, Cols>& m1, const Mat<Rows, Cols>& m2)
{
  Mat<Rows, Cols> res(m1); // copy
  res -= m2;
  return res;
}


/// scale matrix
template <int Rows, int Cols>
Mat<Rows, Cols> operator* (const double factor, const Mat<Rows, Cols>& m)
{
  Mat<Rows, Cols> res(m);
  res *= factor;
  return res;
}

/// output
template <int Rows, int Cols>
std::ostream& operator<< (std::ostream& os, const Mat<Rows, Cols>& m)
{
  m.print (os);
  return os;
}

#endif // __MAT_H