KrylovPreconditioner.h 5.07 KB
Newer Older
1 2 3 4 5 6 7
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
8
 * Authors:
9 10 11 12 13 14 15 16 17
 * Simon Vey, Thomas Witkowski, Andreas Naumann, Simon Praetorius, et al.
 *
 * This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
 * WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
 *
 *
 * This file is part of AMDiS
 *
 * See also license.opensource.txt in the distribution.
18
 *
19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35
 ******************************************************************************/


/** \file KrylovPreconditioner.h */

#ifndef AMDIS_KRYLOV_PRECONDITIONER_H
#define AMDIS_KRYLOV_PRECONDITIONER_H

#include "solver/ITL_Preconditioner.h"
#ifdef HAVE_PARALLEL_MTL4
#include "parallel/PITL_Solver.h"
#endif

namespace AMDiS {

  /**
   * \ingroup Solver
36
   *
37
   * \brief Implements an inner solver for (flexible) krylov subspace solvers
38
   *
39 40 41 42 43 44 45 46
   * Any AMDiS solver can be used as preconditioner for any other AMDiS solver.
   * Initfile-syntax:<br />
   * [...]->solver: any-krylov-solver<br />
   * [...]->solver->left precon: [krylov|solver]<br />
   * [...]->solver->left precon->solver: any-amdis-solver<br />
   * [...]->solver->left precon->solver->parameter1: ...<br />
   * [...]->solver->left precon->solver->parameter2: ...<br />
   * [...]->solver->left precon->solver->parameter3: ...
47
   *
48 49 50 51 52
   * This Krylov-Preconditioner can also be used as preconditioner for the inner solver, so
   * that a sequence of iner-inner-solver can be loaded.
   **/

  template< typename MatrixType, typename VectorType >
53
  struct KrylovPreconditioner : ITL_PreconditionerBase< MatrixType, VectorType >
54
  {
55 56
    typedef ITL_PreconditionerBase< MatrixType, VectorType >  precon_base;
    typedef KrylovPreconditioner< MatrixType, VectorType >    self;
57

58
    class Creator : public CreatorInterfaceName<precon_base>
59 60 61
    {
    public:
      virtual ~Creator() {}
62 63

      precon_base* create() {
64
	return new self(this->name);
65 66
      }
    };
67 68 69 70

    KrylovPreconditioner(std::string name)
      : fullMatrix(NULL),
	solver(NULL),
71
	runner(NULL)
72
    {
Praetorius, Simon's avatar
Praetorius, Simon committed
73 74 75 76 77
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
      std::string backend = "p_";
#else
      std::string backend = "";
#endif
78

Praetorius, Simon's avatar
Praetorius, Simon committed
79 80
#if defined(HAVE_PARALLEL_PETSC) || defined(HAVE_SEQ_PETSC)
      backend += "petsc";
81
#else
Praetorius, Simon's avatar
Praetorius, Simon committed
82
      backend += "mtl";
83
#endif
84

85
      // === read backend-name ===
86
      std::string initFileStr = name + "->solver";
87 88 89
      Parameters::get(initFileStr + "->backend", backend);

      // === read solver-name ===
90
      std::string solverType("0");
91
      Parameters::get(initFileStr, solverType);
92

93 94
      if (backend != "0" && backend != "no" && backend != "")
	solverType = backend + "_" + solverType;
95 96

      LinearSolverCreator *solverCreator =
97
	dynamic_cast<LinearSolverCreator*>(CreatorMap<LinearSolverInterface>::getCreator(solverType, initFileStr));
98 99 100 101
      TEST_EXIT(solverCreator)
	("No valid solver type found in parameter \"%s\"\n", initFileStr.c_str());
      solverCreator->setName(initFileStr);
      solver = solverCreator->create();
102

103
      runner = dynamic_cast<RunnerBase< MatrixType, VectorType >*>(solver->getRunner());
104
    }
105

106 107 108 109
    virtual ~KrylovPreconditioner()
    {
      delete solver;
    }
110

111
    /// Implementation of \ref ITL_PreconditionerBase::init()
112
    virtual void init(const SolverMatrix<Matrix<DOFMatrix*> >& A, const MatrixType& fullMatrix_) override
113 114 115 116
    {
      fullMatrix = &fullMatrix_;
      runner->init(A, fullMatrix_);
    }
117

118
    /// Implementation of \ref PreconditionerInterface::init()
119
    virtual void exit() override
120 121 122
    {
      runner->exit();
    }
123

124
    /// Implementation of \ref ITL_PreconditionerBase::solve()
125
    virtual void solve(const VectorType& b, VectorType& x) const override
126 127 128 129
    {
      initVector(x);
      runner->solve(*fullMatrix, x, b);
    }
130

131
    /// Implementation of \ref ITL_PreconditionerBase::adjoint_solve()
132
    virtual void adjoint_solve(const VectorType& b, VectorType& x) const override
133 134 135 136
    {
      initVector(x);
      runner->adjoint_solve(*fullMatrix, x, b);
    }
137

138
  protected: // methods
139

140 141 142 143 144 145 146 147 148
    template<typename VectorT>
    typename boost::enable_if< mtl::traits::is_distributed<VectorT>, void >::type
    initVector(VectorT& x) const
    {
#ifdef HAVE_PARALLEL_MTL4
      x.init_distribution(col_distribution(*fullMatrix), num_cols(*fullMatrix));
#endif
      set_to_zero(x);
    }
149

150 151 152 153 154 155 156
    template<typename VectorT>
    typename boost::disable_if< mtl::traits::is_distributed<VectorT>, void >::type
    initVector(VectorT& x) const
    {
      x.change_dim(num_cols(*fullMatrix));
      set_to_zero(x);
    }
157

158 159 160
  protected: // member variables

    const MatrixType* fullMatrix;
161

162 163
    LinearSolverInterface* solver;
    RunnerBase< MatrixType, VectorType >* runner;
164
  };
165 166


167 168 169 170 171
#ifdef HAVE_PARALLEL_MTL4
  typedef KrylovPreconditioner<MTLTypes::PMTLMatrix, MTLTypes::PMTLVector> KrylovPreconditionerParallel;
#endif
  typedef KrylovPreconditioner<MTLTypes::MTLMatrix, MTLTypes::MTLVector> KrylovPreconditionerSeq;

172

173 174 175
} // end namespace AMDiS

#endif // AMDIS_KRYLOV_PRECONDITIONER_H