BlockPreconditioner.h 4.6 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

/** \file BlockPreconditioner.h */

#ifndef AMDIS_BLOCK_PRECONDITIONER_H
#define AMDIS_BLOCK_PRECONDITIONER_H

#include "solver/ITL_Preconditioner.h"
#include "solver/SolverMatrix.h"

namespace AMDiS {

  /// Basis preconditioner structure for block-preconditioners
33 34
  template <class MatrixType, class VectorType = MTLTypes::MTLVector>
  struct BlockPreconditioner : public ITL_PreconditionerBase<MatrixType, VectorType>
35
  {
36 37 38
    typedef BlockPreconditioner                             self;
    typedef ITL_PreconditionerBase<MatrixType, VectorType>  super;
    typedef super                                           precon_base;
39 40

    BlockPreconditioner()
41 42
      : A(NULL), fullMatrix(NULL)
    { }
43

44
    virtual ~BlockPreconditioner() {}
45

46 47 48 49 50
    /// extract iranges from BlockMatrix to be used to extract sub-vectors and sub-matrices.
    void init(const SolverMatrix<Matrix<DOFMatrix*> >& A_, const MatrixType& fullMatrix_) override
    {
      A = &A_;
      fullMatrix = &fullMatrix_;
51

52 53 54 55 56 57 58 59 60 61
      BlockMapper mapper(A_);
      rows.resize(mapper.getNumComponents());
      int start = 0;
      for (int i = 0; i < mapper.getNumComponents(); i++) {
	mapper.setRow(i+1);
	int finish = mapper.row(0);
	rows[i].set(start, finish);
	start = finish;
      }
    }
62

63
    virtual void solve(const VectorType& b, VectorType& x) const
64 65 66
    { FUNCNAME("BlockPreconditioner::solve()");
      TEST_EXIT(false)("Must be implemented in derived classes!\n");
    }
67

68
    virtual void adjoint_solve(const VectorType& x, VectorType& y) const
69 70 71
    { FUNCNAME("BlockPreconditioner::adjoint_solve()");
      TEST_EXIT(false)("Must be implemented in derived classes!\n");
    }
72

73 74 75 76
    inline const mtl::irange& getRowRange(size_t i) const
    {
      return rows[i];
    }
77

78 79 80
    inline const mtl::irange& getColRange(size_t i) const
    {
      return rows[i];
81 82 83
    }

    const DOFMatrix::base_matrix_type& getSubMatrix(size_t i, size_t j) const
84 85 86
    {
      return A->getSubMatrix(i, j);
    }
87

88
    template <class SolverType, class RunnerType>
89 90
    static void createSubSolver(std::string param, SolverType*& solver, RunnerType*& runner,
				std::string solverType = "0", std::string preconType = "no",
91
				int max_iter = 100, double tol = 1.e-8)
92 93
    {
      // definition of standard-backends
Praetorius, Simon's avatar
Praetorius, Simon committed
94 95
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
      std::string backend = "p_";
96
#else
Praetorius, Simon's avatar
Praetorius, Simon committed
97 98 99 100 101 102 103
      std::string backend = "";
#endif

#if defined(HAVE_PARALLEL_PETSC) || defined(HAVE_SEQ_PETSC)
      backend += "petsc";
#else
      backend += "mtl";
104
#endif
105

106
      // === read backend-name ===
107
      std::string initFileStr = param + "->solver";
108 109

      // === read solver-name ===
Praetorius, Simon's avatar
Praetorius, Simon committed
110 111 112 113
      std::string tmp = "0";
      Parameters::get(initFileStr, tmp);
      if (tmp == "0") {
	TEST_EXIT(solverType != "0")("You have to specify the parameter %s in the initfile!\n", initFileStr.c_str());
114

Praetorius, Simon's avatar
Praetorius, Simon committed
115 116 117 118 119 120 121 122
	Parameters::set(initFileStr, solverType);
	Parameters::set(initFileStr + "->backend", backend);
	Parameters::set(initFileStr + "->left precon", preconType);
	Parameters::set(initFileStr + "->max iteration", max_iter);
	Parameters::set(initFileStr + "->tolerance", tol);
      } else {
	solverType = tmp;
      }
123

Praetorius, Simon's avatar
Praetorius, Simon committed
124
      Parameters::get(initFileStr + "->backend", backend);
125

126 127
      if (backend != "0" && backend != "no" && backend != "")
	solverType = backend + "_" + solverType;
128 129

      LinearSolverCreator *solverCreator =
130
	dynamic_cast<LinearSolverCreator*>(CreatorMap<LinearSolverInterface>::getCreator(solverType, initFileStr));
131 132 133 134 135
      TEST_EXIT(solverCreator)
	("No valid solver type found in parameter \"%s\"\n", initFileStr.c_str());
      solverCreator->setName(initFileStr);
      solver = dynamic_cast<SolverType*>(solverCreator->create());
      assert(solver != NULL);
136

137 138 139
      runner = dynamic_cast<RunnerType*>(solver->getRunner());
      assert(runner != NULL);
    }
140 141


142 143 144
  protected:
    const SolverMatrix<Matrix<DOFMatrix*> >* A;
    const MatrixType* fullMatrix;
145

146 147
    std::vector<mtl::irange> rows;
  };
148

149 150 151
} // end namespace AMDiS

#endif // AMDIS_BLOCK_PRECONDITIONER_H