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;
49 const HYPRE_Int _dof_offset;
50 const HYPRE_Int _num_owned_dofs;
51 const HYPRE_Int _num_nonzeros;
53 std::vector<HYPRE_Int> _hypre_dof_idx;
55 std::vector<HYPRE_Int> _hypre_num_nze;
57 std::vector<HYPRE_Int> _hypre_row_ptr;
59 std::vector<HYPRE_Int> _hypre_col_idx;
61 std::vector<HYPRE_Real> _hypre_mat_val;
63 std::vector<HYPRE_Real> _hypre_def_val;
65 std::vector<HYPRE_Real> _hypre_cor_val;
66
67 // note: the following members are actually pointers
69 HYPRE_IJMatrix _hypre_matrix;
71 HYPRE_IJVector _hypre_vec_def, _hypre_vec_cor;
73 HYPRE_ParCSRMatrix _hypre_parcsr;
75 HYPRE_ParVector _hypre_par_def, _hypre_par_cor;
76
77 explicit Core(const void* comm, Index my_dof_offset, Index num_owned_dofs, Index num_nonzeros) :
78#ifdef FEAT_HAVE_MPI
79 _comm(*reinterpret_cast<const MPI_Comm*>(comm)),
80#endif
81 _dof_offset(HYPRE_Int(my_dof_offset)),
82 _num_owned_dofs(HYPRE_Int(num_owned_dofs)),
83 _num_nonzeros(HYPRE_Int(num_nonzeros)),
84 _hypre_dof_idx(),
85 _hypre_num_nze(),
86 _hypre_col_idx(),
87 _hypre_matrix(nullptr),
88 _hypre_vec_def(nullptr),
89 _hypre_vec_cor(nullptr),
90 _hypre_parcsr(nullptr),
91 _hypre_par_def(nullptr),
92 _hypre_par_cor(nullptr)
93 {
94#ifndef FEAT_HAVE_MPI
95 (void)comm; // suppress unused parameter warnings
96#endif
97 // set up HYPRE dof indices vector
98 _hypre_dof_idx.resize(std::size_t(_num_owned_dofs));
99 for(HYPRE_Int i(0); i < _num_owned_dofs; ++i)
100 _hypre_dof_idx[std::size_t(i)] = _dof_offset + i;
101
102 // get lower and upper DOF bounds
103 const HYPRE_Int ilower = _hypre_dof_idx.front();
104 const HYPRE_Int iupper = _hypre_dof_idx.back();
105
106 // create HYPRE matrix
107 HYPRE_IJMatrixCreate(_comm, ilower, iupper, ilower, iupper, &_hypre_matrix);
108 HYPRE_IJMatrixSetObjectType(_hypre_matrix, HYPRE_PARCSR);
109 HYPRE_IJMatrixInitialize(_hypre_matrix);
110
111 // create HYPRE vectors
112 HYPRE_IJVectorCreate(_comm, ilower, iupper, &_hypre_vec_def);
113 HYPRE_IJVectorCreate(_comm, ilower, iupper, &_hypre_vec_cor);
114 HYPRE_IJVectorSetObjectType(_hypre_vec_def, HYPRE_PARCSR);
115 HYPRE_IJVectorSetObjectType(_hypre_vec_cor, HYPRE_PARCSR);
116 HYPRE_IJVectorInitialize(_hypre_vec_def);
117 HYPRE_IJVectorInitialize(_hypre_vec_cor);
118
119 // create HYPRE index vectors
120 _hypre_num_nze.resize(std::size_t(_num_owned_dofs));
121 _hypre_row_ptr.resize(std::size_t(_num_owned_dofs+1));
122 _hypre_col_idx.resize(std::size_t(_num_nonzeros));
123 _hypre_mat_val.resize(std::size_t(_num_nonzeros));
124 _hypre_def_val.resize(std::size_t(_num_owned_dofs));
125 _hypre_cor_val.resize(std::size_t(_num_owned_dofs));
126 }
127
128 ~Core()
129 {
130 HYPRE_IJVectorDestroy(_hypre_vec_cor);
131 HYPRE_IJVectorDestroy(_hypre_vec_def);
132 HYPRE_IJMatrixDestroy(_hypre_matrix);
133 }
134
135 void upload_symbolic()
136 {
137 // compute row non-zeros from row pointer
138 std::size_t n = _hypre_num_nze.size();
139 for(std::size_t i(0); i < n; ++i)
140 _hypre_num_nze[i] = _hypre_row_ptr[i+1] - _hypre_row_ptr[i];
141
142 // set matrix structure + dummy values
143 HYPRE_IJMatrixSetValues(_hypre_matrix, _num_owned_dofs, _hypre_num_nze.data(), _hypre_dof_idx.data(),
144 _hypre_col_idx.data(), _hypre_mat_val.data());
145
146 // set vector structures + dummy values
147 HYPRE_IJVectorSetValues(_hypre_vec_def, _num_owned_dofs, _hypre_dof_idx.data(), _hypre_def_val.data());
148 HYPRE_IJVectorSetValues(_hypre_vec_cor, _num_owned_dofs, _hypre_dof_idx.data(), _hypre_cor_val.data());
149
150 // assemble matrix and get ParCSR object pointer
151 HYPRE_IJMatrixAssemble(_hypre_matrix);
152 HYPRE_IJMatrixGetObject(_hypre_matrix, (void**) &_hypre_parcsr);
153
154 // assemble vectors and get Par object pointers
155 HYPRE_IJVectorAssemble(_hypre_vec_def);
156 HYPRE_IJVectorAssemble(_hypre_vec_cor);
157 HYPRE_IJVectorGetObject(_hypre_vec_def, (void **) &_hypre_par_def);
158 HYPRE_IJVectorGetObject(_hypre_vec_cor, (void **) &_hypre_par_cor);
159 }
160
161 void upload_mat_val()
162 {
163 HYPRE_IJMatrixSetValues(_hypre_matrix, _num_owned_dofs, _hypre_num_nze.data(), _hypre_dof_idx.data(),
164 _hypre_col_idx.data(), _hypre_mat_val.data());
165 }
166
167 void upload_vec_def()
168 {
169 HYPRE_IJVectorSetValues(_hypre_vec_def, _num_owned_dofs, _hypre_dof_idx.data(), _hypre_def_val.data());
170 }
171
172 void download_vec_cor()
173 {
174 HYPRE_IJVectorGetValues(_hypre_vec_cor, _num_owned_dofs, _hypre_dof_idx.data(), _hypre_cor_val.data());
175 }
176
177 void format_vec_cor()
178 {
179 HYPRE_ParVectorSetConstantValues(_hypre_par_cor, HYPRE_Real(0));
180 }
181
182 HYPRE_Int* get_row_ptr()
183 {
184 return _hypre_row_ptr.data();
185 }
186
187 HYPRE_Int* get_col_idx()
188 {
189 return _hypre_col_idx.data();
190 }
191
192 HYPRE_Real* get_mat_val()
193 {
194 return _hypre_mat_val.data();
195 }
196
197 HYPRE_Real* get_vec_def()
198 {
199 return _hypre_def_val.data();
200 }
201
202 HYPRE_Real* get_vec_cor()
203 {
204 return _hypre_cor_val.data();
205 }
206 }; // class Core
207
208 void* create_core(const void* comm, Index dof_offset, Index num_owned_dofs, Index num_nonzeros)
209 {
210 return new Core(comm, dof_offset, num_owned_dofs, num_nonzeros);
211 }
212
213 void destroy_core(void* core)
214 {
215 delete reinterpret_cast<Core*>(core);
216 }
217
218 HYPRE_Int* get_row_ptr(void* core)
219 {
220 return reinterpret_cast<Core*>(core)->get_row_ptr();
221 }
222
223 HYPRE_Int* get_col_idx(void* core)
224 {
225 return reinterpret_cast<Core*>(core)->get_col_idx();
226 }
227
228 HYPRE_Real* get_mat_val(void* core)
229 {
230 return reinterpret_cast<Core*>(core)->get_mat_val();
231 }
232
233 HYPRE_Real* get_vec_def(void* core)
234 {
235 return reinterpret_cast<Core*>(core)->get_vec_def();
236 }
237
238 HYPRE_Real* get_vec_cor(void* core)
239 {
240 return reinterpret_cast<Core*>(core)->get_vec_cor();
241 }
242
243 void upload_symbolic(void* core)
244 {
245 reinterpret_cast<Core*>(core)->upload_symbolic();
246 }
247
248 void upload_mat_val(void* core)
249 {
250 reinterpret_cast<Core*>(core)->upload_mat_val();
251 }
252
253 void upload_vec_def(void* core)
254 {
255 reinterpret_cast<Core*>(core)->upload_vec_def();
256 }
257
258 void download_vec_cor(void* core)
259 {
260 reinterpret_cast<Core*>(core)->download_vec_cor();
261 }
262
263 void format_vec_cor(void* core)
264 {
265 reinterpret_cast<Core*>(core)->format_vec_cor();
266 }
267
268 /* *********************************************************************************************************** */
269 /* *********************************************************************************************************** */
270 /* *********************************************************************************************************** */
271
272 void* create_parasails(void* cr, int* iparam, double* dparam)
273 {
274 Core* core = reinterpret_cast<Core*>(cr);
275 HYPRE_Solver* solver = new HYPRE_Solver();
276
277 // create ParaSails preconditioner
278 HYPRE_ParaSailsCreate(core->_comm, solver);
279
280 // set ParaSails parameters
281 HYPRE_ParaSailsSetParams(*solver, dparam[0], iparam[0]);
282 HYPRE_ParaSailsSetFilter(*solver, dparam[1]);
283 HYPRE_ParaSailsSetSym(*solver, iparam[1]);
284
285 // setup ParaSails
286 // according to the documentation, the two vectors are ignored by this function
287 HYPRE_ParaSailsSetup(*solver, core->_hypre_parcsr, core->_hypre_par_def, core->_hypre_par_cor);
288
289 return solver;
290 }
291
292 void destroy_parasails(void* slv)
293 {
294 HYPRE_Solver* solver = reinterpret_cast<HYPRE_Solver*>(slv);
295 HYPRE_ParaSailsDestroy(*solver);
296 delete solver;
297 }
298
299 void solve_parasails(void* cr, void* slv)
300 {
301 Core* core = reinterpret_cast<Core*>(cr);
302 HYPRE_Solver* solver = reinterpret_cast<HYPRE_Solver*>(slv);
303
304 HYPRE_ParaSailsSolve(*solver, core->_hypre_parcsr, core->_hypre_par_def, core->_hypre_par_cor);
305 }
306
307 /* *********************************************************************************************************** */
308 /* *********************************************************************************************************** */
309 /* *********************************************************************************************************** */
310
311 void* create_euclid(void* cr, int* iparam, double* dparam)
312 {
313 Core* core = reinterpret_cast<Core*>(cr);
314 HYPRE_Solver* solver = new HYPRE_Solver();
315
316 // create Euclid preconditioner
317 HYPRE_EuclidCreate(core->_comm, solver);
318
319 // set Euclid parameters
320 HYPRE_EuclidSetLevel(*solver, iparam[0]);
321 HYPRE_EuclidSetSparseA(*solver, dparam[0]);
322
323 // setup Euclid
324 // according to the documentation, the two vectors are ignored by this function
325 HYPRE_EuclidSetup(*solver, core->_hypre_parcsr, core->_hypre_par_def, core->_hypre_par_cor);
326
327 return solver;
328 }
329
330 void destroy_euclid(void* slv)
331 {
332 HYPRE_Solver* solver = reinterpret_cast<HYPRE_Solver*>(slv);
333 HYPRE_EuclidDestroy(*solver);
334 delete solver;
335 }
336
337 void solve_euclid(void* cr, void* slv)
338 {
339 Core* core = reinterpret_cast<Core*>(cr);
340 HYPRE_Solver* solver = reinterpret_cast<HYPRE_Solver*>(slv);
341
342 HYPRE_EuclidSolve(*solver, core->_hypre_parcsr, core->_hypre_par_def, core->_hypre_par_cor);
343 }
344
345 /* *********************************************************************************************************** */
346 /* *********************************************************************************************************** */
347 /* *********************************************************************************************************** */
348
349 void* create_boomeramg(void* cr, int* iparam, double* dparam)
350 {
351 Core* core = reinterpret_cast<Core*>(cr);
352 HYPRE_Solver* solver = new HYPRE_Solver();
353
354 // create BoomerAMG preconditioner
355 HYPRE_BoomerAMGCreate(solver);
356
357 // set BoomerAMG parameters
358 HYPRE_BoomerAMGSetMaxIter(*solver, iparam[0]);
359 HYPRE_BoomerAMGSetTol(*solver, dparam[0]);
360
361 // setup BoomerAMG
362 // according to the documentation, the two vectors are ignored by this function
363 HYPRE_BoomerAMGSetup(*solver, core->_hypre_parcsr, core->_hypre_par_def, core->_hypre_par_cor);
364
365 return solver;
366 }
367
368 void destroy_boomeramg(void* slv)
369 {
370 HYPRE_Solver* solver = reinterpret_cast<HYPRE_Solver*>(slv);
371 HYPRE_BoomerAMGDestroy(*solver);
372 delete solver;
373 }
374
375 void solve_boomeramg(void* cr, void* slv)
376 {
377 Core* core = reinterpret_cast<Core*>(cr);
378 HYPRE_Solver* solver = reinterpret_cast<HYPRE_Solver*>(slv);
379
380 HYPRE_BoomerAMGSolve(*solver, core->_hypre_parcsr, core->_hypre_par_def, core->_hypre_par_cor);
381 }
382 } // namespace Hypre
383 } // namespace Solver
384} // namespace FEAT
385#else
386// insert dummy function to suppress linker warnings
387void dummy_hypre_function() {}
388#endif // FEAT_HAVE_HYPRE
base_header.hpp
FEAT Kernel base header.
FEAT
FEAT namespace.
Definition: adjactor.hpp:12