Arcane  4.2.3.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-2026 */
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 std::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 std::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 if (values[i] != 0.0)
341 rtn += HYPRE_IJMatrixSetValues(m_hypre_ijmatrix, 1, col, &rows[i], &cols[i], &values[i]);
342 }
343 hypreCheck("HYPRE_IJMatrixSetValues", rtn);
344 //HYPRE_IJMatrixSetValues(m_hypre_ijmatrix, nrows, ncols, rows, cols, values);
345 debug() << "[AlephMatrixHypre::AlephMatrixFill] done";
346 }
347
348 /******************************************************************************
349 * isAlreadySolved
350 *****************************************************************************/
351 bool isAlreadySolved(AlephVectorHypre* x,
352 AlephVectorHypre* b,
353 AlephVectorHypre* tmp,
354 Real* residual_norm,
355 AlephParams* params)
356 {
357 HYPRE_ClearAllErrors();
358 const bool convergence_analyse = params->convergenceAnalyse();
359
360 // test the right-hand side of the linear system
361 const Real res0 = b->norm_max();
362
363 if (convergence_analyse)
364 info() << "convergence analysis: max norm of the right-hand side res0: " << res0;
365
366 const Real considered_as_null = params->minRHSNorm();
367 if (res0 < considered_as_null) {
368 HYPRE_ParVectorSetConstantValues(x->m_hypre_parvector, 0.0);
369 residual_norm[0] = res0;
370 if (convergence_analyse)
371 info() << "convergence analysis: the right-hand side of the linear system is less than: " << considered_as_null;
372 return true;
373 }
374
375 if (params->xoUser()) {
376 // we test if b is already a solution to the system within epsilon tolerance
377 //matrix->vectorProduct(b, tmp_vector); tmp_vector->sub(x);
378 //M->vector_multiply(*tmp,*x); // tmp=A*x
379 //tmp->substract(*tmp,*b); // tmp=A*x-b
380
381 // X= alpha* M.B + beta * X (read from HYPRE sources)
382 HYPRE_ParCSRMatrixMatvec(1.0, m_hypre_parmatrix, x->m_hypre_parvector, 0., tmp->m_hypre_parvector);
383 HYPRE_ParVectorAxpy(-1.0, b->m_hypre_parvector, tmp->m_hypre_parvector);
384
385 const Real residu = tmp->norm_max();
386 //info() << "[IAlephHypre::isAlreadySolved] residual="<<residu;
387
388 if (residu < considered_as_null) {
389 if (convergence_analyse) {
390 info() << "convergence analysis: |Ax0-b| is less than " << considered_as_null;
391 info() << "convergence analysis: x0 is already a solution to the system.";
392 }
393 residual_norm[0] = residu;
394 return true;
395 }
396
397 const Real relative_error = residu / res0;
398 if (convergence_analyse)
399 info() << "convergence analysis: initial residual : " << residu
400 << " --- initial relative residual (residu/res0) : " << residu / res0;
401
402 if (relative_error < (params->epsilon())) {
403 if (convergence_analyse)
404 info() << "convergence analysis: X is already a solution to the system";
405 residual_norm[0] = residu;
406 return true;
407 }
408 }
409 return false;
410 }
411
412 /******************************************************************************
413 *****************************************************************************/
414 int AlephMatrixSolve(AlephVector* x,
415 AlephVector* b,
416 AlephVector* t,
417 Integer& nb_iteration,
418 Real* residual_norm,
419 AlephParams* solver_param)
420 {
421 solver_param->setAmgCoarseningMethod(TypesSolver::AMG_COARSENING_AUTO);
422 const String func_name("SolverMatrixHypre::solve");
423 void* object;
424 int ierr = 0;
425
426 auto* ximpl = ARCANE_CHECK_POINTER(dynamic_cast<AlephVectorHypre*>(x->implementation()));
427 auto* bimpl = ARCANE_CHECK_POINTER(dynamic_cast<AlephVectorHypre*>(b->implementation()));
428
429 HYPRE_IJVector solution = ximpl->m_hypre_ijvector;
430 HYPRE_IJVector RHS = bimpl->m_hypre_ijvector;
431 //HYPRE_IJVector tmp = (dynamic_cast<AlephVectorHypre*> (t->implementation()))->m_hypre_ijvector;
432
433 HYPRE_IJMatrixGetObject(m_hypre_ijmatrix, &object);
434 HYPRE_ParCSRMatrix M = (HYPRE_ParCSRMatrix)object;
435 HYPRE_IJVectorGetObject(solution, &object);
436 HYPRE_ParVector X = (HYPRE_ParVector)object;
437 HYPRE_IJVectorGetObject(RHS, &object);
438 HYPRE_ParVector B = (HYPRE_ParVector)object;
439 //HYPRE_IJVectorGetObject(tmp,&object);
440 //HYPRE_ParVector T = (HYPRE_ParVector)object;
441
442 auto* ximpl2 = ARCANE_CHECK_POINTER(dynamic_cast<AlephVectorHypre*>(x->implementation()));
443 auto* bimpl2 = ARCANE_CHECK_POINTER(dynamic_cast<AlephVectorHypre*>(b->implementation()));
444 auto* timpl2 = ARCANE_CHECK_POINTER(dynamic_cast<AlephVectorHypre*>(t->implementation()));
445 if (isAlreadySolved(ximpl2, bimpl2, timpl2, residual_norm, solver_param)) {
446 ItacRegion(isAlreadySolved, AlephMatrixHypre);
447 debug() << "[AlephMatrixHypre::AlephMatrixSolve] isAlreadySolved !";
448 nb_iteration = 0;
449 return 0;
450 }
451
452 TypesSolver::ePreconditionerMethod preconditioner_method = solver_param->precond();
453 // TypesSolver::ePreconditionerMethod preconditioner_method = TypesSolver::NONE;
454 TypesSolver::eSolverMethod solver_method = solver_param->method();
455
456 // declaration and initialization of the solver
457 HYPRE_Solver solver = 0;
458
459 switch (solver_method) {
460 case TypesSolver::PCG:
461 initSolverPCG(solver_param, solver);
462 break;
463 case TypesSolver::BiCGStab:
464 initSolverBiCGStab(solver_param, solver);
465 break;
466 case TypesSolver::GMRES:
467 initSolverGMRES(solver_param, solver);
468 break;
469 default:
470 throw ArgumentException(func_name, "unknown solver");
471 }
472
473 // declaration and initialization of the preconditioner
474 HYPRE_Solver precond = 0;
475
476 switch (preconditioner_method) {
477 case TypesSolver::NONE:
478 break;
479 case TypesSolver::DIAGONAL:
480 setDiagonalPreconditioner(solver_method, solver, precond);
481 break;
482 case TypesSolver::ILU:
483 setILUPreconditioner(solver_method, solver, precond);
484 break;
485 case TypesSolver::SPAIstat:
486 setSpaiStatPreconditioner(solver_method, solver, solver_param, precond);
487 break;
488 case TypesSolver::AMG:
489 setAMGPreconditioner(solver_method, solver, solver_param, precond);
490 break;
491 case TypesSolver::AINV:
492 throw ArgumentException(func_name, "AINV preconditioning unavailable");
493 case TypesSolver::SPAIdyn:
494 throw ArgumentException(func_name, "SPAIdyn preconditioning unavailable");
495 case TypesSolver::ILUp:
496 throw ArgumentException(func_name, "ILUp preconditioning unavailable");
497 case TypesSolver::IC:
498 throw ArgumentException(func_name, "IC preconditioning unavailable");
499 case TypesSolver::POLY:
500 throw ArgumentException(func_name, "POLY preconditioning unavailable");
501 default:
502 throw ArgumentException(func_name, "unknown preconditioner");
503 }
504
505 // solving the algebraic system
506 HYPRE_Int iteration = 0;
507 double residue = 0.0;
508
509 switch (solver_method) {
510 case TypesSolver::PCG:
511 ierr = solvePCG(solver_param, solver, M, B, X, iteration, residue);
512 break;
513 case TypesSolver::BiCGStab:
514 ierr = solveBiCGStab(solver, M, B, X, iteration, residue);
515 break;
516 case TypesSolver::GMRES:
517 ierr = solveGMRES(solver, M, B, X, iteration, residue);
518 break;
519 default:
520 ierr = -3;
521 return ierr;
522 }
523 nb_iteration = static_cast<Integer>(iteration);
524 residual_norm[0] = static_cast<Real>(residue);
525
526 /* for(int i=0;i<8;++i){
527 int idx[1];
528 double valx[1]={-1.};
529 //double valb[1]={-1.};
530 idx[0]=i;
531 HYPRE_IJVectorGetValues(solution, 1, idx, valx);
532 debug()<<"[AlephMatrixHypre::AlephMatrixSolve] X["<<i<<"]="<<valx[0];
533 //(static_cast<AlephVectorHypre*> (x->implementation()))->AlephVectorGet(valx,idx,1);
534 //debug()<<"[AlephMatrixHypre::AlephMatrixSolve] x["<<i<<"]="<<valx[0];
535 //HYPRE_IJVectorGetValues(RHS, 1, idx, valb);
536 //debug()<<"[AlephMatrixHypre::AlephMatrixSolve] B["<<i<<"]="<<valb[0];
537 }
538*/
539 switch (preconditioner_method) {
540 case TypesSolver::NONE:
541 break;
542 case TypesSolver::DIAGONAL:
543 break;
544 case TypesSolver::ILU:
545 HYPRE_ParCSRPilutDestroy(precond);
546 break;
547 case TypesSolver::SPAIstat:
548 HYPRE_ParCSRParaSailsDestroy(precond);
549 break;
550 case TypesSolver::AMG:
551 HYPRE_BoomerAMGDestroy(precond);
552 break;
553 default:
554 throw ArgumentException(func_name, "unknown preconditioner");
555 }
556
557 if (iteration == solver_param->maxIter() && solver_param->stopErrorStrategy()) {
558 info() << "\n============================================================";
559 info() << "\nThis error is returned after " << iteration << "\n";
560 info() << "\nMaximum number of solver iterations reached.";
561 info() << "\nIt is possible to ask the code not to consider this error.";
562 info() << "\nSee the dataset documentation regarding the solver service.";
563 info() << "\n======================================================";
564 throw Exception("AlephMatrixHypre::Solve", "Maximum number of solver iterations reached");
565 }
566 return ierr;
567 }
568
569 /******************************************************************************
570 *****************************************************************************/
571 void writeToFile(const String filename)
572 {
573 HYPRE_IJMatrixPrint(m_hypre_ijmatrix, filename.localstr());
574 }
575
576 /******************************************************************************
577 *****************************************************************************/
578 void initSolverPCG(const AlephParams* solver_param, HYPRE_Solver& solver)
579 {
580 const String func_name = "SolverMatrixHypre::initSolverPCG";
581 double epsilon = solver_param->epsilon();
582 int max_it = solver_param->maxIter();
583 int output_level = solver_param->getOutputLevel();
584
585 HYPRE_ParCSRPCGCreate(MPI_COMM_SUB, &solver);
586 HYPRE_ParCSRPCGSetMaxIter(solver, max_it);
587 HYPRE_ParCSRPCGSetTol(solver, epsilon);
588 HYPRE_ParCSRPCGSetTwoNorm(solver, 1);
589 HYPRE_ParCSRPCGSetPrintLevel(solver, output_level);
590 }
591
592 /******************************************************************************
593 *****************************************************************************/
594 void initSolverBiCGStab(const AlephParams* solver_param, HYPRE_Solver& solver)
595 {
596 const String func_name = "SolverMatrixHypre::initSolverBiCGStab";
597 double epsilon = solver_param->epsilon();
598 int max_it = solver_param->maxIter();
599 int output_level = solver_param->getOutputLevel();
600
601 HYPRE_ParCSRBiCGSTABCreate(MPI_COMM_SUB, &solver);
602 HYPRE_ParCSRBiCGSTABSetMaxIter(solver, max_it);
603 HYPRE_ParCSRBiCGSTABSetTol(solver, epsilon);
604 HYPRE_ParCSRBiCGSTABSetPrintLevel(solver, output_level);
605 }
606
607 /******************************************************************************
608 *****************************************************************************/
609 void initSolverGMRES(const AlephParams* solver_param, HYPRE_Solver& solver)
610 {
611 const String func_name = "SolverMatrixHypre::initSolverGMRES";
612 double epsilon = solver_param->epsilon();
613 int max_it = solver_param->maxIter();
614 int output_level = solver_param->getOutputLevel();
615
616 HYPRE_ParCSRGMRESCreate(MPI_COMM_SUB, &solver);
617 const int krylov_dim = 20; // dimension Krylov space for GMRES
618 HYPRE_ParCSRGMRESSetKDim(solver, krylov_dim);
619 HYPRE_ParCSRGMRESSetMaxIter(solver, max_it);
620 HYPRE_ParCSRGMRESSetTol(solver, epsilon);
621 HYPRE_ParCSRGMRESSetPrintLevel(solver, output_level);
622 }
623
624 /******************************************************************************
625 *****************************************************************************/
626 void setDiagonalPreconditioner(const TypesSolver::eSolverMethod solver_method,
627 const HYPRE_Solver& solver,
628 HYPRE_Solver& precond)
629 {
630 const String func_name = "SolverMatrixHypre::setDiagonalPreconditioner";
631 switch (solver_method) {
632 case TypesSolver::PCG:
633 HYPRE_ParCSRPCGSetPrecond(solver,
634 HYPRE_ParCSRDiagScale,
635 HYPRE_ParCSRDiagScaleSetup,
636 precond);
637 break;
638 case TypesSolver::BiCGStab:
639 HYPRE_ParCSRBiCGSTABSetPrecond(solver,
640 HYPRE_ParCSRDiagScale,
641 HYPRE_ParCSRDiagScaleSetup,
642 precond);
643 break;
644 case TypesSolver::GMRES:
645 HYPRE_ParCSRGMRESSetPrecond(solver,
646 HYPRE_ParCSRDiagScale,
647 HYPRE_ParCSRDiagScaleSetup,
648 precond);
649 break;
650 default:
651 throw ArgumentException(func_name, "unknown solver for 'Diagonal' preconditioner");
652 }
653 }
654
655 /******************************************************************************
656 *****************************************************************************/
657 void setILUPreconditioner(const TypesSolver::eSolverMethod solver_method,
658 const HYPRE_Solver& solver,
659 HYPRE_Solver& precond)
660 {
661 const String func_name = "SolverMatrixHypre::setILUPreconditioner";
662 switch (solver_method) {
663 case TypesSolver::PCG:
664 throw ArgumentException(func_name, "PCG solver unavailable with 'ILU' preconditioner");
665 break;
666 case TypesSolver::BiCGStab:
667 HYPRE_ParCSRPilutCreate(MPI_COMM_SUB, &precond);
668 HYPRE_ParCSRBiCGSTABSetPrecond(solver,
669 HYPRE_ParCSRPilutSolve,
670 HYPRE_ParCSRPilutSetup,
671 precond);
672 break;
673 case TypesSolver::GMRES:
674 HYPRE_ParCSRPilutCreate(MPI_COMM_SUB,
675 &precond);
676 HYPRE_ParCSRGMRESSetPrecond(solver,
677 HYPRE_ParCSRPilutSolve,
678 HYPRE_ParCSRPilutSetup,
679 precond);
680 break;
681 default:
682 throw ArgumentException(func_name, "unknown solver for ILU preconditioner\n");
683 }
684 }
685
686 /******************************************************************************
687 *****************************************************************************/
688 void setSpaiStatPreconditioner(const TypesSolver::eSolverMethod solver_method,
689 const HYPRE_Solver& solver,
690 const AlephParams* solver_param,
691 HYPRE_Solver& precond)
692 {
693 HYPRE_ParCSRParaSailsCreate(MPI_COMM_SUB, &precond);
694 double alpha = solver_param->alpha();
695 int gamma = solver_param->gamma();
696 if (alpha < 0.0)
697 alpha = 0.1; // default value for the tolerance parameter
698 if (gamma == -1)
699 gamma = 1; // default value for the fill-in parameter
700 HYPRE_ParCSRParaSailsSetParams(precond, alpha, gamma);
701 switch (solver_method) {
702 case TypesSolver::PCG:
703 HYPRE_ParCSRPCGSetPrecond(solver, HYPRE_ParCSRParaSailsSolve, HYPRE_ParCSRParaSailsSetup, precond);
704 break;
705 case TypesSolver::BiCGStab:
706 throw ArgumentException("AlephMatrixHypre::setSpaiStatPreconditioner", "solveur 'BiCGStab' invalide pour preconditionnement 'SPAIstat'");
707 break;
708 case TypesSolver::GMRES:
709 // non-symmetric matrix
710 HYPRE_ParCSRParaSailsSetSym(precond, 0);
711 HYPRE_ParCSRGMRESSetPrecond(solver, HYPRE_ParaSailsSolve, HYPRE_ParaSailsSetup, precond);
712 break;
713 default:
714 throw ArgumentException("AlephMatrixHypre::setSpaiStatPreconditioner", "solveur inconnu pour preconditionnement 'SPAIstat'\n");
715 break;
716 }
717 }
718
719 /******************************************************************************
720 *****************************************************************************/
721 void setAMGPreconditioner(const TypesSolver::eSolverMethod solver_method,
722 const HYPRE_Solver& solver,
723 const AlephParams* solver_param,
724 HYPRE_Solver& precond)
725 {
726 // defaults for BoomerAMG from hypre example -- lc
727 // TODO : options and defaults values must be completed
728 double trunc_factor = 0.1; // set AMG interpolation truncation factor = val
729 int cycle_type = solver_param->getAmgCycle(); // set AMG cycles (1=V, 2=W, etc.)
730 int coarsen_type = solver_param->amgCoarseningMethod();
731 // Ruge coarsening (local) if <val> == 1
732 int relax_default = 3; // relaxation type <val> :
733 // 0=Weighted Jacobi
734 // 1=Gauss-Seidel (very slow!)
735 // 3=Hybrid Jacobi/Gauss-Seidel
736 int num_sweep = 1; // Use <val> sweeps on each level (here 1)
737 int hybrid = 1; // no switch in coarsening if -1
738 int measure_type = 1; // use global measures
739 double max_row_sum = 1.0; // set AMG maximum row sum threshold for dependency weakening
740
741 int max_levels = 50; // 25; // maximum number of AMG levels
742 const int gamma = solver_param->gamma();
743 if (gamma != -1)
744 max_levels = gamma; // use the dataset value
745
746 double strong_threshold = 0.1; // 0.25; // set AMG threshold Theta = val
747 const double alpha = solver_param->alpha();
748 if (alpha > 0.0)
749 strong_threshold = alpha; // use the dataset value
750 // news
751 Integer output_level = solver_param->getOutputLevel();
752
753 HYPRE_Int* num_grid_sweeps = _allocHypre<HYPRE_Int>(4);
754 HYPRE_Int* grid_relax_type = _allocHypre<HYPRE_Int>(4);
755 HYPRE_Int** grid_relax_points = _allocHypre<HYPRE_Int*>(4);
756 double* relax_weight = _allocHypre<double>(max_levels);
757
758 for (int i = 0; i < max_levels; i++)
759 relax_weight[i] = 1.0;
760
761 if (coarsen_type == 5) {
762 /* fine grid */
763 num_grid_sweeps[0] = 3;
764 grid_relax_type[0] = relax_default;
765 grid_relax_points[0] = _allocHypre<HYPRE_Int>(3);
766 grid_relax_points[0][0] = -2;
767 grid_relax_points[0][1] = -1;
768 grid_relax_points[0][2] = 1;
769
770 /* down cycle */
771 num_grid_sweeps[1] = 4;
772 grid_relax_type[1] = relax_default;
773 grid_relax_points[1] = _callocHypre<HYPRE_Int>(4);
774 grid_relax_points[1][0] = -1;
775 grid_relax_points[1][1] = 1;
776 grid_relax_points[1][2] = -2;
777 grid_relax_points[1][3] = -2;
778
779 /* up cycle */
780 num_grid_sweeps[2] = 4;
781 grid_relax_type[2] = relax_default;
782 grid_relax_points[2] = _allocHypre<HYPRE_Int>(4);
783 grid_relax_points[2][0] = -2;
784 grid_relax_points[2][1] = -2;
785 grid_relax_points[2][2] = 1;
786 grid_relax_points[2][3] = -1;
787 }
788 else {
789 /* fine grid */
790 num_grid_sweeps[0] = 2 * num_sweep;
791 grid_relax_type[0] = relax_default;
792 grid_relax_points[0] = _allocHypre<HYPRE_Int>(2 * num_sweep);
793 for (int i = 0; i < 2 * num_sweep; i += 2) {
794 grid_relax_points[0][i] = -1;
795 grid_relax_points[0][i + 1] = 1;
796 }
797
798 /* down cycle */
799 num_grid_sweeps[1] = 2 * num_sweep;
800 grid_relax_type[1] = relax_default;
801 grid_relax_points[1] = _allocHypre<HYPRE_Int>(2 * num_sweep);
802 for (int i = 0; i < 2 * num_sweep; i += 2) {
803 grid_relax_points[1][i] = -1;
804 grid_relax_points[1][i + 1] = 1;
805 }
806
807 /* up cycle */
808 num_grid_sweeps[2] = 2 * num_sweep;
809 grid_relax_type[2] = relax_default;
810 grid_relax_points[2] = _allocHypre<HYPRE_Int>(2 * num_sweep);
811 for (int i = 0; i < 2 * num_sweep; i += 2) {
812 grid_relax_points[2][i] = -1;
813 grid_relax_points[2][i + 1] = 1;
814 }
815 }
816
817 /* coarsest grid */
818 num_grid_sweeps[3] = 1;
819 grid_relax_type[3] = 9;
820 grid_relax_points[3] = _allocHypre<HYPRE_Int>(1);
821 grid_relax_points[3][0] = 0;
822
823 // end of default setting
824
825 HYPRE_BoomerAMGCreate(&precond);
826 HYPRE_BoomerAMGSetPrintLevel(precond, output_level);
827 HYPRE_BoomerAMGSetCoarsenType(precond, (hybrid * coarsen_type));
828 HYPRE_BoomerAMGSetMeasureType(precond, measure_type);
829 HYPRE_BoomerAMGSetStrongThreshold(precond, strong_threshold);
830 HYPRE_BoomerAMGSetTruncFactor(precond, trunc_factor);
831 HYPRE_BoomerAMGSetMaxIter(precond, 1);
832 HYPRE_BoomerAMGSetCycleType(precond, cycle_type);
833 HYPRE_BoomerAMGSetNumGridSweeps(precond, num_grid_sweeps);
834 HYPRE_BoomerAMGSetGridRelaxType(precond, grid_relax_type);
835 HYPRE_BoomerAMGSetRelaxWeight(precond, relax_weight);
836 HYPRE_BoomerAMGSetGridRelaxPoints(precond, grid_relax_points);
837 HYPRE_BoomerAMGSetTol(precond, 0.0);
838 HYPRE_BoomerAMGSetMaxLevels(precond, max_levels);
839 HYPRE_BoomerAMGSetMaxRowSum(precond, max_row_sum);
840
841 switch (solver_method) {
842 case TypesSolver::PCG:
843 HYPRE_ParCSRPCGSetPrecond(solver,
844 HYPRE_BoomerAMGSolve,
845 HYPRE_BoomerAMGSetup,
846 precond);
847 break;
848 case TypesSolver::BiCGStab:
849 HYPRE_ParCSRBiCGSTABSetPrecond(solver,
850 HYPRE_BoomerAMGSolve,
851 HYPRE_BoomerAMGSetup,
852 precond);
853 break;
854 case TypesSolver::GMRES:
855 HYPRE_ParCSRGMRESSetPrecond(solver,
856 HYPRE_BoomerAMGSolve,
857 HYPRE_BoomerAMGSetup,
858 precond);
859 break;
860 default:
861 throw ArgumentException("AlephMatrixHypre::setAMGPreconditioner", "solveur inconnu pour preconditionnement 'AMG'\n");
862 }
863 }
864
865 /******************************************************************************
866 *****************************************************************************/
867 bool solvePCG(const AlephParams* solver_param,
868 HYPRE_Solver& solver,
869 HYPRE_ParCSRMatrix& M,
870 HYPRE_ParVector& B,
871 HYPRE_ParVector& X,
872 HYPRE_Int& iteration,
873 double& residue)
874 {
875 const String func_name = "SolverMatrixHypre::solvePCG";
876 const bool xo = solver_param->xoUser();
877 bool error = false;
878
879 if (!xo)
880 HYPRE_ParVectorSetConstantValues(X, 0.0);
881 HYPRE_ParCSRPCGSetup(solver, M, B, X);
882 HYPRE_ParCSRPCGSolve(solver, M, B, X);
883 HYPRE_ParCSRPCGGetNumIterations(solver, &iteration);
884 HYPRE_ParCSRPCGGetFinalRelativeResidualNorm(solver, &residue);
885
886 HYPRE_Int converged = 0;
887 HYPRE_PCGGetConverged(solver, &converged);
888 error |= (!converged);
889
890 HYPRE_ParCSRPCGDestroy(solver);
891
892 return !error;
893 }
894
895 /******************************************************************************
896 *****************************************************************************/
897 bool solveBiCGStab(HYPRE_Solver& solver,
898 HYPRE_ParCSRMatrix& M,
899 HYPRE_ParVector& B,
900 HYPRE_ParVector& X,
901 HYPRE_Int& iteration,
902 double& residue)
903 {
904 const String func_name = "SolverMatrixHypre::solveBiCGStab";
905 bool error = false;
906 HYPRE_ParVectorSetRandomValues(X, 775);
907 HYPRE_ParCSRBiCGSTABSetup(solver, M, B, X);
908 HYPRE_ParCSRBiCGSTABSolve(solver, M, B, X);
909 HYPRE_ParCSRBiCGSTABGetNumIterations(solver, &iteration);
910 HYPRE_ParCSRBiCGSTABGetFinalRelativeResidualNorm(solver, &residue);
911
912 HYPRE_Int converged = 0;
913 hypre_BiCGSTABGetConverged(solver, &converged);
914 error |= (!converged);
915
916 HYPRE_ParCSRBiCGSTABDestroy(solver);
917
918 return !error;
919 }
920
921 /******************************************************************************
922 *****************************************************************************/
923 bool solveGMRES(HYPRE_Solver& solver,
924 HYPRE_ParCSRMatrix& M,
925 HYPRE_ParVector& B,
926 HYPRE_ParVector& X,
927 HYPRE_Int& iteration,
928 double& residue)
929 {
930 const String func_name = "SolverMatrixHypre::solveGMRES";
931 bool error = false;
932 HYPRE_ParCSRGMRESSetup(solver, M, B, X);
933 HYPRE_ParCSRGMRESSolve(solver, M, B, X);
934 HYPRE_ParCSRGMRESGetNumIterations(solver, &iteration);
935 HYPRE_ParCSRGMRESGetFinalRelativeResidualNorm(solver, &residue);
936
937 HYPRE_Int converged = 0;
938 HYPRE_GMRESGetConverged(solver, &converged);
939 error |= (!converged);
940
941 HYPRE_ParCSRGMRESDestroy(solver);
942 return !error;
943 }
944
945 private:
946
947 HYPRE_IJMatrix m_hypre_ijmatrix = nullptr;
948 HYPRE_ParCSRMatrix m_hypre_parmatrix = nullptr;
949};
950
951/*---------------------------------------------------------------------------*/
952/*---------------------------------------------------------------------------*/
953
954class HypreAlephFactoryImpl
955: public AbstractService
956, public IAlephFactoryImpl
957{
958 public:
959
960 HypreAlephFactoryImpl(const ServiceBuildInfo& sbi)
961 : AbstractService(sbi)
962 {}
963 ~HypreAlephFactoryImpl()
964 {
965 for (auto* v : m_IAlephVectors)
966 delete v;
967 for (auto* v : m_IAlephMatrixs)
968 delete v;
969 }
970
971 public:
972
973 void initialize() override
974 {
975 // NOTE: Starting from 2.29, we can use
976 // HYPRE_Initialize() and test if the initialization
977 // has already been done via HYPRE_Initialized().
978#if HYPRE_RELEASE_NUMBER >= 22900
979 if (!HYPRE_Initialized()) {
980 info() << "Initializing HYPRE";
981 HYPRE_Initialize();
982 }
983#elif HYPRE_RELEASE_NUMBER >= 22700
984 info() << "Initializing HYPRE";
985 HYPRE_Init();
986#endif
987
988#if HYPRE_RELEASE_NUMBER >= 22700
989 HYPRE_SetMemoryLocation(HYPRE_MEMORY_HOST);
990 HYPRE_SetExecutionPolicy(HYPRE_EXEC_HOST);
991#endif
992 }
993
994 IAlephTopology* createTopology(ITraceMng* tm,
995 AlephKernel* kernel,
996 Integer index,
997 Integer nb_row_size) override
998 {
999 ARCANE_UNUSED(tm);
1000 ARCANE_UNUSED(kernel);
1001 ARCANE_UNUSED(index);
1002 ARCANE_UNUSED(nb_row_size);
1003 return NULL;
1004 }
1005
1006 IAlephVector* createVector(ITraceMng* tm,
1007 AlephKernel* kernel,
1008 Integer index) override
1009 {
1010 IAlephVector* new_vector = new AlephVectorHypre(tm, kernel, index);
1011 m_IAlephVectors.add(new_vector);
1012 return new_vector;
1013 }
1014
1015 IAlephMatrix* createMatrix(ITraceMng* tm,
1016 AlephKernel* kernel,
1017 Integer index) override
1018 {
1019 IAlephMatrix* new_matrix = new AlephMatrixHypre(tm, kernel, index);
1020 m_IAlephMatrixs.add(new_matrix);
1021 return new_matrix;
1022 }
1023
1024 private:
1025
1026 UniqueArray<IAlephVector*> m_IAlephVectors;
1027 UniqueArray<IAlephMatrix*> m_IAlephMatrixs;
1028};
1029
1030/*---------------------------------------------------------------------------*/
1031/*---------------------------------------------------------------------------*/
1032
1033ARCANE_REGISTER_APPLICATION_FACTORY(HypreAlephFactoryImpl, IAlephFactoryImpl, HypreAlephFactory);
1034
1035/*---------------------------------------------------------------------------*/
1036/*---------------------------------------------------------------------------*/
1037
1038} // namespace Arcane
1039
1040/*---------------------------------------------------------------------------*/
1041/*---------------------------------------------------------------------------*/
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:815
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:511
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:974
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