Arcane  4.1.12.0
Developer documentation
Loading...
Searching...
No Matches
AlephHypre.cc
1// -*- tab-width: 2; indent-tabs-mode: nil; coding: utf-8-with-signature -*-
2//-----------------------------------------------------------------------------
3// Copyright 2000-2026 CEA (www.cea.fr) IFPEN (www.ifpenergiesnouvelles.com)
4// See the top-level COPYRIGHT file for details.
5// SPDX-License-Identifier: Apache-2.0
6//-----------------------------------------------------------------------------
7/*---------------------------------------------------------------------------*/
8/* AlephHypre.cc (C) 2000-2025 */
9/* */
10/* Hypre implementation of Aleph. */
11/*---------------------------------------------------------------------------*/
12/*---------------------------------------------------------------------------*/
13
14#define HAVE_MPI
15#define MPI_COMM_SUB (*(MPI_Comm*)(m_kernel->subParallelMng(m_index)->getMPICommunicator()))
16#define OMPI_SKIP_MPICXX
17#ifndef MPICH_SKIP_MPICXX
18#define MPICH_SKIP_MPICXX
19#endif
20#include <HYPRE.h>
21#include <HYPRE_utilities.h>
22#include <HYPRE_IJ_mv.h>
23#include <HYPRE_parcsr_mv.h>
24#include <HYPRE_parcsr_ls.h>
25#include <_hypre_parcsr_mv.h>
26
27#if HYPRE_RELEASE_NUMBER >= 30000
28#include <_hypre_krylov.h>
29#else
30#include <krylov.h>
31#endif
32
33#ifndef ItacRegion
34#define ItacRegion(a, x)
35#endif
36
37#include "arcane/aleph/AlephArcane.h"
38
39// The HYPRE_BigInt type only exists starting from Hypre 2.16.0
40#if HYPRE_RELEASE_NUMBER < 21600
41using HYPRE_BigInt = HYPRE_Int;
42#endif
43
44/*---------------------------------------------------------------------------*/
45/*---------------------------------------------------------------------------*/
46
47namespace Arcane
48{
49
50/*---------------------------------------------------------------------------*/
51/*---------------------------------------------------------------------------*/
52
53/*
54 * NOTE: Starting from Hypre version 2.14 (maybe a little earlier),
55 * hypre_TAlloc() and hypre_CAlloc() take a 3rd argument specifying
56 * which peripheral the memory is allocated on (GPU or CPU). There is no
57 * simple way to know the Hypre version from the .h files, but since
58 * HYPRE_MEMORY_DEVICE and HYPRE_MEMORY_HOST are macros, we can test
59 * their existence to determine whether to call hypre_TAlloc() and
60 * hypre_CAlloc() with 2 or 3 arguments.
61 */
62namespace
63{
64 inline void
65 check(const char* hypre_func, HYPRE_Int error_code)
66 {
67 if (error_code == 0)
68 return;
69 char buf[8192];
70 HYPRE_DescribeError(error_code, buf);
71 cout << "\nXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX"
72 << "\nHYPRE ERROR in function "
73 << hypre_func
74 << "\nError_code=" << error_code
75 << "\nMessage=" << buf
76 << "\nXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX"
77 << "\n"
78 << std::flush;
79 throw Exception("HYPRE Check", hypre_func);
80 }
81
82 template <typename T>
83 inline T*
84 _allocHypre(Integer size)
85 {
86 size_t s = sizeof(T) * size;
87 return reinterpret_cast<T*>(hypre_TAlloc(char, s, HYPRE_MEMORY_HOST));
88 }
89
90 template <typename T>
91 inline T*
92 _callocHypre(Integer size)
93 {
94 size_t s = sizeof(T) * size;
95 return reinterpret_cast<T*>(hypre_CTAlloc(char, s, HYPRE_MEMORY_HOST));
96 }
97
98 inline void
99 hypreCheck(const char* hypre_func, HYPRE_Int error_code)
100 {
101 check(hypre_func, error_code);
102 HYPRE_Int r = HYPRE_GetError();
103 if (r != 0)
104 cout << "HYPRE GET ERROR r=" << r
105 << " error_code=" << error_code << " func=" << hypre_func << '\n';
106 }
107
108} // namespace
109
110/******************************************************************************
111 AlephVectorHypre
112 *****************************************************************************/
113class AlephVectorHypre
114: public IAlephVector
115{
116 public:
117
118 AlephVectorHypre(ITraceMng* tm, AlephKernel* kernel, Integer index)
119 : IAlephVector(tm, kernel, index)
120 , jSize(0)
121 , jUpper(0)
122 , jLower(-1)
123 {
124 debug() << "[AlephVectorHypre::AlephVectorHypre] new SolverVectorHypre";
125 }
126 ~AlephVectorHypre()
127 {
128 if (m_hypre_ijvector)
129 HYPRE_IJVectorDestroy(m_hypre_ijvector);
130 }
131
132 public:
133
134 /******************************************************************************
135 * The Create() routine creates an empty vector object that lives on the comm communicator. This is
136 * a collective call, with each process passing its own index extents, jLower and jupper. The names
137 * of these extent parameters begin with a j because we typically think of matrix-vector multiplies
138 * as the fundamental operation involving both matrices and vectors. For matrix-vector multiplies,
139 * the vector partitioning should match the column partitioning of the matrix (which also uses the j
140 * notation). For linear system solves, these extents will typically match the row partitioning of the
141 * matrix as well.
142 *****************************************************************************/
143 void AlephVectorCreate(void)
144 {
145 debug() << "[AlephVectorHypre::AlephVectorCreate] HYPRE VectorCreate";
146 void* object;
147
148 for (int iCpu = 0; iCpu < m_kernel->size(); ++iCpu) {
149 if (m_kernel->rank() != m_kernel->solverRanks(m_index)[iCpu])
150 continue;
151 debug() << "[AlephVectorHypre::AlephVectorCreate] adding contibution of core #" << iCpu;
152 if (jLower == -1)
153 jLower = m_kernel->topology()->gathered_nb_row(iCpu);
154 jUpper = m_kernel->topology()->gathered_nb_row(iCpu + 1) - 1;
155 }
156
157 debug() << "[AlephVectorHypre::AlephVectorCreate] jLower=" << jLower << ", jupper=" << jUpper;
158 // Update the local buffer size for later max norm calculation, for example
159 jSize = jUpper - jLower + 1;
160
161 hypreCheck("IJVectorCreate",
162 HYPRE_IJVectorCreate(MPI_COMM_SUB,
163 jLower,
164 jUpper,
165 &m_hypre_ijvector));
166
167 debug() << "[AlephVectorHypre::AlephVectorCreate] HYPRE IJVectorSetObjectType";
168 hypreCheck("IJVectorSetObjectType", HYPRE_IJVectorSetObjectType(m_hypre_ijvector, HYPRE_PARCSR));
169
170 debug() << "[AlephVectorHypre::AlephVectorCreate] HYPRE IJVectorInitialize";
171 hypreCheck("HYPRE_IJVectorInitialize", HYPRE_IJVectorInitialize(m_hypre_ijvector));
172
173 HYPRE_IJVectorGetObject(m_hypre_ijvector, &object);
174 m_hypre_parvector = (HYPRE_ParVector)object;
175 debug() << "[AlephVectorHypre::AlephVectorCreate] done";
176 }
177
178 /******************************************************************************
179 *****************************************************************************/
180 void AlephVectorSet(const double* bfr_val, const AlephInt* bfr_idx, Integer size)
181 {
182 debug() << "[AlephVectorHypre::AlephVectorSet] size=" << size;
183 hypreCheck("IJVectorSetValues", HYPRE_IJVectorSetValues(m_hypre_ijvector, size, bfr_idx, bfr_val));
184 }
185
186 /******************************************************************************
187 *****************************************************************************/
188 int AlephVectorAssemble(void)
189 {
190 debug() << "[AlephVectorHypre::AlephVectorAssemble]";
191 hypreCheck("IJVectorAssemble", HYPRE_IJVectorAssemble(m_hypre_ijvector));
192 return 0;
193 }
194
195 /******************************************************************************
196 *****************************************************************************/
197 void AlephVectorGet(double* bfr_val, const AlephInt* bfr_idx, Integer size)
198 {
199 //HYPRE_Int* hypre_bfr_idx = static
200 debug() << "[AlephVectorHypre::AlephVectorGet] size=" << size;
201 hypreCheck("HYPRE_IJVectorGetValues", HYPRE_IJVectorGetValues(m_hypre_ijvector, size, bfr_idx, bfr_val));
202 }
203
204 /******************************************************************************
205 * norm_max
206 *****************************************************************************/
207 Real norm_max()
208 {
209 Real normInf = 0.0;
210 UniqueArray<HYPRE_BigInt> bfr_idx(jSize);
211 UniqueArray<double> bfr_val(jSize);
212
213 for (HYPRE_Int i = 0; i < jSize; ++i)
214 bfr_idx[i] = jLower + i;
215
216 hypreCheck("HYPRE_IJVectorGetValues", HYPRE_IJVectorGetValues(m_hypre_ijvector, jSize, bfr_idx.data(), bfr_val.data()));
217 for (HYPRE_Int i = 0; i < jSize; ++i) {
218 const Real abs_val = math::abs(bfr_val[i]);
219 if (abs_val > normInf)
220 normInf = abs_val;
221 }
222 normInf = m_kernel->subParallelMng(m_index)->reduce(Parallel::ReduceMax, normInf);
223 return normInf;
224 }
225
226 /******************************************************************************
227 *****************************************************************************/
228 void writeToFile(const String filename)
229 {
230 String filename_idx = filename; // + "_" + (int)m_kernel->subDomain()->commonVariables().globalIteration();
231 debug() << "[AlephVectorHypre::writeToFile]";
232 hypreCheck("HYPRE_IJVectorPrint",
233 HYPRE_IJVectorPrint(m_hypre_ijvector, filename_idx.localstr()));
234 }
235
236 public:
237
238 HYPRE_IJVector m_hypre_ijvector = nullptr;
239 HYPRE_ParVector m_hypre_parvector = nullptr;
240 HYPRE_Int jSize;
241 HYPRE_Int jUpper;
242 HYPRE_Int jLower;
243};
244
245/******************************************************************************
246 AlephMatrixHypre
247*****************************************************************************/
248class AlephMatrixHypre
249: public IAlephMatrix
250{
251 public:
252
253 /******************************************************************************
254 AlephMatrixHypre
255 *****************************************************************************/
256 AlephMatrixHypre(ITraceMng* tm, AlephKernel* kernel, Integer index)
257 : IAlephMatrix(tm, kernel, index)
258 , m_hypre_ijmatrix(0)
259 {
260 debug() << "[AlephMatrixHypre] new AlephMatrixHypre";
261 }
262
263 ~AlephMatrixHypre()
264 {
265 debug() << "[~AlephMatrixHypre]";
266 if (m_hypre_ijmatrix)
267 HYPRE_IJMatrixDestroy(m_hypre_ijmatrix);
268 }
269
270 public:
271
272 /******************************************************************************
273 * Each submatrix Ap is "owned" by a single process and its first and last row numbers are
274 * given by the global indices ilower and iupper in the Create() call below.
275 *******************************************************************************
276 * The Create() routine creates an empty matrix object that lives on the comm communicator. This
277 * is a collective call (i.e., must be called on all processes from a common synchronization point),
278 * with each process passing its own row extents, ilower and iupper. The row partitioning must be
279 * contiguous, i.e., iupper for process i must equal ilower-1 for process i+1. Note that this allows
280 * matrices to have 0- or 1-based indexing. The parameters jlower and jupper define a column
281 * partitioning, and should match ilower and iupper when solving square linear systems.
282 *****************************************************************************/
283 void AlephMatrixCreate(void)
284 {
285 debug() << "[AlephMatrixHypre::AlephMatrixCreate] HYPRE MatrixCreate idx:" << m_index;
286 void* object;
287 AlephInt ilower = -1;
288 AlephInt iupper = 0;
289 for (int iCpu = 0; iCpu < m_kernel->size(); ++iCpu) {
290 if (m_kernel->rank() != m_kernel->solverRanks(m_index)[iCpu])
291 continue;
292 if (ilower == -1)
293 ilower = m_kernel->topology()->gathered_nb_row(iCpu);
294 iupper = m_kernel->topology()->gathered_nb_row(iCpu + 1) - 1;
295 }
296 debug() << "[AlephMatrixHypre::AlephMatrixCreate] ilower=" << ilower << ", iupper=" << iupper;
297
298 AlephInt jlower = ilower; //0;
299 AlephInt jupper = iupper; //m_kernel->topology()->gathered_nb_row(m_kernel->size())-1;
300 debug() << "[AlephMatrixHypre::AlephMatrixCreate] jlower=" << jlower << ", jupper=" << jupper;
301
302 hypreCheck("HYPRE_IJMatrixCreate",
303 HYPRE_IJMatrixCreate(MPI_COMM_SUB,
304 ilower, iupper,
305 jlower, jupper,
306 &m_hypre_ijmatrix));
307
308 debug() << "[AlephMatrixHypre::AlephMatrixCreate] HYPRE IJMatrixSetObjectType";
309 HYPRE_IJMatrixSetObjectType(m_hypre_ijmatrix, HYPRE_PARCSR);
310 debug() << "[AlephMatrixHypre::AlephMatrixCreate] HYPRE IJMatrixSetRowSizes";
311 HYPRE_IJMatrixSetRowSizes(m_hypre_ijmatrix, m_kernel->topology()->gathered_nb_row_elements().data());
312 debug() << "[AlephMatrixHypre::AlephMatrixCreate] HYPRE IJMatrixInitialize";
313 HYPRE_IJMatrixInitialize(m_hypre_ijmatrix);
314 HYPRE_IJMatrixGetObject(m_hypre_ijmatrix, &object);
315 m_hypre_parmatrix = (HYPRE_ParCSRMatrix)object;
316 }
317
318 /******************************************************************************
319 *****************************************************************************/
320 void AlephMatrixSetFilled(bool) {}
321
322 /******************************************************************************
323 *****************************************************************************/
324 int AlephMatrixAssemble(void)
325 {
326 debug() << "[AlephMatrixHypre::AlephMatrixAssemble]";
327 hypreCheck("HYPRE_IJMatrixAssemble",
328 HYPRE_IJMatrixAssemble(m_hypre_ijmatrix));
329 return 0;
330 }
331
332 /******************************************************************************
333 *****************************************************************************/
334 void AlephMatrixFill(int size, HYPRE_Int* rows, HYPRE_Int* cols, double* values)
335 {
336 debug() << "[AlephMatrixHypre::AlephMatrixFill] size=" << size;
337 HYPRE_Int rtn = 0;
338 HYPRE_Int col[1] = { 1 };
339 for (int i = 0; i < size; i++) {
340 rtn += HYPRE_IJMatrixSetValues(m_hypre_ijmatrix, 1, col, &rows[i], &cols[i], &values[i]);
341 }
342 hypreCheck("HYPRE_IJMatrixSetValues", rtn);
343 //HYPRE_IJMatrixSetValues(m_hypre_ijmatrix, nrows, ncols, rows, cols, values);
344 debug() << "[AlephMatrixHypre::AlephMatrixFill] done";
345 }
346
347 /******************************************************************************
348 * isAlreadySolved
349 *****************************************************************************/
350 bool isAlreadySolved(AlephVectorHypre* x,
351 AlephVectorHypre* b,
352 AlephVectorHypre* tmp,
353 Real* residual_norm,
354 AlephParams* params)
355 {
356 HYPRE_ClearAllErrors();
357 const bool convergence_analyse = params->convergenceAnalyse();
358
359 // test the right-hand side of the linear system
360 const Real res0 = b->norm_max();
361
362 if (convergence_analyse)
363 info() << "convergence analysis: max norm of the right-hand side res0: " << res0;
364
365 const Real considered_as_null = params->minRHSNorm();
366 if (res0 < considered_as_null) {
367 HYPRE_ParVectorSetConstantValues(x->m_hypre_parvector, 0.0);
368 residual_norm[0] = res0;
369 if (convergence_analyse)
370 info() << "convergence analysis: the right-hand side of the linear system is less than: " << considered_as_null;
371 return true;
372 }
373
374 if (params->xoUser()) {
375 // we test if b is already a solution to the system within epsilon tolerance
376 //matrix->vectorProduct(b, tmp_vector); tmp_vector->sub(x);
377 //M->vector_multiply(*tmp,*x); // tmp=A*x
378 //tmp->substract(*tmp,*b); // tmp=A*x-b
379
380 // X= alpha* M.B + beta * X (read from HYPRE sources)
381 HYPRE_ParCSRMatrixMatvec(1.0, m_hypre_parmatrix, x->m_hypre_parvector, 0., tmp->m_hypre_parvector);
382 HYPRE_ParVectorAxpy(-1.0, b->m_hypre_parvector, tmp->m_hypre_parvector);
383
384 const Real residu = tmp->norm_max();
385 //info() << "[IAlephHypre::isAlreadySolved] residual="<<residu;
386
387 if (residu < considered_as_null) {
388 if (convergence_analyse) {
389 info() << "convergence analysis: |Ax0-b| is less than " << considered_as_null;
390 info() << "convergence analysis: x0 is already a solution to the system.";
391 }
392 residual_norm[0] = residu;
393 return true;
394 }
395
396 const Real relative_error = residu / res0;
397 if (convergence_analyse)
398 info() << "convergence analysis: initial residual : " << residu
399 << " --- initial relative residual (residu/res0) : " << residu / res0;
400
401 if (relative_error < (params->epsilon())) {
402 if (convergence_analyse)
403 info() << "convergence analysis: X is already a solution to the system";
404 residual_norm[0] = residu;
405 return true;
406 }
407 }
408 return false;
409 }
410
411 /******************************************************************************
412 *****************************************************************************/
413 int AlephMatrixSolve(AlephVector* x,
414 AlephVector* b,
415 AlephVector* t,
416 Integer& nb_iteration,
417 Real* residual_norm,
418 AlephParams* solver_param)
419 {
420 solver_param->setAmgCoarseningMethod(TypesSolver::AMG_COARSENING_AUTO);
421 const String func_name("SolverMatrixHypre::solve");
422 void* object;
423 int ierr = 0;
424
425 auto* ximpl = ARCANE_CHECK_POINTER(dynamic_cast<AlephVectorHypre*>(x->implementation()));
426 auto* bimpl = ARCANE_CHECK_POINTER(dynamic_cast<AlephVectorHypre*>(b->implementation()));
427
428 HYPRE_IJVector solution = ximpl->m_hypre_ijvector;
429 HYPRE_IJVector RHS = bimpl->m_hypre_ijvector;
430 //HYPRE_IJVector tmp = (dynamic_cast<AlephVectorHypre*> (t->implementation()))->m_hypre_ijvector;
431
432 HYPRE_IJMatrixGetObject(m_hypre_ijmatrix, &object);
433 HYPRE_ParCSRMatrix M = (HYPRE_ParCSRMatrix)object;
434 HYPRE_IJVectorGetObject(solution, &object);
435 HYPRE_ParVector X = (HYPRE_ParVector)object;
436 HYPRE_IJVectorGetObject(RHS, &object);
437 HYPRE_ParVector B = (HYPRE_ParVector)object;
438 //HYPRE_IJVectorGetObject(tmp,&object);
439 //HYPRE_ParVector T = (HYPRE_ParVector)object;
440
441 auto* ximpl2 = ARCANE_CHECK_POINTER(dynamic_cast<AlephVectorHypre*>(x->implementation()));
442 auto* bimpl2 = ARCANE_CHECK_POINTER(dynamic_cast<AlephVectorHypre*>(b->implementation()));
443 auto* timpl2 = ARCANE_CHECK_POINTER(dynamic_cast<AlephVectorHypre*>(t->implementation()));
444 if (isAlreadySolved(ximpl2, bimpl2, timpl2, residual_norm, solver_param)) {
445 ItacRegion(isAlreadySolved, AlephMatrixHypre);
446 debug() << "[AlephMatrixHypre::AlephMatrixSolve] isAlreadySolved !";
447 nb_iteration = 0;
448 return 0;
449 }
450
451 TypesSolver::ePreconditionerMethod preconditioner_method = solver_param->precond();
452 // TypesSolver::ePreconditionerMethod preconditioner_method = TypesSolver::NONE;
453 TypesSolver::eSolverMethod solver_method = solver_param->method();
454
455 // declaration and initialization of the solver
456 HYPRE_Solver solver = 0;
457
458 switch (solver_method) {
459 case TypesSolver::PCG:
460 initSolverPCG(solver_param, solver);
461 break;
462 case TypesSolver::BiCGStab:
463 initSolverBiCGStab(solver_param, solver);
464 break;
465 case TypesSolver::GMRES:
466 initSolverGMRES(solver_param, solver);
467 break;
468 default:
469 throw ArgumentException(func_name, "unknown solver");
470 }
471
472 // declaration and initialization of the preconditioner
473 HYPRE_Solver precond = 0;
474
475 switch (preconditioner_method) {
476 case TypesSolver::NONE:
477 break;
478 case TypesSolver::DIAGONAL:
479 setDiagonalPreconditioner(solver_method, solver, precond);
480 break;
481 case TypesSolver::ILU:
482 setILUPreconditioner(solver_method, solver, precond);
483 break;
484 case TypesSolver::SPAIstat:
485 setSpaiStatPreconditioner(solver_method, solver, solver_param, precond);
486 break;
487 case TypesSolver::AMG:
488 setAMGPreconditioner(solver_method, solver, solver_param, precond);
489 break;
490 case TypesSolver::AINV:
491 throw ArgumentException(func_name, "AINV preconditioning unavailable");
492 case TypesSolver::SPAIdyn:
493 throw ArgumentException(func_name, "SPAIdyn preconditioning unavailable");
494 case TypesSolver::ILUp:
495 throw ArgumentException(func_name, "ILUp preconditioning unavailable");
496 case TypesSolver::IC:
497 throw ArgumentException(func_name, "IC preconditioning unavailable");
498 case TypesSolver::POLY:
499 throw ArgumentException(func_name, "POLY preconditioning unavailable");
500 default:
501 throw ArgumentException(func_name, "unknown preconditioner");
502 }
503
504 // solving the algebraic system
505 HYPRE_Int iteration = 0;
506 double residue = 0.0;
507
508 switch (solver_method) {
509 case TypesSolver::PCG:
510 ierr = solvePCG(solver_param, solver, M, B, X, iteration, residue);
511 break;
512 case TypesSolver::BiCGStab:
513 ierr = solveBiCGStab(solver, M, B, X, iteration, residue);
514 break;
515 case TypesSolver::GMRES:
516 ierr = solveGMRES(solver, M, B, X, iteration, residue);
517 break;
518 default:
519 ierr = -3;
520 return ierr;
521 }
522 nb_iteration = static_cast<Integer>(iteration);
523 residual_norm[0] = static_cast<Real>(residue);
524
525 /* for(int i=0;i<8;++i){
526 int idx[1];
527 double valx[1]={-1.};
528 //double valb[1]={-1.};
529 idx[0]=i;
530 HYPRE_IJVectorGetValues(solution, 1, idx, valx);
531 debug()<<"[AlephMatrixHypre::AlephMatrixSolve] X["<<i<<"]="<<valx[0];
532 //(static_cast<AlephVectorHypre*> (x->implementation()))->AlephVectorGet(valx,idx,1);
533 //debug()<<"[AlephMatrixHypre::AlephMatrixSolve] x["<<i<<"]="<<valx[0];
534 //HYPRE_IJVectorGetValues(RHS, 1, idx, valb);
535 //debug()<<"[AlephMatrixHypre::AlephMatrixSolve] B["<<i<<"]="<<valb[0];
536 }
537*/
538 switch (preconditioner_method) {
539 case TypesSolver::NONE:
540 break;
541 case TypesSolver::DIAGONAL:
542 break;
543 case TypesSolver::ILU:
544 HYPRE_ParCSRPilutDestroy(precond);
545 break;
546 case TypesSolver::SPAIstat:
547 HYPRE_ParCSRParaSailsDestroy(precond);
548 break;
549 case TypesSolver::AMG:
550 HYPRE_BoomerAMGDestroy(precond);
551 break;
552 default:
553 throw ArgumentException(func_name, "unknown preconditioner");
554 }
555
556 if (iteration == solver_param->maxIter() && solver_param->stopErrorStrategy()) {
557 info() << "\n============================================================";
558 info() << "\nThis error is returned after " << iteration << "\n";
559 info() << "\nMaximum number of solver iterations reached.";
560 info() << "\nIt is possible to ask the code not to consider this error.";
561 info() << "\nSee the dataset documentation regarding the solver service.";
562 info() << "\n======================================================";
563 throw Exception("AlephMatrixHypre::Solve", "Maximum number of solver iterations reached");
564 }
565 return ierr;
566 }
567
568 /******************************************************************************
569 *****************************************************************************/
570 void writeToFile(const String filename)
571 {
572 HYPRE_IJMatrixPrint(m_hypre_ijmatrix, filename.localstr());
573 }
574
575 /******************************************************************************
576 *****************************************************************************/
577 void initSolverPCG(const AlephParams* solver_param, HYPRE_Solver& solver)
578 {
579 const String func_name = "SolverMatrixHypre::initSolverPCG";
580 double epsilon = solver_param->epsilon();
581 int max_it = solver_param->maxIter();
582 int output_level = solver_param->getOutputLevel();
583
584 HYPRE_ParCSRPCGCreate(MPI_COMM_SUB, &solver);
585 HYPRE_ParCSRPCGSetMaxIter(solver, max_it);
586 HYPRE_ParCSRPCGSetTol(solver, epsilon);
587 HYPRE_ParCSRPCGSetTwoNorm(solver, 1);
588 HYPRE_ParCSRPCGSetPrintLevel(solver, output_level);
589 }
590
591 /******************************************************************************
592 *****************************************************************************/
593 void initSolverBiCGStab(const AlephParams* solver_param, HYPRE_Solver& solver)
594 {
595 const String func_name = "SolverMatrixHypre::initSolverBiCGStab";
596 double epsilon = solver_param->epsilon();
597 int max_it = solver_param->maxIter();
598 int output_level = solver_param->getOutputLevel();
599
600 HYPRE_ParCSRBiCGSTABCreate(MPI_COMM_SUB, &solver);
601 HYPRE_ParCSRBiCGSTABSetMaxIter(solver, max_it);
602 HYPRE_ParCSRBiCGSTABSetTol(solver, epsilon);
603 HYPRE_ParCSRBiCGSTABSetPrintLevel(solver, output_level);
604 }
605
606 /******************************************************************************
607 *****************************************************************************/
608 void initSolverGMRES(const AlephParams* solver_param, HYPRE_Solver& solver)
609 {
610 const String func_name = "SolverMatrixHypre::initSolverGMRES";
611 double epsilon = solver_param->epsilon();
612 int max_it = solver_param->maxIter();
613 int output_level = solver_param->getOutputLevel();
614
615 HYPRE_ParCSRGMRESCreate(MPI_COMM_SUB, &solver);
616 const int krylov_dim = 20; // dimension Krylov space for GMRES
617 HYPRE_ParCSRGMRESSetKDim(solver, krylov_dim);
618 HYPRE_ParCSRGMRESSetMaxIter(solver, max_it);
619 HYPRE_ParCSRGMRESSetTol(solver, epsilon);
620 HYPRE_ParCSRGMRESSetPrintLevel(solver, output_level);
621 }
622
623 /******************************************************************************
624 *****************************************************************************/
625 void setDiagonalPreconditioner(const TypesSolver::eSolverMethod solver_method,
626 const HYPRE_Solver& solver,
627 HYPRE_Solver& precond)
628 {
629 const String func_name = "SolverMatrixHypre::setDiagonalPreconditioner";
630 switch (solver_method) {
631 case TypesSolver::PCG:
632 HYPRE_ParCSRPCGSetPrecond(solver,
633 HYPRE_ParCSRDiagScale,
634 HYPRE_ParCSRDiagScaleSetup,
635 precond);
636 break;
637 case TypesSolver::BiCGStab:
638 HYPRE_ParCSRBiCGSTABSetPrecond(solver,
639 HYPRE_ParCSRDiagScale,
640 HYPRE_ParCSRDiagScaleSetup,
641 precond);
642 break;
643 case TypesSolver::GMRES:
644 HYPRE_ParCSRGMRESSetPrecond(solver,
645 HYPRE_ParCSRDiagScale,
646 HYPRE_ParCSRDiagScaleSetup,
647 precond);
648 break;
649 default:
650 throw ArgumentException(func_name, "unknown solver for 'Diagonal' preconditioner");
651 }
652 }
653
654 /******************************************************************************
655 *****************************************************************************/
656 void setILUPreconditioner(const TypesSolver::eSolverMethod solver_method,
657 const HYPRE_Solver& solver,
658 HYPRE_Solver& precond)
659 {
660 const String func_name = "SolverMatrixHypre::setILUPreconditioner";
661 switch (solver_method) {
662 case TypesSolver::PCG:
663 throw ArgumentException(func_name, "PCG solver unavailable with 'ILU' preconditioner");
664 break;
665 case TypesSolver::BiCGStab:
666 HYPRE_ParCSRPilutCreate(MPI_COMM_SUB, &precond);
667 HYPRE_ParCSRBiCGSTABSetPrecond(solver,
668 HYPRE_ParCSRPilutSolve,
669 HYPRE_ParCSRPilutSetup,
670 precond);
671 break;
672 case TypesSolver::GMRES:
673 HYPRE_ParCSRPilutCreate(MPI_COMM_SUB,
674 &precond);
675 HYPRE_ParCSRGMRESSetPrecond(solver,
676 HYPRE_ParCSRPilutSolve,
677 HYPRE_ParCSRPilutSetup,
678 precond);
679 break;
680 default:
681 throw ArgumentException(func_name, "unknown solver for ILU preconditioner\n");
682 }
683 }
684
685 /******************************************************************************
686 *****************************************************************************/
687 void setSpaiStatPreconditioner(const TypesSolver::eSolverMethod solver_method,
688 const HYPRE_Solver& solver,
689 const AlephParams* solver_param,
690 HYPRE_Solver& precond)
691 {
692 HYPRE_ParCSRParaSailsCreate(MPI_COMM_SUB, &precond);
693 double alpha = solver_param->alpha();
694 int gamma = solver_param->gamma();
695 if (alpha < 0.0)
696 alpha = 0.1; // default value for the tolerance parameter
697 if (gamma == -1)
698 gamma = 1; // default value for the fill-in parameter
699 HYPRE_ParCSRParaSailsSetParams(precond, alpha, gamma);
700 switch (solver_method) {
701 case TypesSolver::PCG:
702 HYPRE_ParCSRPCGSetPrecond(solver, HYPRE_ParCSRParaSailsSolve, HYPRE_ParCSRParaSailsSetup, precond);
703 break;
704 case TypesSolver::BiCGStab:
705 throw ArgumentException("AlephMatrixHypre::setSpaiStatPreconditioner", "solveur 'BiCGStab' invalide pour preconditionnement 'SPAIstat'");
706 break;
707 case TypesSolver::GMRES:
708 // non-symmetric matrix
709 HYPRE_ParCSRParaSailsSetSym(precond, 0);
710 HYPRE_ParCSRGMRESSetPrecond(solver, HYPRE_ParaSailsSolve, HYPRE_ParaSailsSetup, precond);
711 break;
712 default:
713 throw ArgumentException("AlephMatrixHypre::setSpaiStatPreconditioner", "solveur inconnu pour preconditionnement 'SPAIstat'\n");
714 break;
715 }
716 }
717
718 /******************************************************************************
719 *****************************************************************************/
720 void setAMGPreconditioner(const TypesSolver::eSolverMethod solver_method,
721 const HYPRE_Solver& solver,
722 const AlephParams* solver_param,
723 HYPRE_Solver& precond)
724 {
725 // defaults for BoomerAMG from hypre example -- lc
726 // TODO : options and defaults values must be completed
727 double trunc_factor = 0.1; // set AMG interpolation truncation factor = val
728 int cycle_type = solver_param->getAmgCycle(); // set AMG cycles (1=V, 2=W, etc.)
729 int coarsen_type = solver_param->amgCoarseningMethod();
730 // Ruge coarsening (local) if <val> == 1
731 int relax_default = 3; // relaxation type <val> :
732 // 0=Weighted Jacobi
733 // 1=Gauss-Seidel (very slow!)
734 // 3=Hybrid Jacobi/Gauss-Seidel
735 int num_sweep = 1; // Use <val> sweeps on each level (here 1)
736 int hybrid = 1; // no switch in coarsening if -1
737 int measure_type = 1; // use global measures
738 double max_row_sum = 1.0; // set AMG maximum row sum threshold for dependency weakening
739
740 int max_levels = 50; // 25; // maximum number of AMG levels
741 const int gamma = solver_param->gamma();
742 if (gamma != -1)
743 max_levels = gamma; // use the dataset value
744
745 double strong_threshold = 0.1; // 0.25; // set AMG threshold Theta = val
746 const double alpha = solver_param->alpha();
747 if (alpha > 0.0)
748 strong_threshold = alpha; // use the dataset value
749 // news
750 Integer output_level = solver_param->getOutputLevel();
751
752 HYPRE_Int* num_grid_sweeps = _allocHypre<HYPRE_Int>(4);
753 HYPRE_Int* grid_relax_type = _allocHypre<HYPRE_Int>(4);
754 HYPRE_Int** grid_relax_points = _allocHypre<HYPRE_Int*>(4);
755 double* relax_weight = _allocHypre<double>(max_levels);
756
757 for (int i = 0; i < max_levels; i++)
758 relax_weight[i] = 1.0;
759
760 if (coarsen_type == 5) {
761 /* fine grid */
762 num_grid_sweeps[0] = 3;
763 grid_relax_type[0] = relax_default;
764 grid_relax_points[0] = _allocHypre<HYPRE_Int>(3);
765 grid_relax_points[0][0] = -2;
766 grid_relax_points[0][1] = -1;
767 grid_relax_points[0][2] = 1;
768
769 /* down cycle */
770 num_grid_sweeps[1] = 4;
771 grid_relax_type[1] = relax_default;
772 grid_relax_points[1] = _callocHypre<HYPRE_Int>(4);
773 grid_relax_points[1][0] = -1;
774 grid_relax_points[1][1] = 1;
775 grid_relax_points[1][2] = -2;
776 grid_relax_points[1][3] = -2;
777
778 /* up cycle */
779 num_grid_sweeps[2] = 4;
780 grid_relax_type[2] = relax_default;
781 grid_relax_points[2] = _allocHypre<HYPRE_Int>(4);
782 grid_relax_points[2][0] = -2;
783 grid_relax_points[2][1] = -2;
784 grid_relax_points[2][2] = 1;
785 grid_relax_points[2][3] = -1;
786 }
787 else {
788 /* fine grid */
789 num_grid_sweeps[0] = 2 * num_sweep;
790 grid_relax_type[0] = relax_default;
791 grid_relax_points[0] = _allocHypre<HYPRE_Int>(2 * num_sweep);
792 for (int i = 0; i < 2 * num_sweep; i += 2) {
793 grid_relax_points[0][i] = -1;
794 grid_relax_points[0][i + 1] = 1;
795 }
796
797 /* down cycle */
798 num_grid_sweeps[1] = 2 * num_sweep;
799 grid_relax_type[1] = relax_default;
800 grid_relax_points[1] = _allocHypre<HYPRE_Int>(2 * num_sweep);
801 for (int i = 0; i < 2 * num_sweep; i += 2) {
802 grid_relax_points[1][i] = -1;
803 grid_relax_points[1][i + 1] = 1;
804 }
805
806 /* up cycle */
807 num_grid_sweeps[2] = 2 * num_sweep;
808 grid_relax_type[2] = relax_default;
809 grid_relax_points[2] = _allocHypre<HYPRE_Int>(2 * num_sweep);
810 for (int i = 0; i < 2 * num_sweep; i += 2) {
811 grid_relax_points[2][i] = -1;
812 grid_relax_points[2][i + 1] = 1;
813 }
814 }
815
816 /* coarsest grid */
817 num_grid_sweeps[3] = 1;
818 grid_relax_type[3] = 9;
819 grid_relax_points[3] = _allocHypre<HYPRE_Int>(1);
820 grid_relax_points[3][0] = 0;
821
822 // end of default setting
823
824 HYPRE_BoomerAMGCreate(&precond);
825 HYPRE_BoomerAMGSetPrintLevel(precond, output_level);
826 HYPRE_BoomerAMGSetCoarsenType(precond, (hybrid * coarsen_type));
827 HYPRE_BoomerAMGSetMeasureType(precond, measure_type);
828 HYPRE_BoomerAMGSetStrongThreshold(precond, strong_threshold);
829 HYPRE_BoomerAMGSetTruncFactor(precond, trunc_factor);
830 HYPRE_BoomerAMGSetMaxIter(precond, 1);
831 HYPRE_BoomerAMGSetCycleType(precond, cycle_type);
832 HYPRE_BoomerAMGSetNumGridSweeps(precond, num_grid_sweeps);
833 HYPRE_BoomerAMGSetGridRelaxType(precond, grid_relax_type);
834 HYPRE_BoomerAMGSetRelaxWeight(precond, relax_weight);
835 HYPRE_BoomerAMGSetGridRelaxPoints(precond, grid_relax_points);
836 HYPRE_BoomerAMGSetTol(precond, 0.0);
837 HYPRE_BoomerAMGSetMaxLevels(precond, max_levels);
838 HYPRE_BoomerAMGSetMaxRowSum(precond, max_row_sum);
839
840 switch (solver_method) {
841 case TypesSolver::PCG:
842 HYPRE_ParCSRPCGSetPrecond(solver,
843 HYPRE_BoomerAMGSolve,
844 HYPRE_BoomerAMGSetup,
845 precond);
846 break;
847 case TypesSolver::BiCGStab:
848 HYPRE_ParCSRBiCGSTABSetPrecond(solver,
849 HYPRE_BoomerAMGSolve,
850 HYPRE_BoomerAMGSetup,
851 precond);
852 break;
853 case TypesSolver::GMRES:
854 HYPRE_ParCSRGMRESSetPrecond(solver,
855 HYPRE_BoomerAMGSolve,
856 HYPRE_BoomerAMGSetup,
857 precond);
858 break;
859 default:
860 throw ArgumentException("AlephMatrixHypre::setAMGPreconditioner", "solveur inconnu pour preconditionnement 'AMG'\n");
861 }
862 }
863
864 /******************************************************************************
865 *****************************************************************************/
866 bool solvePCG(const AlephParams* solver_param,
867 HYPRE_Solver& solver,
868 HYPRE_ParCSRMatrix& M,
869 HYPRE_ParVector& B,
870 HYPRE_ParVector& X,
871 HYPRE_Int& iteration,
872 double& residue)
873 {
874 const String func_name = "SolverMatrixHypre::solvePCG";
875 const bool xo = solver_param->xoUser();
876 bool error = false;
877
878 if (!xo)
879 HYPRE_ParVectorSetConstantValues(X, 0.0);
880 HYPRE_ParCSRPCGSetup(solver, M, B, X);
881 HYPRE_ParCSRPCGSolve(solver, M, B, X);
882 HYPRE_ParCSRPCGGetNumIterations(solver, &iteration);
883 HYPRE_ParCSRPCGGetFinalRelativeResidualNorm(solver, &residue);
884
885 HYPRE_Int converged = 0;
886 HYPRE_PCGGetConverged(solver, &converged);
887 error |= (!converged);
888
889 HYPRE_ParCSRPCGDestroy(solver);
890
891 return !error;
892 }
893
894 /******************************************************************************
895 *****************************************************************************/
896 bool solveBiCGStab(HYPRE_Solver& solver,
897 HYPRE_ParCSRMatrix& M,
898 HYPRE_ParVector& B,
899 HYPRE_ParVector& X,
900 HYPRE_Int& iteration,
901 double& residue)
902 {
903 const String func_name = "SolverMatrixHypre::solveBiCGStab";
904 bool error = false;
905 HYPRE_ParVectorSetRandomValues(X, 775);
906 HYPRE_ParCSRBiCGSTABSetup(solver, M, B, X);
907 HYPRE_ParCSRBiCGSTABSolve(solver, M, B, X);
908 HYPRE_ParCSRBiCGSTABGetNumIterations(solver, &iteration);
909 HYPRE_ParCSRBiCGSTABGetFinalRelativeResidualNorm(solver, &residue);
910
911 HYPRE_Int converged = 0;
912 hypre_BiCGSTABGetConverged(solver, &converged);
913 error |= (!converged);
914
915 HYPRE_ParCSRBiCGSTABDestroy(solver);
916
917 return !error;
918 }
919
920 /******************************************************************************
921 *****************************************************************************/
922 bool solveGMRES(HYPRE_Solver& solver,
923 HYPRE_ParCSRMatrix& M,
924 HYPRE_ParVector& B,
925 HYPRE_ParVector& X,
926 HYPRE_Int& iteration,
927 double& residue)
928 {
929 const String func_name = "SolverMatrixHypre::solveGMRES";
930 bool error = false;
931 HYPRE_ParCSRGMRESSetup(solver, M, B, X);
932 HYPRE_ParCSRGMRESSolve(solver, M, B, X);
933 HYPRE_ParCSRGMRESGetNumIterations(solver, &iteration);
934 HYPRE_ParCSRGMRESGetFinalRelativeResidualNorm(solver, &residue);
935
936 HYPRE_Int converged = 0;
937 HYPRE_GMRESGetConverged(solver, &converged);
938 error |= (!converged);
939
940 HYPRE_ParCSRGMRESDestroy(solver);
941 return !error;
942 }
943
944 private:
945
946 HYPRE_IJMatrix m_hypre_ijmatrix = nullptr;
947 HYPRE_ParCSRMatrix m_hypre_parmatrix = nullptr;
948};
949
950/*---------------------------------------------------------------------------*/
951/*---------------------------------------------------------------------------*/
952
953class HypreAlephFactoryImpl
954: public AbstractService
955, public IAlephFactoryImpl
956{
957 public:
958
959 HypreAlephFactoryImpl(const ServiceBuildInfo& sbi)
960 : AbstractService(sbi)
961 {}
962 ~HypreAlephFactoryImpl()
963 {
964 for (auto* v : m_IAlephVectors)
965 delete v;
966 for (auto* v : m_IAlephMatrixs)
967 delete v;
968 }
969
970 public:
971
972 void initialize() override
973 {
974 // NOTE: Starting from 2.29, we can use
975 // HYPRE_Initialize() and test if the initialization
976 // has already been done via HYPRE_Initialized().
977#if HYPRE_RELEASE_NUMBER >= 22900
978 if (!HYPRE_Initialized()) {
979 info() << "Initializing HYPRE";
980 HYPRE_Initialize();
981 }
982#elif HYPRE_RELEASE_NUMBER >= 22700
983 info() << "Initializing HYPRE";
984 HYPRE_Init();
985#endif
986
987#if HYPRE_RELEASE_NUMBER >= 22700
988 HYPRE_SetMemoryLocation(HYPRE_MEMORY_HOST);
989 HYPRE_SetExecutionPolicy(HYPRE_EXEC_HOST);
990#endif
991 }
992
993 IAlephTopology* createTopology(ITraceMng* tm,
994 AlephKernel* kernel,
995 Integer index,
996 Integer nb_row_size) override
997 {
998 ARCANE_UNUSED(tm);
999 ARCANE_UNUSED(kernel);
1000 ARCANE_UNUSED(index);
1001 ARCANE_UNUSED(nb_row_size);
1002 return NULL;
1003 }
1004
1005 IAlephVector* createVector(ITraceMng* tm,
1006 AlephKernel* kernel,
1007 Integer index) override
1008 {
1009 IAlephVector* new_vector = new AlephVectorHypre(tm, kernel, index);
1010 m_IAlephVectors.add(new_vector);
1011 return new_vector;
1012 }
1013
1014 IAlephMatrix* createMatrix(ITraceMng* tm,
1015 AlephKernel* kernel,
1016 Integer index) override
1017 {
1018 IAlephMatrix* new_matrix = new AlephMatrixHypre(tm, kernel, index);
1019 m_IAlephMatrixs.add(new_matrix);
1020 return new_matrix;
1021 }
1022
1023 private:
1024
1025 UniqueArray<IAlephVector*> m_IAlephVectors;
1026 UniqueArray<IAlephMatrix*> m_IAlephMatrixs;
1027};
1028
1029/*---------------------------------------------------------------------------*/
1030/*---------------------------------------------------------------------------*/
1031
1032ARCANE_REGISTER_APPLICATION_FACTORY(HypreAlephFactoryImpl, IAlephFactoryImpl, HypreAlephFactory);
1033
1034/*---------------------------------------------------------------------------*/
1035/*---------------------------------------------------------------------------*/
1036
1037} // namespace Arcane
1038
1039/*---------------------------------------------------------------------------*/
1040/*---------------------------------------------------------------------------*/
ARCANE_CHECK_POINTER
#define ARCANE_CHECK_POINTER(ptr)
Macro returning the pointer ptr if it is not null or throwing an exception if it is null.
Definition ArcaneGlobal.h:819
ARCANE_REGISTER_APPLICATION_FACTORY
#define ARCANE_REGISTER_APPLICATION_FACTORY(aclass, ainterface, aname)
Registers a factory service for the class aclass.
Definition ServiceFactory.h:636
Arcane::AbstractService::AbstractService
AbstractService(const ServiceBuildInfo &)
Constructor from a ServiceBuildInfo.
Definition AbstractService.cc:30
Arcane::AlephParams
Parameters of a linear system.
Definition AlephParams.h:32
Arcane::AlephVector
Vector of a linear system.
Definition AlephVector.h:33
Arcane::Array::data
const T * data() const
Access to the root of the array without any protection.
Definition arccore/src/common/arccore/common/Array.h:447
Arcane::IAlephFactoryImpl
Interface of an implementation factory for Aleph.
Definition AlephInterface.h:188
Arcane::ServiceBuildInfo
Structure containing the information to create a service.
Definition ServiceBuildInfo.h:199
Arcane::String::localstr
const char * localstr() const
Returns the conversion of the instance into UTF-8 encoding.
Definition String.cc:229
Arcane::TraceAccessor::debug
TraceMessageDbg debug(Trace::eDebugLevel=Trace::Medium) const
Flow for a debug message.
Definition arccore/src/trace/arccore/trace/TraceAccessor.h:122
Arcane::TraceAccessor::info
TraceMessage info() const
Flow for an information message.
Definition TraceAccessor.cc:101
Arcane::TraceAccessor::error
TraceMessage error() const
Flow for an error message.
Definition TraceAccessor.cc:201
Arcane::UniqueArray
1D data vector with value semantics (STL style).
Definition arccore/src/common/arccore/common/Array.h:890
Arcane::MessagePassing::ReduceMax
@ ReduceMax
Maximum of values.
Definition MessagePassingGlobal.h:113
Arcane
-- tab-width: 2; indent-tabs-mode: nil; coding: utf-8-with-signature --
Definition arcane/src/arcane/accelerator/AcceleratorGlobal.h:37
Arcane::Integer
Int32 Integer
Type representing an integer.
Definition ArccoreGlobal.h:293
Arcane::Real
double Real
Type representing a real number.
Definition ArccoreGlobal.h:275
Arcane::AlephInt
int AlephInt
Default type for indexing rows and columns of matrices and vectors.
Definition AlephGlobal.h:50