linear_matrix_solver.h

Go to the documentation of this file.

// This file is part of HermesCommon
//
// Copyright (c) 2009 hp-FEM group at the University of Nevada, Reno (UNR).
// Email: hpfem-group@unr.edu, home page: http://www.hpfem.org/.
//
// Hermes2D is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published
// by the Free Software Foundation; either version 2 of the License,
// or (at your option) any later version.
//
// Hermes2D is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with Hermes2D; if not, write to the Free Software
// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
#ifndef __HERMES_COMMON_LINEAR_MATRIX_SOLVER_H_
#define __HERMES_COMMON_LINEAR_MATRIX_SOLVER_H_

#include "precond.h"
#include "exceptions.h"
#include "algebra/cs_matrix.h"
#include "algebra/vector.h"
#include "mixins.h"
#include "algebra/algebra_mixins.h"

using namespace Hermes::Algebra;

namespace Hermes
{
  namespace Solvers
  {
    using namespace Hermes::Preconditioners;
    enum MatrixStructureReuseScheme
    {
      HERMES_CREATE_STRUCTURE_FROM_SCRATCH,              
      HERMES_REUSE_MATRIX_REORDERING,             
      HERMES_REUSE_MATRIX_REORDERING_AND_SCALING, 
      HERMES_REUSE_MATRIX_STRUCTURE_COMPLETELY       
    };

    template <typename Scalar> class DirectSolver;
    template <typename Scalar> class LoopSolver;
    template <typename Scalar> class IterSolver;
    template <typename Scalar> class AMGSolver;

    template <typename Scalar>
    class HERMES_API LinearMatrixSolver : public Hermes::Mixins::Loggable, public Hermes::Mixins::TimeMeasurable
    {
    public:
      LinearMatrixSolver(SparseMatrix<Scalar>* matrix, Hermes::Algebra::Vector<Scalar>* rhs);
      DirectSolver<Scalar>* as_DirectSolver() const;
      LoopSolver<Scalar>* as_LoopSolver() const;
      IterSolver<Scalar>* as_IterSolver() const;
      AMGSolver<Scalar>* as_AMGSolver() const;
      virtual ~LinearMatrixSolver();
      virtual void free() = 0;

      virtual void solve() = 0;

      virtual void solve(Scalar* initial_guess) = 0;

      Scalar *get_sln_vector();

      double get_time();

      virtual int get_matrix_size() = 0;

      virtual MatrixStructureReuseScheme get_used_reuse_scheme() const;

      virtual void set_reuse_scheme(MatrixStructureReuseScheme reuse_scheme);

      virtual void set_reuse_scheme();

      virtual void use_node_wise_ordering(unsigned int num_pdes);
      virtual void use_equations_wise_ordering();

      virtual double get_residual_norm();

      SparseMatrix<Scalar>* get_matrix();
      Vector<Scalar>* get_rhs();

    protected:
      MatrixStructureReuseScheme reuse_scheme;

      SparseMatrix<Scalar>* general_matrix;
      Vector<Scalar>* general_rhs;

      Scalar *sln;

      double time;

      unsigned int n_eq;

      bool node_wise_ordering;
    };

    template <typename Scalar>
    class HERMES_API ExternalSolver : public LinearMatrixSolver<Scalar>, public Algebra::Mixins::MatrixRhsOutput < Scalar >
    {
    public:
      typedef ExternalSolver<Scalar>* (*creation)(CSCMatrix<Scalar> *m, SimpleVector<Scalar> *rhs);
      static creation create_external_solver;

      ExternalSolver(CSCMatrix<Scalar> *m, SimpleVector<Scalar> *rhs);
      virtual void solve() { throw Exceptions::MethodNotOverridenException("ExternalSolver::solve()."); };
      virtual void solve(Scalar* initial_guess) { throw Exceptions::MethodNotOverridenException("ExternalSolver::solve()."); };
      virtual int get_matrix_size() { return this->m->get_size(); };
      CSCMatrix<Scalar> *get_matrix() { return this->m; }
      SimpleVector<Scalar> *get_rhs() { return this->rhs; }

    protected:
      CSCMatrix<Scalar> *m;
      SimpleVector<Scalar> *rhs;
    };

    template <typename Scalar>
    class HERMES_API SimpleExternalSolver : public ExternalSolver < Scalar >
    {
    public:
      SimpleExternalSolver(CSCMatrix<Scalar> *m, SimpleVector<Scalar> *rhs);
      void solve();
      void solve(Scalar* initial_guess);

    protected:
      virtual std::string command() = 0;
    };

    template <typename Scalar>
    class HERMES_API DirectSolver : public LinearMatrixSolver < Scalar >
    {
    public:
      DirectSolver(SparseMatrix<Scalar>* matrix, Vector<Scalar>* rhs);
      virtual void solve() = 0;
      virtual void solve(Scalar* initial_guess);

      virtual double get_residual_norm() { return 0.; };
    };

    enum LoopSolverToleranceType
    {
      AbsoluteTolerance = 0,
      RelativeTolerance = 1,
      DivergenceTolerance = 2
    };

    template <typename Scalar>
    class HERMES_API LoopSolver : public LinearMatrixSolver < Scalar >
    {
    public:
      LoopSolver(SparseMatrix<Scalar>* matrix, Vector<Scalar>* rhs);

      virtual int get_num_iters() = 0;

      virtual double get_residual_norm() = 0;

      virtual void set_tolerance(double tol);

      virtual void set_tolerance(double tolerance, LoopSolverToleranceType toleranceType);

      virtual void set_max_iters(int iters);

    protected:
      int max_iters;
      double tolerance;
      LoopSolverToleranceType toleranceType;
    };

    enum IterSolverType
    {
      CG = 0,
      GMRES = 1,
      BiCGStab = 2,
      CR = 3,
      IDR = 4
    };

    template <typename Scalar>
    class HERMES_API IterSolver : public virtual LoopSolver < Scalar >
    {
    public:
      IterSolver(SparseMatrix<Scalar>* matrix, Vector<Scalar>* rhs);

      virtual void set_precond(Precond<Scalar> *pc) = 0;

      void set_solver_type(IterSolverType iterSolverType);

    protected:
      bool precond_yes;

      // Paralution solver type.
      IterSolverType iterSolverType;
    };

    template <typename Scalar>
    class HERMES_API AMGSolver : public virtual LoopSolver < Scalar >
    {
    public:
      AMGSolver(SparseMatrix<Scalar>* matrix, Vector<Scalar>* rhs);

      virtual void set_smoother(IterSolverType solverType, PreconditionerType preconditionerType);

    protected:
      IterSolverType smootherSolverType;
      PreconditionerType smootherPreconditionerType;
    };

    template<typename Scalar>
    HERMES_API LinearMatrixSolver<Scalar>*
      create_linear_solver(Matrix<Scalar>* matrix, Vector<Scalar>* rhs, bool use_direct_solver = false);
  }
} // End of documentation group Solvers.
#endif