/**************************************************************************
**       Title: slow_code
**    $RCSfile: $
**   $Revision:$$Name:$
**       $Date: 27.6.2006$
**   Copyright: GPL $Author: Bernard Haasdonk$
** Description: program, which demonstrates various efficiency-degrading
**              programming methods
**************************************************************************/

#include <vector>
#include <fstream>
#include <iostream>
#include <time.h>

using namespace std;

/*************************************************************/
// pure virtual interface class providing some basis functions
/*************************************************************/
class MatrixInterface
{
public:
  MatrixInterface() 
        {}  
  ~MatrixInterface() 
        {}
//  virtual 
  virtual void resize(int nrows, int ncols) = 0;
  virtual int nrows() const = 0;
  virtual int ncols() const = 0;
};

/*************************************************************/
// implementation class with interface methods plus some more
/*************************************************************/
class SlowMatrix: public MatrixInterface
{
public:
  typedef vector<double> RowType;  
  
//  SlowMatrix()
//        {
          // cout << "entering constructor \n";          
          // do nothing
//        }
  
  virtual 
  void resize(int nrows, int ncols)
        {
          // cout << "entering resize \n";          
          data_.resize(nrows);
          for (int i=0;i<nrows;i++)
              data_[i].resize(ncols);                   
        }
     
  virtual 
  int nrows() const
        {
          // cout << "entering nrows \n";          
          return data_.size();
        }
  
  virtual 
  int ncols() const
        {
          // cout << "entering ncols \n";          
          if(nrows()>0)
              return data_[0].size();  
          else 
              return 0;
        }

  RowType& operator[](int rownr)
        {
          // cout << "entering operator[] \n";          
          assert(rownr < data_.size());
          return data_[rownr];
        }

  const RowType& operator[](int rownr) const
        {
          // cout << "entering const operator[] \n";          
          assert(rownr < data_.size());
          return data_[rownr];
        }

  
  SlowMatrix& operator = (const SlowMatrix &rhs)
        { 
          // cout << "entering operator= \n";          
          data_.resize(rhs.nrows());
          for (int ri=0;ri< rhs.nrows(); ri++)
              data_[ri] = rhs[ri];
          return *this;
        }
  
  SlowMatrix& operator = (const double val)
        {
          // cout << "entering setValue \n";          
          for (int ri=0;ri<nrows();ri++)
              for (int ci=0;ci<ncols();ci++)
              {
//                cout << "setting entry " << ri << "," << ci <<"\n";
                data_[ri][ci]= val;
              } 
        }
  
  const SlowMatrix operator*(SlowMatrix& m) 
        {
          // cout << "entering operator* \n";          
          assert(ncols()==m.nrows());
          SlowMatrix result;
          result.resize(nrows(),m.ncols());
          result = 0.0;
          
          for (int c=0;c<m.ncols();c++) {
            for (int r=0;r<nrows();r++) {
              for (int i=0;i<ncols();i++) {
                result[r][c]+=data_[r][i]*m[i][c];
              }
            }
          }
          return result;
        }
 

  const SlowMatrix operator*(const double val) const
        {
          // cout << "entering operator* \n";          
          SlowMatrix result;
          result.resize(nrows(),ncols());
          
          for (int c=0;c<ncols();c++) 
              for (int r=0;r<nrows();r++) 
                  result[r][c]+=data_[r][c]*  val;
          return result;
        }
 
  const SlowMatrix operator+(SlowMatrix& m) 
        {
          // cout << "entering operator+ \n";          
          assert(ncols()==m.ncols());
          assert(nrows()==m.nrows());
          SlowMatrix result;
          result.resize(nrows(),ncols());
          
          for (int c=0;c<m.ncols();c++) 
              for (int r=0;r<nrows();r++) 
                  result[r][c]+=data_[r][c]+m[r][c];
          
          return result;
        }
 
  friend std::ostream& operator<<(std::ostream &os,
                                  const SlowMatrix &mat);  
  
private:
  vector< vector< double> > data_;

};


/*************************************************************/
// friend output function
/*************************************************************/
std::ostream& operator <<(std::ostream &os,
                          const SlowMatrix &mat)
{
  os << "number of rows: "<< mat.nrows() << " \n";
  os << "number of columns: "<< mat.ncols() << " \n";
  
  os << "[\n";
  for (int ri=0;ri<mat.nrows();ri++)
  {
    for (int ci=0;ci<mat.ncols();ci++)
        os << mat[ri][ci] << ", ";
    os << "\n";    
  }
  os << "]\n";
  return os;
}

/*************************************************************/
// main program doing some (useless) matrix operations...
/*************************************************************/
int main(int argc, char* argv[])
{
  // the following parameters can be chosen freely:
  const int s = 100;       // size of the matrix
  const double eps = 1e-4; // off-diagonal element value
  const int loops = 10;    // number of arithmetic loops

  clock_t start, end;

  start = clock();    
  
  SlowMatrix mat1;
  mat1.resize(s,s);  
  mat1 = eps;
  // set diagonal
  for (int i=0; i<s;i++)
      mat1[i][i] = 0.5;

  SlowMatrix mat2;
  mat2 = mat1;
  
  // do some matrix-arithmetic operations
  for (int l=0;l< loops; l++)
  {
    SlowMatrix tmp1 = mat2 * mat2;
    SlowMatrix tmp2 = mat2 * mat1;
    SlowMatrix tmp3 = mat1 * mat2;
    SlowMatrix tmp4 = mat1 * mat1;
    mat2 = tmp1;
    mat2 = mat2 + tmp2;
    mat2 = mat2 + tmp3;
    mat2 = mat2 + tmp4;
    mat2 = mat2 * 0.5;
  }
  
  cout << mat2; 

  end = clock();
  cout << "time required = " << 
      (end-start)*1.0/CLOCKS_PER_SEC << " sec.\n";
}