paralution_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_PARALUTION_SOLVER_H_
#define __HERMES_COMMON_PARALUTION_SOLVER_H_
#include "config.h"
#ifdef WITH_PARALUTION
#include "solvers/linear_matrix_solver.h"
#include "algebra/cs_matrix.h"

#include "paralution.hpp"
#include "solvers/precond.h"

using namespace Hermes::Algebra;

namespace Hermes
{
  namespace Algebra
  {
    using namespace Hermes::Solvers;

    enum ParalutionMatrixType
    {
      ParalutionMatrixTypeCSR,
      // ParalutionMatrixTypeBCSR,
      // ParalutionMatrixTypeMCSR,
      // ParalutionMatrixTypeCOO,
      // ParalutionMatrixTypeDIA,
      // ParalutionMatrixTypeELL,
      // ParalutionMatrixTypeHYB,
      // ParalutionMatrixTypeDENSE
    };

    template <typename Scalar>
    class HERMES_API ParalutionMatrix : public CSRMatrix < Scalar >
    {
    public:
      ParalutionMatrix(ParalutionMatrixType type = ParalutionMatrixTypeCSR);
      virtual ~ParalutionMatrix();

      virtual void free();
      virtual void zero();
      virtual void alloc();

      paralution::LocalMatrix<Scalar>& get_paralutionMatrix();

    private:
      paralution::LocalMatrix<Scalar> paralutionMatrix;
      ParalutionMatrixType paralutionMatrixType;

      // Friends.
      template<typename T> friend SparseMatrix<T>*  create_matrix();
    };

    template <typename Scalar>
    class HERMES_API ParalutionVector : public SimpleVector < Scalar >
    {
    public:
      ParalutionVector();
      ParalutionVector(unsigned int size);
      virtual ~ParalutionVector();
      virtual void alloc(unsigned int ndofs);
      virtual void free();
      virtual void zero();

      paralution::LocalVector<Scalar>* get_paralutionVector();

    private:
      paralution::LocalVector<Scalar>* paralutionVector;
    };
  }

  namespace Preconditioners
  {
    template <typename Scalar>
    class ParalutionPrecond : public Hermes::Preconditioners::Precond < Scalar >
    {
    public:
      ParalutionPrecond(PreconditionerType preconditionerType);

      paralution::Preconditioner<paralution::LocalMatrix<Scalar>, paralution::LocalVector<Scalar>, Scalar>& get_paralutionPreconditioner();
      static paralution::Preconditioner<paralution::LocalMatrix<Scalar>, paralution::LocalVector<Scalar>, Scalar>* return_paralutionPreconditioner(PreconditionerType preconditionerType);
    private:
      paralution::Preconditioner<paralution::LocalMatrix<Scalar>, paralution::LocalVector<Scalar>, Scalar>* paralutionPreconditioner;

#if __PARALUTION_VER >= 500
      paralution::FSAI<paralution::LocalMatrix<Scalar>, paralution::LocalVector<Scalar>, Scalar>* saddlePoint_p_k;
      paralution::SPAI<paralution::LocalMatrix<Scalar>, paralution::LocalVector<Scalar>, Scalar>* saddlePoint_p_s;
      paralution::MultiColoredILU<paralution::LocalMatrix<Scalar>, paralution::LocalVector<Scalar>, Scalar >* mcilu_p;
#endif
    };
  }

  namespace Solvers
  {
    class ParalutionInitialization
    {
    public:
      static void init_paralution()
      {
        paralution::init_paralution();
        paralution::set_omp_threads_paralution(HermesCommonApi.get_integral_param_value(numThreads));
        paralution::info_paralution();
      }

      static void set_threads_paralution()
      {
        if (HermesCommonApi.get_integral_param_value(matrixSolverType) == SOLVER_PARALUTION_AMG || HermesCommonApi.get_integral_param_value(matrixSolverType) == SOLVER_PARALUTION_ITERATIVE)
        {
          paralution::set_omp_threads_paralution(HermesCommonApi.get_integral_param_value(numThreads));
          paralution::info_paralution();
        }
      }

      static void deinit_paralution()
      {
        paralution::stop_paralution();
      }
    };

    template <typename Scalar>
    class HERMES_API AbstractParalutionLinearMatrixSolver : public virtual LoopSolver < Scalar >
    {
    public:
      virtual void solve(Scalar* initial_guess);
      virtual void solve();

      virtual int get_num_iters();

      virtual double get_residual_norm();

      virtual int get_matrix_size();

      virtual void free();

    protected:
      AbstractParalutionLinearMatrixSolver();
      AbstractParalutionLinearMatrixSolver(ParalutionMatrix<Scalar> *m, ParalutionVector<Scalar> *rhs);
      virtual ~AbstractParalutionLinearMatrixSolver();

      void reset_internal_solver();

      virtual void init_internal_solver() = 0;

      ParalutionMatrix<Scalar> *matrix;

      ParalutionVector<Scalar> *rhs;

      // Store num_iters.
      int num_iters;

      // Store final_residual.
      double final_residual;

      paralution::IterativeLinearSolver<paralution::LocalMatrix<Scalar>, paralution::LocalVector<Scalar>, Scalar>* paralutionSolver;

    private:
      void presolve_init();

      template<typename T> friend LinearMatrixSolver<T>* create_linear_solver(Matrix<T>* matrix, Vector<T>* rhs, bool use_direct_solver);
      template<typename T> friend class AMGParalutionLinearMatrixSolver;
    };

    template <typename Scalar>
    class HERMES_API IterativeParalutionLinearMatrixSolver : public AbstractParalutionLinearMatrixSolver<Scalar>, public virtual IterSolver < Scalar >
    {
    public:
      IterativeParalutionLinearMatrixSolver();
      IterativeParalutionLinearMatrixSolver(ParalutionMatrix<Scalar> *m, ParalutionVector<Scalar> *rhs);
      virtual ~IterativeParalutionLinearMatrixSolver();

      void set_solver_type(IterSolverType iterSolverType);

      virtual void init_internal_solver();

      void free();

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

      // Linear Solver creation.
      static paralution::IterativeLinearSolver<paralution::LocalMatrix<Scalar>, paralution::LocalVector<Scalar>, Scalar>* return_paralutionSolver(IterSolverType type);

    private:
      Preconditioners::ParalutionPrecond<Scalar> *preconditioner;
    };

    template <typename Scalar>
    class HERMES_API AMGParalutionLinearMatrixSolver : public AbstractParalutionLinearMatrixSolver<Scalar>, public virtual AMGSolver < Scalar >
    {
    public:
      AMGParalutionLinearMatrixSolver(ParalutionMatrix<Scalar> *m, ParalutionVector<Scalar> *rhs);
      virtual ~AMGParalutionLinearMatrixSolver();

      virtual void set_smoother(IterSolverType solverType, PreconditionerType preconditionerType);

      virtual void init_internal_solver();
    };
  }
}
#endif
#endif