FEAT 3
Finite Element Analysis Toolbox
Loading...
Searching...
No Matches
hypre.cpp
1// FEAT3: Finite Element Analysis Toolbox, Version 3
2// Copyright (C) 2010 by Stefan Turek & the FEAT group
3// FEAT3 is released under the GNU General Public License version 3,
4// see the file 'copyright.txt' in the top level directory for details.
5
6// includes, FEAT
7#include <kernel/base_header.hpp>
8
9#ifdef FEAT_HAVE_HYPRE
10
11// preprocessor macro sanity checks
12#ifdef FEAT_HAVE_MPI
13# ifdef HYPRE_SEQUENTIAL
14# error FEAT_HAVE_MPI and HYPRE_SEQUENTIAL are both defined
15# endif
16# ifndef HYPRE_HAVE_MPI
17# define HYPRE_HAVE_MPI 1
18# endif
19#else
20# ifdef HYPRE_HAVE_MPI
21# error FEAT_HAVE_MPI is not defined, but HYPRE_HAVE_MPI is defined
22# endif
23# ifndef HYPRE_SEQUENTIAL
24# define HYPRE_SEQUENTIAL 1
25# endif
26#endif // FEAT_HAVE_MPI
27
28// includes, HYPRE third-party lib
29#include <HYPRE.h>
30#include <HYPRE_parcsr_ls.h>
31
32// includes, system
33#include <vector>
34
35namespace FEAT
36{
37 namespace Solver
38 {
39 namespace Hypre
40 {
44 class Core
45 {
46 public:
47 // Note: if compiled without MPI, HYPRE defines its own MPI_Comm dummy type
48 MPI_Comm _comm;
50 std::vector<HYPRE_Int> _hypre_dof_idx;
52 std::vector<HYPRE_Int> _hypre_num_nze;
54 std::vector<HYPRE_Int> _hypre_col_idx;
55
56 // note: the following members are actually pointers
58 HYPRE_IJMatrix _hypre_matrix;
60 HYPRE_IJVector _hypre_vec_def, _hypre_vec_cor;
62 HYPRE_ParCSRMatrix _hypre_parcsr;
64 HYPRE_ParVector _hypre_par_def, _hypre_par_cor;
65
66
67 template<typename IT_>
68 explicit Core(const void* comm, Index my_dof_offset, Index num_owned_dofs, const IT_* row_ptr, const IT_* col_idx) :
69#ifdef FEAT_HAVE_MPI
70 _comm(*reinterpret_cast<const MPI_Comm*>(comm)),
71#endif
72 _hypre_dof_idx(),
73 _hypre_num_nze(),
74 _hypre_col_idx(),
75 _hypre_matrix(nullptr),
76 _hypre_vec_def(nullptr),
77 _hypre_vec_cor(nullptr),
78 _hypre_parcsr(nullptr),
79 _hypre_par_def(nullptr),
80 _hypre_par_cor(nullptr)
81 {
82#ifndef FEAT_HAVE_MPI
83 (void)comm; // suppress unused parameter warnings
84#endif
85 // get dof offset and count
86 const HYPRE_Int dof_offset = HYPRE_Int(my_dof_offset);
87 const HYPRE_Int num_owned = HYPRE_Int(num_owned_dofs);
88
89 // set up HYPRE dof indices vector
90 _hypre_dof_idx.resize((std::size_t)num_owned);
91 for(HYPRE_Int i(0); i < num_owned; ++i)
92 _hypre_dof_idx[std::size_t(i)] = dof_offset + i;
93
94 // get lower and upper DOF bounds
95 const HYPRE_Int ilower = _hypre_dof_idx.front();
96 const HYPRE_Int iupper = _hypre_dof_idx.back();
97
98 // create HYPRE matrix
99 HYPRE_IJMatrixCreate(_comm, ilower, iupper, ilower, iupper, &_hypre_matrix);
100 HYPRE_IJMatrixSetObjectType(_hypre_matrix, HYPRE_PARCSR);
101 HYPRE_IJMatrixInitialize(_hypre_matrix);
102
103 // create HYPRE vectors
104 HYPRE_IJVectorCreate(_comm, ilower, iupper, &_hypre_vec_def);
105 HYPRE_IJVectorCreate(_comm, ilower, iupper, &_hypre_vec_cor);
106 HYPRE_IJVectorSetObjectType(_hypre_vec_def, HYPRE_PARCSR);
107 HYPRE_IJVectorSetObjectType(_hypre_vec_cor, HYPRE_PARCSR);
108 HYPRE_IJVectorInitialize(_hypre_vec_def);
109 HYPRE_IJVectorInitialize(_hypre_vec_cor);
110
111 // get matrix dimensions
112 const HYPRE_Int num_nze = HYPRE_Int(row_ptr[num_owned_dofs]);
113
114 // create HYPRE index vectors
115 _hypre_num_nze.resize((std::size_t)num_owned);
116 _hypre_col_idx.resize((std::size_t)num_nze);
117
118 // loop over all owned matrix rows
119 for(HYPRE_Int i(0); i < num_owned; ++i)
120 {
121 _hypre_num_nze[std::size_t(i)] = HYPRE_Int(row_ptr[i+1] - row_ptr[i]);
122 }
123 for(HYPRE_Int i(0); i < num_nze; ++i)
124 {
125 _hypre_col_idx[std::size_t(i)] = HYPRE_Int(col_idx[i]);
126 }
127
128 // vector for initial vector values
129 std::vector<HYPRE_Real> va((std::size_t)num_nze, 0.0);
130 std::vector<HYPRE_Real> vx((std::size_t)num_owned, 0.0);
131
132 // set matrix structure + dummy values
133 HYPRE_IJMatrixSetValues(_hypre_matrix, num_owned, _hypre_num_nze.data(), _hypre_dof_idx.data(),
134 _hypre_col_idx.data(), va.data());
135
136 // set vector structures + dummy values
137 HYPRE_IJVectorSetValues(_hypre_vec_def, num_owned, _hypre_dof_idx.data(), vx.data());
138 HYPRE_IJVectorSetValues(_hypre_vec_cor, num_owned, _hypre_dof_idx.data(), vx.data());
139
140 // assemble matrix and get ParCSR object pointer
141 HYPRE_IJMatrixAssemble(_hypre_matrix);
142 HYPRE_IJMatrixGetObject(_hypre_matrix, (void**) &_hypre_parcsr);
143
144 // assemble vectors and get Par object pointers
145 HYPRE_IJVectorAssemble(_hypre_vec_def);
146 HYPRE_IJVectorAssemble(_hypre_vec_cor);
147 HYPRE_IJVectorGetObject(_hypre_vec_def, (void **) &_hypre_par_def);
148 HYPRE_IJVectorGetObject(_hypre_vec_cor, (void **) &_hypre_par_cor);
149 }
150
151 ~Core()
152 {
153 HYPRE_IJVectorDestroy(_hypre_vec_cor);
154 HYPRE_IJVectorDestroy(_hypre_vec_def);
155 HYPRE_IJMatrixDestroy(_hypre_matrix);
156 }
157
158 void set_matrix_values(const double* vals)
159 {
160 HYPRE_IJMatrixSetValues(_hypre_matrix, HYPRE_Int(_hypre_dof_idx.size()), _hypre_num_nze.data(),
161 _hypre_dof_idx.data(), _hypre_col_idx.data(), vals);
162 }
163
164 void set_vec_cor_values(const double* vals)
165 {
166 if(vals != nullptr)
167 HYPRE_IJVectorSetValues(_hypre_vec_cor, HYPRE_Int(_hypre_dof_idx.size()), _hypre_dof_idx.data(), vals);
168 else
169 HYPRE_ParVectorSetConstantValues(_hypre_par_cor, HYPRE_Real(0));
170 }
171
172 void set_vec_def_values(const double* vals)
173 {
174 if(vals != nullptr)
175 HYPRE_IJVectorSetValues(_hypre_vec_def, HYPRE_Int(_hypre_dof_idx.size()), _hypre_dof_idx.data(), vals);
176 else
177 HYPRE_ParVectorSetConstantValues(_hypre_par_def, HYPRE_Real(0));
178 }
179
180 void get_vec_cor_values(double* vals) const
181 {
182 HYPRE_IJVectorGetValues(_hypre_vec_cor, HYPRE_Int(_hypre_dof_idx.size()), _hypre_dof_idx.data(), vals);
183 }
184
185 void get_vec_def_values(double* vals) const
186 {
187 HYPRE_IJVectorGetValues(_hypre_vec_cor, HYPRE_Int(_hypre_dof_idx.size()), _hypre_dof_idx.data(), vals);
188 }
189 }; // class Core
190
191 void* create_core(const void* comm, Index dof_offset, Index num_owned_dofs,
192 const unsigned int* row_ptr, const unsigned int* col_idx)
193 {
194 return new Core(comm, dof_offset, num_owned_dofs, row_ptr, col_idx);
195 }
196
197 void* create_core(const void* comm, Index dof_offset, Index num_owned_dofs,
198 const unsigned long* row_ptr, const unsigned long* col_idx)
199 {
200 return new Core(comm, dof_offset, num_owned_dofs, row_ptr, col_idx);
201 }
202
203 void* create_core(const void* comm, Index dof_offset, Index num_owned_dofs,
204 const unsigned long long* row_ptr, const unsigned long long* col_idx)
205 {
206 return new Core(comm, dof_offset, num_owned_dofs, row_ptr, col_idx);
207 }
208
209 void destroy_core(void* core)
210 {
211 delete reinterpret_cast<Core*>(core);
212 }
213
214 void set_matrix_values(void* core, const double* vals)
215 {
216 reinterpret_cast<Core*>(core)->set_matrix_values(vals);
217 }
218
219 void set_vec_cor_values(void* core, const double* vals)
220 {
221 reinterpret_cast<Core*>(core)->set_vec_cor_values(vals);
222 }
223
224 void set_vec_def_values(void* core, const double* vals)
225 {
226 reinterpret_cast<Core*>(core)->set_vec_def_values(vals);
227 }
228
229 void get_vec_cor_values(const void* core, double* vals)
230 {
231 reinterpret_cast<const Core*>(core)->get_vec_cor_values(vals);
232 }
233
234 void get_vec_def_values(const void* core, double* vals)
235 {
236 reinterpret_cast<const Core*>(core)->get_vec_def_values(vals);
237 }
238
239 /* *********************************************************************************************************** */
240 /* *********************************************************************************************************** */
241 /* *********************************************************************************************************** */
242
243 void* create_parasails(void* cr, int* iparam, double* dparam)
244 {
245 Core* core = reinterpret_cast<Core*>(cr);
246 HYPRE_Solver* solver = new HYPRE_Solver();
247
248 // create ParaSails preconditioner
249 HYPRE_ParaSailsCreate(core->_comm, solver);
250
251 // set ParaSails parameters
252 HYPRE_ParaSailsSetParams(*solver, dparam[0], iparam[0]);
253 HYPRE_ParaSailsSetFilter(*solver, dparam[1]);
254 HYPRE_ParaSailsSetSym(*solver, iparam[1]);
255
256 // setup ParaSails
257 // according to the documentation, the two vectors are ignored by this function
258 HYPRE_ParaSailsSetup(*solver, core->_hypre_parcsr, core->_hypre_par_def, core->_hypre_par_cor);
259
260 return solver;
261 }
262
263 void destroy_parasails(void* slv)
264 {
265 HYPRE_Solver* solver = reinterpret_cast<HYPRE_Solver*>(slv);
266 HYPRE_ParaSailsDestroy(*solver);
267 delete solver;
268 }
269
270 void solve_parasails(void* cr, void* slv)
271 {
272 Core* core = reinterpret_cast<Core*>(cr);
273 HYPRE_Solver* solver = reinterpret_cast<HYPRE_Solver*>(slv);
274
275 HYPRE_ParaSailsSolve(*solver, core->_hypre_parcsr, core->_hypre_par_def, core->_hypre_par_cor);
276 }
277
278 /* *********************************************************************************************************** */
279 /* *********************************************************************************************************** */
280 /* *********************************************************************************************************** */
281
282 void* create_euclid(void* cr, int* iparam, double* dparam)
283 {
284 Core* core = reinterpret_cast<Core*>(cr);
285 HYPRE_Solver* solver = new HYPRE_Solver();
286
287 // create Euclid preconditioner
288 HYPRE_EuclidCreate(core->_comm, solver);
289
290 // set Euclid parameters
291 HYPRE_EuclidSetLevel(*solver, iparam[0]);
292 HYPRE_EuclidSetSparseA(*solver, dparam[0]);
293
294 // setup Euclid
295 // according to the documentation, the two vectors are ignored by this function
296 HYPRE_EuclidSetup(*solver, core->_hypre_parcsr, core->_hypre_par_def, core->_hypre_par_cor);
297
298 return solver;
299 }
300
301 void destroy_euclid(void* slv)
302 {
303 HYPRE_Solver* solver = reinterpret_cast<HYPRE_Solver*>(slv);
304 HYPRE_EuclidDestroy(*solver);
305 delete solver;
306 }
307
308 void solve_euclid(void* cr, void* slv)
309 {
310 Core* core = reinterpret_cast<Core*>(cr);
311 HYPRE_Solver* solver = reinterpret_cast<HYPRE_Solver*>(slv);
312
313 HYPRE_EuclidSolve(*solver, core->_hypre_parcsr, core->_hypre_par_def, core->_hypre_par_cor);
314 }
315
316 /* *********************************************************************************************************** */
317 /* *********************************************************************************************************** */
318 /* *********************************************************************************************************** */
319
320 void* create_boomeramg(void* cr, int* iparam, double* dparam)
321 {
322 Core* core = reinterpret_cast<Core*>(cr);
323 HYPRE_Solver* solver = new HYPRE_Solver();
324
325 // create BoomerAMG preconditioner
326 HYPRE_BoomerAMGCreate(solver);
327
328 // set BoomerAMG parameters
329 HYPRE_BoomerAMGSetMaxIter(*solver, iparam[0]);
330 HYPRE_BoomerAMGSetTol(*solver, dparam[0]);
331
332 // setup BoomerAMG
333 // according to the documentation, the two vectors are ignored by this function
334 HYPRE_BoomerAMGSetup(*solver, core->_hypre_parcsr, core->_hypre_par_def, core->_hypre_par_cor);
335
336 return solver;
337 }
338
339 void destroy_boomeramg(void* slv)
340 {
341 HYPRE_Solver* solver = reinterpret_cast<HYPRE_Solver*>(slv);
342 HYPRE_BoomerAMGDestroy(*solver);
343 delete solver;
344 }
345
346 void solve_boomeramg(void* cr, void* slv)
347 {
348 Core* core = reinterpret_cast<Core*>(cr);
349 HYPRE_Solver* solver = reinterpret_cast<HYPRE_Solver*>(slv);
350
351 HYPRE_BoomerAMGSolve(*solver, core->_hypre_parcsr, core->_hypre_par_def, core->_hypre_par_cor);
352 }
353 } // namespace Hypre
354 } // namespace Solver
355} // namespace FEAT
356#else
357// insert dummy function to suppress linker warnings
358void dummy_hypre_function() {}
359#endif // FEAT_HAVE_HYPRE
base_header.hpp
FEAT Kernel base header.
FEAT
FEAT namespace.
Definition: adjactor.hpp:12