#ifdef HAVE_DUNE

#define DUNE_MINIMAL_DEBUG_LEVEL 4

#include <dune/common/fmatrix.hh>
#include <dune/istl/bvector.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/ilu.hh>
#include <dune/istl/operators.hh>
#include <dune/istl/solvers.hh>
#include <dune/istl/preconditioners.hh>
#include <dune/istl/matrix.hh>

#include "DuneSolver.h"
#include "DOFVector.h"

namespace AMDiS {

  template<>
  int DuneSolver< DOFVector<double> >::solveSystem(MatVecMultiplier< DOFVector<double> > *mv,
 						   DOFVector<double> *x,
 						   DOFVector<double> *b)
  {
    FUNCNAME("DuneSolver::solveSystem()");

    TEST_EXIT(x->getSize() == b->getSize())("Vectors x and b must have the same size!");

    StandardMatVec<DOFMatrix, DOFVector<double> > *stdMatVec = 
      dynamic_cast<StandardMatVec<DOFMatrix, DOFVector<double> > *>(mv);
  
    DOFMatrix *m = stdMatVec->getMatrix();
    int nRows = m->getFESpace()->getAdmin()->getUsedSize();

    DuneMatrix duneMatrix(nRows, nRows, DuneMatrix::random);
    DuneVector duneVecX(nRows);
    DuneVector duneVecB(nRows);

    mapDOFMatrix(m, &duneMatrix);
    mapDOFVector(x, &duneVecX);
    mapDOFVector(b, &duneVecB);
      
    MSG("solving system with DUNE ...\n");
    
    Dune::MatrixAdapter<DuneMatrix, DuneVector, DuneVector> op(duneMatrix);
    Dune::InverseOperatorResult opRes;
    Dune::Preconditioner<DuneVector, DuneVector> *dunePreconditioner = 
      getPreconditioner<DuneMatrix, DuneVector>(&duneMatrix);
    Dune::InverseOperator<DuneVector, DuneVector> *duneSolver = 
      getSolver<DuneMatrix, DuneVector>(&op, dunePreconditioner);

    duneSolver->apply(duneVecX, duneVecB, opRes);

    DuneVector::Iterator duneVecIt;    
    DOFVector<double>::Iterator vecXIter(x, USED_DOFS);
    for (vecXIter.reset(), duneVecIt = duneVecX.begin(); 
	 !vecXIter.end(); 
	 ++vecXIter, ++duneVecIt) {
      *vecXIter = *duneVecIt;
    }

    delete dunePreconditioner;
    delete duneSolver;
   
    return opRes.iterations;
  }


  template<>
  int DuneSolver<SystemVector>::solveSystem(MatVecMultiplier<SystemVector> *mv,
					    SystemVector *x, SystemVector *b)
  {
    FUNCNAME("DuneSolver::solveSystem()");

    // Calculate residual reduction which should be achieved by the dune solver.
    *p = *x;
    *p *= -1.0;
    mv->matVec(NoTranspose, *p, *r);
    *r += *b;
    double reduction = this->tolerance / norm(r);

    typedef Dune::Matrix<DuneMatrix> SystemDuneMatrix;
    typedef Dune::BlockVector<DuneVector> SystemDuneVector;

    TEST_EXIT(x->getSize() == b->getSize())("Vectors x and b must have the same size!");

    StandardMatVec<Matrix<DOFMatrix*>, SystemVector> *stdMatVec = 
      dynamic_cast<StandardMatVec<Matrix<DOFMatrix*>, SystemVector> *>(mv);
    Matrix<DOFMatrix*> *m = stdMatVec->getMatrix();

    // Number of systems.
    int nComponents = m->getSize();  
    // Size of one DOF matrix.
    int oneMatrixSize = ((*m)[0][0])->getFESpace()->getAdmin()->getUsedSize();

    SystemDuneMatrix duneMatrix(nComponents, nComponents);
    SystemDuneVector duneVecX(nComponents);
    SystemDuneVector duneVecB(nComponents);

    for (int i = 0; i < nComponents; i++) {
      for (int j = 0; j < nComponents; j++) {
	if ((*m)[i][j]) {
	  duneMatrix[i][j].setSize(oneMatrixSize, oneMatrixSize);
	  duneMatrix[i][j].setBuildMode(DuneMatrix::random);

	  mapDOFMatrix((*m)[i][j], &duneMatrix[i][j]);
	} else {
	  duneMatrix[i][j].setSize(0, 0);
	}
      }

      duneVecX[i].resize(oneMatrixSize, false);      
      mapDOFVector(x->getDOFVector(i), &duneVecX[i]);      

      duneVecB[i].resize(oneMatrixSize, false);
      mapDOFVector(b->getDOFVector(i), &duneVecB[i]);
    }

    MSG("solving system with DUNE ...\n");
    
    Dune::MatrixAdapter<SystemDuneMatrix, SystemDuneVector, SystemDuneVector> op(duneMatrix);
    Dune::InverseOperatorResult opRes;
    Dune::Preconditioner<SystemDuneVector, SystemDuneVector> *dunePreconditioner =
      getSystemPreconditioner<SystemDuneMatrix, SystemDuneVector>(&duneMatrix);
    Dune::InverseOperator<SystemDuneVector, SystemDuneVector> *duneSolver =
      getSolver<SystemDuneMatrix, SystemDuneVector>(&op, dunePreconditioner);

    duneSolver->apply(duneVecX, duneVecB, reduction, opRes);

    for (int i = 0; i < nComponents; i++) {
      DuneVector::Iterator duneVecIt;    
      DOFVector<double>::Iterator vecXIter(x->getDOFVector(i), USED_DOFS);
      for (vecXIter.reset(), duneVecIt = duneVecX[i].begin(); 
	   !vecXIter.end(); 
	   ++vecXIter, ++duneVecIt) {
	*vecXIter = *duneVecIt;
      }
    }

    // Calculate and print the residual.
    *p = *x;
    *p *= -1.0;
    mv->matVec(NoTranspose, *p, *r);
    *r += *b;
    
    this->residual = norm(r);

    MSG("Residual: %e\n", this->residual);

    delete dunePreconditioner;
    delete duneSolver;
   
    return opRes.iterations;
  }
}

#endif // HAVE_DUNE