// start with including some headers
#include "config.h"

#include<iostream>               // for input/output to shell
#include<fstream>                // for input/output to files
#include<vector>                 // STL vector class
#include<complex>

#include<math.h>                 // Yes, we do some math here
#include<stdio.h>                // There is nothing better than sprintf
#include<sys/times.h>            // for timing measurements

#include<dune/istl/istlexception.hh>
#include<dune/istl/basearray.hh>
#include<dune/common/fvector.hh>
#include<dune/common/fmatrix.hh>
#include<dune/istl/bvector.hh>
#include<dune/istl/vbvector.hh>
#include<dune/istl/bcrsmatrix.hh>
#include<dune/istl/io.hh>
#include<dune/istl/gsetc.hh>
#include<dune/istl/ilu.hh>
#include<dune/istl/operators.hh>
#include<dune/istl/solvers.hh>
#include<dune/istl/preconditioners.hh>
#include<dune/istl/scalarproducts.hh>

// a simple stop watch
class Timer
{
public:
  Timer ()
  {
	struct tms buf;
	cstart = times(&buf);
  }

  void start ()
  {
	struct tms buf;
	cstart = times(&buf);
  }

  double stop ()
  {
	struct tms buf;
	cend = times(&buf);
	return ((double)(cend-cstart))/100.0;
  }

  double gettime ()
  {
	return ((double)(cend-cstart))/100.0;
  }

private:
  clock_t cstart,cend;  
};


void test_BCRSMatrix ()
{
  const int N=5, M =4,M2 =5, L = 4;
  typedef double ValueType;
  typedef Dune::BCRSMatrix<Dune::FieldMatrix<ValueType,1,1> > EmpMatrixTopoElemOrder;
  typedef Dune::BCRSMatrix<EmpMatrixTopoElemOrder>            EmpMatrixTopoElem;


  Dune::FieldMatrix<ValueType,1,1> F;
  F = 42.0;


  EmpMatrixTopoElemOrder E(N,N,N*N,EmpMatrixTopoElemOrder::row_wise);
  for (EmpMatrixTopoElemOrder::CreateIterator i=E.createbegin(); i!=E.createend(); ++i)
    for (int j=0; j<N; ++j)
    i.insert(j);

  for (EmpMatrixTopoElemOrder::RowIterator i=E.begin(); i!=E.end(); ++i)
  for (EmpMatrixTopoElemOrder::ColIterator j=(*i).begin(); j!=(*i).end(); ++j)
    *j = F;

  std::cout << "============================================." << std::endl;
  printmatrix(std::cout,E,"E","row",10,2);
  std::cout << "=============================================."<< std::endl;


  EmpMatrixTopoElem D(M,M,2,EmpMatrixTopoElem::row_wise);
  for (EmpMatrixTopoElem::CreateIterator i=D.createbegin(); i!=D.createend(); ++i)
    for (int j=0; j<M; ++j)
    {
      if ((i.index() == 1) && (j == 1))
      {
        i.insert(j);
      }
    }
    
  for (EmpMatrixTopoElem::RowIterator i=D.begin(); i!=D.end(); ++i)
  for (EmpMatrixTopoElem::ColIterator j=(*i).begin(); j!=(*i).end(); ++j)
    *j = E;

  std::cout << "============================================." << std::endl;
  printmatrix(std::cout,D,"D","row",10,2);
  std::cout << "=============================================."<< std::endl;


}


int main (int argc , char ** argv)
{
  try {
 	test_BCRSMatrix();
  }
  catch (Dune::ISTLError& error)
	{
	  std::cout << error << std::endl;
	}
  catch (Dune::Exception& error)
	{
	  std::cout << error << std::endl;
	}
  catch (const std::bad_alloc& e)
	{
	  std::cout << "memory exhausted" << std::endl;
	}
  catch (...)
	{
	  std::cout << "unknown exception caught" << std::endl;
	}

  return 0;
}
