// ============================================================================
// ==                                                                        ==
// == AMDiS - Adaptive multidimensional simulations                          ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
// ==  TU Dresden                                                            ==
// ==                                                                        ==
// ==  Institut f�r Wissenschaftliches Rechnen                               ==
// ==  Zellescher Weg 12-14                                                  ==
// ==  01069 Dresden                                                         ==
// ==  germany                                                               ==
// ==                                                                        ==
// ============================================================================
// ==                                                                        ==
// ==  https://gforge.zih.tu-dresden.de/projects/amdis/                      ==
// ==                                                                        ==
// ============================================================================

/** \file UmfPackSolver.h */

#ifndef AMDIS_UMFPACKSOLVER_H
#define AMDIS_UMFPACKSOLVER_H

#if defined HAVE_UMFPACK && defined MTL_HAS_UMFPACK

#include <iostream>
#include <boost/numeric/mtl/operation/two_norm.hpp>
#include <boost/numeric/mtl/interface/umfpack_solve.hpp>
#include "OEMSolver.h"

namespace AMDiS {

  /**
   * \ingroup Solver
   * 
   * \brief
   * Wrapper for the external UMFPACK solver:
   *   http://www.cise.ufl.edu/research/sparse/umfpack/
   *
   * This is a direct solver for large sparse matrices.
   */

  class UmfPackSolver : public OEMSolver
  {
  public:
    /// Creator class used in the OEMSolverMap.
    class Creator : public OEMSolverCreator
    {
    public:
      virtual ~Creator() {}
      
      /// Returns a new UmfPackSolver object.
      OEMSolver* create() 
      { 
	return new UmfPackSolver(this->name); 
      }
    };

    /// Constructor
    UmfPackSolver(std::string name) 
      : OEMSolver(name), 
	solver(0),
	store_symbolic(0),
	symmetric_strategy(0)
    {
      GET_PARAMETER(0, name + "->store symbolic", "%d", &store_symbolic);
      GET_PARAMETER(0, name + "->symmetric strategy", "%d", &symmetric_strategy);
    }
    
    /// Destructor
    ~UmfPackSolver() 
    { 
      if (solver) 
	delete solver;
    }

    /// Solves the system directly
    int solveSystem(const DOFMatrix::base_matrix_type& A,
		    mtl::dense_vector<value_type>& x,
		    const mtl::dense_vector<value_type>& b)
    {
      if (!solver) {
	if (!symmetric_strategy)
	  solver = new mtl::matrix::umfpack::solver<matrix_type>(A);
	else
	  solver = new mtl::matrix::umfpack::solver<matrix_type>(A, UMFPACK_STRATEGY_SYMMETRIC);
      } else {
	if (!multipleRhs)
 	  if (store_symbolic)
 	    solver->update_numeric();
 	  else
	    solver->update();		
     }
     
      int code = (*solver)(x, b);
      mtl::dense_vector<value_type> r(b); 
      r -= A * x; 
      residual = two_norm(r);
      std::cout << "UmfPackSolver: ||b-Ax|| = " << residual << "\n";

      return code;
    }

    
  private:
    mtl::matrix::umfpack::solver<matrix_type> *solver;

    int store_symbolic;

    int symmetric_strategy;
  };
  
}

#endif // HAVE_UMFPACK

#endif // AMDIS_UMFPACKSOLVER_H