HypreSolver.h 5.86 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
/******************************************************************************
 *
 * 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 HypreSolver.h */


#ifndef AMDIS_HYPRE_SOLVER_H
#define AMDIS_HYPRE_SOLVER_H

#ifdef MTL_HAS_HYPRE

#include "solver/MTL4Solver.h"
#include "solver/ITL_Preconditioner.h"
#include "MTL4Types.h"
#include "solver/itl/hypre.hpp"

#include <boost/numeric/itl/itl.hpp>
#include <boost/numeric/mtl/mtl.hpp>

namespace AMDiS {

40
  struct Hypre_Runner : public MTL4Runner< MTLTypes::MTLMatrix, MTLTypes::MTLVector >
41
42
43
44
45
  {       
    typedef MTLTypes::MTLMatrix MatrixType;
    typedef MTLTypes::MTLVector VectorType;
    typedef MTL4Runner< MatrixType, VectorType > super;
    
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
    /** Interface to the HYPRE BoomerAMG solver [...]
     * Parameters provided by AMDiS:
     * 
     * [solver]->cycle mode: 
     * 	1...V-cycle
     *	2...W-cycle
     * 
     * [solver]->interpolation type:
     *  0...classical modified interpolation
     *	1...LS interpolation (for use with GSMG)
     * 	2...classical modified interpolation for hyperbolic PDEs
     * 	3...direct interpolation (with separation of weights)
     * 	4...multipass interpolation
     * 	5...multipass interpolation (with separation of weights)
     * 	6...extended+i interpolation
     * 	7...extended+i (if no common C neighbor) interpolation
     * 	8...standard interpolation
     * 	9...standard interpolation (with separation of weights)
     * 	10..classical block interpolation (for use with nodal systems version only)
     * 	11..classical block interpolation (for use with nodal systems version only)
     * 		with diagonalized diagonal blocks
     * 	12..FF interpolation
     * 	13..FF1 interpolation
     * 	14..extended interpolation
     * 
     *  [solver]->info: 
     * 	0...no printout (default) 
     * 	1...print setup information 
     * 	2...print solve information 
     * 	3...print both setup and solve information
     * 
     *  [solver]->relax type:
     * 	0...Jacobi
     * 	1...Gauss-Seidel, sequential (very slow!)
     * 	2...Gauss-Seidel, interior points in parallel, boundary sequential (slow!)
     * 	3...hybrid Gauss-Seidel or SOR, forward solve
     * 	4...hybrid Gauss-Seidel or SOR, backward solve
     * 	5...hybrid chaotic Gauss-Seidel (works only with OpenMP)
     * 	6...hybrid symmetric Gauss-Seidel or SSOR
     * 	9...Gaussian elimination (only on coarsest level)
     * */
87
88
    Hypre_Runner(LinearSolver* oemPtr)
    : oem(*oemPtr),
89
90
      useTransposed(false),
      solverCreated(false)
91
    { 
92
      int cycleMode = -1, interpolation = -1, relaxation = -1;
93
94
      Parameters::get(oem.getName() + "->cycle mode", cycleMode);
      Parameters::get(oem.getName() + "->interpolation type", interpolation);
95
96
      Parameters::get(oem.getName() + "->relax type", relaxation);
            
97
      config.maxIter(oem.getMaxIterations());
98
99
      config.interpolationType(interpolation);
      config.relaxType(relaxation);
100
      config.cycle(cycleMode);
101
      config.tolerance(oem.getRelative());
102
103
104
105
106
      config.printLevel(oem.getInfo());
    }
    
    ~Hypre_Runner()
    {
107
      exit();
108
109
    }      
    
110
    /// Implementation of \ref MTL4Runner::init()
111
    void init(const typename super::BlockMatrix& A, const MatrixType& mtlMatrix) override
112
    {
113
114
115
116
117
118
119
120
121
122
123
      setTransposed(typename MatrixType::orientation());
      // TODO: copy matrix directly from DOFMatrix to HYPRE matrix (?)
      hypreMatrix.init(mtlMatrix);
      HYPRE_IJMatrixGetObject(hypreMatrix, (void**) &matrix);
      HYPRE_BoomerAMGCreate(&solver);
      
      mtl::dense_vector<double> swap(1, 0.0);
      mtl::HypreParVector x(swap);
      HYPRE_BoomerAMGSetup(solver, matrix, x, x);
      
      solverCreated = true;
124
125
    }
        
126
    /// Implementation of \ref MTL4Runner::solve()
127
    int solve(const MatrixType& A , VectorType& mtlX, const VectorType& mtlB) override
128
    {      
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
      mtl::HypreParVector x(mtlX);
      mtl::HypreParVector b(mtlB);
      config(solver);
      int error = 0;
      if(useTransposed)
	error = HYPRE_BoomerAMGSolveT(solver, matrix, b, x);
      else
	error = HYPRE_BoomerAMGSolve(solver, matrix, b, x);
      mtl::convert(x.getHypreVector(), mtlX);
      
      int num_iter = 0;
      HYPRE_BoomerAMGGetNumIterations(solver, &num_iter);
      oem.setIterations(num_iter);
      
      double rel_resid = 0.0;
      HYPRE_BoomerAMGGetFinalRelativeResidualNorm(solver, &rel_resid);
      oem.setRelativeResidual(rel_resid);
      
      oem.setErrorCode(error);
148
149
150
      return error;
    }
    
151
152
    /// Implementation of \ref OEMRunner::exit()
    void exit()
153
    {
154
155
156
157
158
159
160
161
162
163
164
165
166
167
      if (solverCreated)
	HYPRE_BoomerAMGDestroy(solver);
      solverCreated = false;
    }
    
  private:
    inline void setTransposed(mtl::row_major)
    { 
      useTransposed = false;
    }
    
    inline void setTransposed(mtl::col_major)
    { 
      useTransposed = true;
168
169
170
171
    }
    
  protected:
    LinearSolver& oem;
172
173
174
175
176
177
    
  private:
    HYPRE_Solver solver;
    HYPRE_ParCSRMatrix matrix;
    mtl::HypreMatrix hypreMatrix;
    
178
    mtl::AMGConfigurator config;
179
180
181
    
    bool useTransposed;
    bool solverCreated;
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
  };
  

  /**
   * \ingroup Solver
   * \class AMDiS::HypreSolver
   * \brief
   * Wrapper for the external HYPRE-AMG solver
   */
  typedef MTL4Solver< MTLTypes::MTLMatrix, MTLTypes::MTLVector, Hypre_Runner > HypreSolver;

} // namespace AMDiS

#endif // MTL_HAS_HYPRE

#endif // AMDIS_HYPRE_SOLVER_H