BlockPreconditioner.h 4.15 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134
/******************************************************************************
 *
 * AMDiS - Adaptive multidimensional simulations
 *
 * Copyright (C) 2013 Dresden University of Technology. All Rights Reserved.
 * Web: https://fusionforge.zih.tu-dresden.de/projects/amdis
 *
 * Authors: 
 * 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.
 * 
 ******************************************************************************/
 

/** \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
  template<typename MatrixType>
  struct BlockPreconditioner : ITL_BasePreconditioner<MatrixType, MTLTypes::MTLVector>
  {
    typedef SolverMatrix<Matrix<DOFMatrix*> > BlockMatrix;
    typedef ITL_BasePreconditioner<MatrixType, MTLTypes::MTLVector>  super;
    typedef BlockPreconditioner<MatrixType>                          self;
    
    BlockPreconditioner() 
      : A(NULL), fullMatrix(NULL)
    { }
    
    virtual ~BlockPreconditioner() {}
    
    /// 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_;
      
      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;
      }
    }
  
    virtual void solve(const MTLTypes::MTLVector& b, MTLTypes::MTLVector& x) const
    { FUNCNAME("BlockPreconditioner::solve()");
      TEST_EXIT(false)("Must be implemented in derived classes!\n");
    }
    
    virtual void adjoint_solve(const MTLTypes::MTLVector& x, MTLTypes::MTLVector& y) const
    { FUNCNAME("BlockPreconditioner::adjoint_solve()");
      TEST_EXIT(false)("Must be implemented in derived classes!\n");
    }
    
    inline const mtl::irange& getRowRange(size_t i) const
    {
      return rows[i];
    }
    
    inline const mtl::irange& getColRange(size_t i) const
    {
      return rows[i];
    } 
    
    const DOFMatrix::base_matrix_type& getSubMatrix(size_t i, size_t j) const 
    {
      return A->getSubMatrix(i, j);
    }
    
    template<typename SolverType, typename RunnerType>
    static void createSubSolver(std::string param, SolverType*& solver, RunnerType*& runner, std::string solverType = "0")
    {
      // definition of standard-backends
#if defined HAVE_PARALLEL_PETSC
      std::string backend("p_petsc");
#elif defined HAVE_PARALLEL_MTL
      std::string backend("p_mtl");
#elif defined HAVE_PETSC
      std::string backend("petsc");
#else
      std::string backend("mtl");
#endif
    
      // === read backend-name ===
      std::string initFileStr = param + "->solver";    
      Parameters::get(initFileStr + "->backend", backend);

      // === read solver-name ===
      Parameters::get(initFileStr, solverType);
      
      if (backend != "0" && backend != "no" && backend != "")
	solverType = backend + "_" + solverType;
    
      LinearSolverCreator *solverCreator = 
	dynamic_cast<LinearSolverCreator*>(CreatorMap<LinearSolver>::getCreator(solverType, initFileStr));
      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);
      
      runner = dynamic_cast<RunnerType*>(solver->getRunner());
      assert(runner != NULL);
    }
    
    
  protected:
    const SolverMatrix<Matrix<DOFMatrix*> >* A;
    const MatrixType* fullMatrix;
    
    std::vector<mtl::irange> rows;
  };
  
} // end namespace AMDiS

#endif // AMDIS_BLOCK_PRECONDITIONER_H