GEOS
HypreUtils.hpp
Go to the documentation of this file.
1 /*
2  * ------------------------------------------------------------------------------------------------------------
3  * SPDX-License-Identifier: LGPL-2.1-only
4  *
5  * Copyright (c) 2016-2024 Lawrence Livermore National Security LLC
6  * Copyright (c) 2018-2024 TotalEnergies
7  * Copyright (c) 2018-2024 The Board of Trustees of the Leland Stanford Junior University
8  * Copyright (c) 2023-2024 Chevron
9  * Copyright (c) 2019- GEOS/GEOSX Contributors
10  * All rights reserved
11  *
12  * See top level LICENSE, COPYRIGHT, CONTRIBUTORS, NOTICE, and ACKNOWLEDGEMENTS files for details.
13  * ------------------------------------------------------------------------------------------------------------
14  */
15 
20 #ifndef GEOS_LINEARALGEBRA_INTERFACES_HYPREUTILS_HPP_
21 #define GEOS_LINEARALGEBRA_INTERFACES_HYPREUTILS_HPP_
22 
23 #include "common/DataTypes.hpp"
24 #include "common/GEOS_RAJA_Interface.hpp"
25 
26 #include "codingUtilities/Utilities.hpp"
28 
29 #include <HYPRE_krylov.h>
30 #include <HYPRE_parcsr_ls.h>
31 
32 #if GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_CUDA || GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_HIP
34 #define GEOS_HYPRE_DEVICE GEOS_DEVICE
35 #else
37 #define GEOS_HYPRE_DEVICE
38 #endif
39 
40 namespace geos
41 {
42 
50 {
52  using SetupFunc = HYPRE_Int (*)( HYPRE_Solver, HYPRE_ParCSRMatrix, HYPRE_ParVector, HYPRE_ParVector );
53 
55  using SolveFunc = HYPRE_Int (*)( HYPRE_Solver, HYPRE_ParCSRMatrix, HYPRE_ParVector, HYPRE_ParVector );
56 
58  using DestroyFunc = HYPRE_Int (*)( HYPRE_Solver );
59 
60  HYPRE_Solver ptr{};
64 };
65 
69 namespace hypre
70 {
71 
76 constexpr HYPRE_MemoryLocation getMemoryLocation( LvArray::MemorySpace const space )
77 {
78  switch( space )
79  {
80  case hostMemorySpace: return HYPRE_MEMORY_HOST;
81 #if GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_CUDA
82  case LvArray::MemorySpace::cuda: return HYPRE_MEMORY_DEVICE;
83 #endif
84 #if GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_HIP
85  case LvArray::MemorySpace::hip: return HYPRE_MEMORY_DEVICE;
86 #endif
87  default: return HYPRE_MEMORY_HOST;
88  }
89 }
90 
95 constexpr LvArray::MemorySpace getLvArrayMemorySpace( HYPRE_MemoryLocation const location )
96 {
97  switch( location )
98  {
99  case HYPRE_MEMORY_HOST: return hostMemorySpace;
100 #if GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_CUDA
101  case HYPRE_MEMORY_DEVICE: return parallelDeviceMemorySpace;
102 #endif
103 #if GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_HIP
104  case HYPRE_MEMORY_DEVICE: return parallelDeviceMemorySpace;
105 #endif
106  default: return hostMemorySpace;
107  }
108 }
109 
110 #if GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_CUDA || GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_HIP
111 
113 using execPolicy = parallelDevicePolicy<>;
115 constexpr LvArray::MemorySpace memorySpace = parallelDeviceMemorySpace;
117 constexpr HYPRE_MemoryLocation memoryLocation = HYPRE_MEMORY_DEVICE;
118 
119 #else
120 
122 using execPolicy = parallelHostPolicy;
124 constexpr LvArray::MemorySpace memorySpace = hostMemorySpace;
126 constexpr HYPRE_MemoryLocation memoryLocation = HYPRE_MEMORY_HOST;
127 
128 #endif
129 
130 // Check matching requirements on index/value types between GEOS and Hypre
131 
132 // WARNING. We don't have consistent types between HYPRE_Int and localIndex.
133 // Decision needs to be made either to use bigint option, or change
134 // localIndex to int. We are getting away with this because we do not
135 // pass ( localIndex * ) to hypre except when it is on the GPU, in
136 // which case we are using int for localIndex.
137 #if GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_CUDA || GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_HIP
138 static_assert( sizeof( HYPRE_Int ) == sizeof( geos::localIndex ),
139  "HYPRE_Int and geos::localIndex must have the same size" );
140 static_assert( std::is_signed< HYPRE_Int >::value == std::is_signed< geos::localIndex >::value,
141  "HYPRE_Int and geos::localIndex must both be signed or unsigned" );
142 #endif
143 
150 inline void checkDeviceErrors( char const * msg, char const * file, int const line )
151 {
152 #if GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_CUDA
153  cudaError_t const err = cudaGetLastError();
154  GEOS_ERROR_IF( err != cudaSuccess, GEOS_FMT( "Previous CUDA errors found: {} ({} at {}:{})", msg, cudaGetErrorString( err ), file, line ) );
155 #elif GEOS_USE_HYPRE_DEVICE == GEOS_USE_HYPRE_HIP
156  hipError_t const err = hipGetLastError();
157  GEOS_UNUSED_VAR( msg, file, line ); // on crusher geos_error_if ultimately resolves to an assert, which drops the content on release
158  // builds
159  GEOS_ERROR_IF( err != hipSuccess, GEOS_FMT( "Previous HIP errors found: {} ({} at {}:{})", msg, hipGetErrorString( err ), file, line ) );
160 #else
161  GEOS_UNUSED_VAR( msg, file, line );
162 #endif
163 }
164 
169 #define GEOS_HYPRE_CHECK_DEVICE_ERRORS( msg ) ::geos::hypre::checkDeviceErrors( msg, __FILE__, __LINE__ )
170 
171 static_assert( sizeof( HYPRE_BigInt ) == sizeof( geos::globalIndex ),
172  "HYPRE_BigInt and geos::globalIndex must have the same size" );
173 
174 static_assert( std::is_signed< HYPRE_BigInt >::value == std::is_signed< geos::globalIndex >::value,
175  "HYPRE_BigInt and geos::globalIndex must both be signed or unsigned" );
176 
177 static_assert( std::is_same< HYPRE_Real, geos::real64 >::value,
178  "HYPRE_Real and geos::real64 must be the same type" );
179 
185 inline HYPRE_BigInt * toHypreBigInt( geos::globalIndex * const index )
186 {
187  return reinterpret_cast< HYPRE_BigInt * >(index);
188 }
189 
195 inline HYPRE_BigInt const * toHypreBigInt( geos::globalIndex const * const index )
196 {
197  return reinterpret_cast< HYPRE_BigInt const * >(index);
198 }
199 
209 HYPRE_Vector parVectorToVectorAll( HYPRE_ParVector const vec );
210 
218 HYPRE_Vector parVectorToVector( HYPRE_ParVector const vec, int const targetRank );
219 
226 HYPRE_Int dummySetup( HYPRE_Solver,
227  HYPRE_ParCSRMatrix,
228  HYPRE_ParVector,
229  HYPRE_ParVector );
230 
239 HYPRE_Int SuperLUDistSolve( HYPRE_Solver solver,
240  HYPRE_ParCSRMatrix A,
241  HYPRE_ParVector b,
242  HYPRE_ParVector x );
243 
249 HYPRE_Int SuperLUDistDestroy( HYPRE_Solver solver );
250 
257 HYPRE_Int relaxationCreate( HYPRE_Solver & solver,
258  HYPRE_Int const type );
259 
268 HYPRE_Int relaxationSetup( HYPRE_Solver solver,
269  HYPRE_ParCSRMatrix A,
270  HYPRE_ParVector b,
271  HYPRE_ParVector x );
272 
281 HYPRE_Int relaxationSolve( HYPRE_Solver solver,
282  HYPRE_ParCSRMatrix A,
283  HYPRE_ParVector b,
284  HYPRE_ParVector x );
285 
291 HYPRE_Int relaxationDestroy( HYPRE_Solver solver );
292 
300 HYPRE_Int chebyshevCreate( HYPRE_Solver & solver,
301  HYPRE_Int const order,
302  HYPRE_Int const eigNumIter );
303 
312 HYPRE_Int chebyshevSetup( HYPRE_Solver solver,
313  HYPRE_ParCSRMatrix A,
314  HYPRE_ParVector b,
315  HYPRE_ParVector x );
316 
325 HYPRE_Int chebyshevSolve( HYPRE_Solver solver,
326  HYPRE_ParCSRMatrix A,
327  HYPRE_ParVector b,
328  HYPRE_ParVector x );
329 
335 HYPRE_Int chebyshevDestroy( HYPRE_Solver solver );
336 
343 {
345  {
348  };
349  return findOption( typeMap, type, "multigrid cycle", "HyprePreconditioner" );
350 }
351 
358 {
360  {
369  };
370  return findOption( typeMap, type, "multigrid relaxation", "HyprePreconditioner" );
371 }
372 
379 {
381  {
393  };
394  return findOption( typeMap, type, "multigrid interpolation", "HyprePreconditioner" );
395 }
396 
403 {
405  {
415  };
416  return findOption( typeMap, type, "multigrid aggressive interpolation", "HyprePreconditioner" );
417 }
418 
425 {
427  {
430  };
431  return findOption( typeMap, type, "ILU", "HyprePreconditioner" );
432 }
433 
440 {
442  {
454  };
455  return findOption( typeMap, type, "multigrid coarse solver", "HyprePreconditioner" );
456 }
457 
464 {
466  {
473  };
474  return findOption( typeMap, type, "multigrid coarsening", "HyprePreconditioner" );
475 }
476 
483 {
485  {
493  };
494  return findOption( typeMap, type, "relaxation", "HyprePreconditioner" );
495 }
496 
503 {
505  {
508  };
509  return findOption( typeMap, type, "ILU", "HyprePreconditioner" );
510 }
511 
516 enum class AMGCoarseningType : HYPRE_Int
517 {
518  CLJP = 0,
519  Ruge_Stueben = 3,
520  Falgout = 6,
521  CLJPDebug = 7,
522  PMIS = 8,
523  PMISDebug = 9,
524  HMIS = 10,
525  CGC = 21,
526  CGC_E = 22
527 };
528 
533 enum class MGRInterpolationType : HYPRE_Int
534 {
535  injection = 0,
536  l1jacobi = 1,
537  jacobi = 2,
539  approximateInverse = 4,
540  blockJacobi = 12
541 };
542 
548 enum class MGRRestrictionType : HYPRE_Int
549 {
550  injection = 0,
551  jacobi = 2,
552  approximateInverse = 3,
553  blockJacobi = 12,
554  cprLike = 13,
555  blockColLumped = 14
556 };
557 
562 enum class MGRCoarseGridMethod : HYPRE_Int
563 {
564  galerkin = 0,
565  nonGalerkin = 1,
567  cprLikeDiag = 2,
569  cprLikeBlockDiag = 3,
571  approximateInverse = 4
573 };
574 
579 enum class MGRFRelaxationType : HYPRE_Int
580 {
581  none = -1,
583  amgVCycle = 2,
587  jacobi = 7,
589  gsElim = 9,
590  l1forwardGaussSeidel = 13,
591  l1backwardGaussSeidel = 14,
592  l1jacobi = 18,
593  gsElimWPivoting = 99,
594  gsElimWInverse = 199
595 };
596 
601 enum class MGRGlobalSmootherType : HYPRE_Int
602 {
603  none = -1,
604  blockJacobi = 0,
605  blockGaussSeidel = 1,
606  jacobi = 2,
607  ilu0 = 16
608 };
609 
610 } // namespace hypre
611 
612 } // namespace geos
613 
614 #endif /*GEOS_LINEARALGEBRA_INTERFACES_HYPREUTILS_HPP_*/
#define GEOS_UNUSED_VAR(...)
Mark an unused variable and silence compiler warnings.
Definition: GeosxMacros.hpp:84
#define GEOS_ERROR_IF(EXP, msg)
Conditionally raise a hard error and terminate the program.
Definition: Logger.hpp:142
Base template for ordered and unordered maps.
HYPRE_Int getAMGInterpolationType(LinearSolverParameters::AMG::InterpType const &type)
Returns hypre's identifier of the AMG interpolation type.
Definition: HypreUtils.hpp:378
MGRCoarseGridMethod
This enum class specifies the strategy for level coarse grid computation in MGR.
Definition: HypreUtils.hpp:563
@ galerkin
Galerkin coarse grid computation using RAP.
MGRRestrictionType
This enum class specifies the strategy for computing the level restriction operator in MGR.
Definition: HypreUtils.hpp:549
@ blockColLumped
Block column-lumped approximation.
@ approximateInverse
Approximate inverse.
@ cprLike
CPR-like restriction.
HYPRE_Int chebyshevSetup(HYPRE_Solver solver, HYPRE_ParCSRMatrix A, HYPRE_ParVector b, HYPRE_ParVector x)
Setup a Chebyshev smoother.
HYPRE_Int getAMGAggressiveInterpolationType(LinearSolverParameters::AMG::AggInterpType const &type)
Returns hypre's identifier of the AMG aggressive interpolation type.
Definition: HypreUtils.hpp:402
HYPRE_Vector parVectorToVectorAll(HYPRE_ParVector const vec)
Gather a parallel vector on every rank.
HYPRE_Int getRelaxationType(LinearSolverParameters::PreconditionerType const type)
Returns hypre's identifier of the relaxation preconditioner type.
Definition: HypreUtils.hpp:482
HYPRE_BigInt * toHypreBigInt(geos::globalIndex *const index)
Converts a non-const array from GEOS globalIndex type to HYPRE_BigInt.
Definition: HypreUtils.hpp:185
HYPRE_Int relaxationSolve(HYPRE_Solver solver, HYPRE_ParCSRMatrix A, HYPRE_ParVector b, HYPRE_ParVector x)
Solve with a relaxation-based smoother.
constexpr HYPRE_MemoryLocation memoryLocation
Memory location used by hypre matrix/vector objects.
Definition: HypreUtils.hpp:117
HYPRE_Int getILUType(LinearSolverParameters::AMG::SmootherType const type)
Returns hypre's identifier of the AMG ILU smoother type.
Definition: HypreUtils.hpp:424
MGRFRelaxationType
This enum class specifies the F-relaxation type.
Definition: HypreUtils.hpp:580
@ amgVCycle
Full AMG VCycle solver.
@ none
no F-relaxation if performed
@ l1forwardGaussSeidel
Gauss-Seidel, forward solve
@ forwardHybridGaussSeidel
hybrid Gauss-Seidel or SOR, forward solve
@ gsElimWInverse
Direct Inversion with Gaussian Elimination (OK for larger systems)
@ backwardHybridGaussSeidel
hybrid Gauss-Seidel or SOR, backward solve
@ hybridSymmetricGaussSeidel
hybrid symmetric Gauss-Seidel or SSOR
@ l1hybridSymmetricGaussSeidel
-scaled hybrid symmetric Gauss-Seidel
@ l1backwardGaussSeidel
Gauss-Seidel, backward solve
@ singleVCycleSmoother
V(1,0) cycle smoother.
@ gsElimWPivoting
Gaussian Elimination with pivoting direct solver (for small systems)
@ gsElim
Gaussian Elimination direct solver (for small systems)
HYPRE_Int getAMGRelaxationType(LinearSolverParameters::AMG::SmootherType const &type)
Returns hypre's identifier of the AMG smoother type.
Definition: HypreUtils.hpp:357
HYPRE_Int dummySetup(HYPRE_Solver, HYPRE_ParCSRMatrix, HYPRE_ParVector, HYPRE_ParVector)
Dummy function that does nothing but conform to hypre's signature for preconditioner setup/apply func...
AMGCoarseningType
This enum class specifies the AMG parallel coarsening algorithm.
Definition: HypreUtils.hpp:517
@ CLJP
Parallel coarsening algorithm using independent sets.
@ Falgout
Uses Ruge_Stueben first, followed by CLJP.
@ Ruge_Stueben
Classical Ruge-Stueben coarsening on each processor.
@ PMIS
Parallel coarsening algorithm using independent sets.
@ CGC_E
Coarsening by M. Griebel, B. Metsch and A. Schweitzer.
@ HMIS
Uses one pass Ruge-Stueben on each processor independently, followed by PMIS.
@ CLJPDebug
Using a fixed random vector, for debugging purposes only.
@ PMISDebug
Using a fixed random vector, for debugging purposes only.
@ CGC
Coarsening by M. Griebel, B. Metsch and A. Schweitzer.
HYPRE_Int relaxationSetup(HYPRE_Solver solver, HYPRE_ParCSRMatrix A, HYPRE_ParVector b, HYPRE_ParVector x)
Setup a relaxation-based smoother.
HYPRE_Int relaxationCreate(HYPRE_Solver &solver, HYPRE_Int const type)
Create a relaxation-based smoother.
constexpr HYPRE_MemoryLocation getMemoryLocation(LvArray::MemorySpace const space)
Definition: HypreUtils.hpp:76
HYPRE_Int getAMGCoarseType(LinearSolverParameters::AMG::CoarseType const &type)
Returns hypre's identifier of the AMG coarse solver type.
Definition: HypreUtils.hpp:439
constexpr LvArray::MemorySpace getLvArrayMemorySpace(HYPRE_MemoryLocation const location)
Definition: HypreUtils.hpp:95
HYPRE_Int chebyshevDestroy(HYPRE_Solver solver)
Destroy a Chebyshev smoother.
HYPRE_Int getAMGCycleType(LinearSolverParameters::AMG::CycleType const &type)
Returns hypre's identifier of the AMG cycle type.
Definition: HypreUtils.hpp:342
HYPRE_Vector parVectorToVector(HYPRE_ParVector const vec, int const targetRank)
Gather a parallel vector onto a single rank.
HYPRE_Int relaxationDestroy(HYPRE_Solver solver)
Destroy a relaxation-based smoother.
HYPRE_Int getAMGCoarseningType(LinearSolverParameters::AMG::CoarseningType const &type)
Returns hypre's identifier of the AMG coarsening type.
Definition: HypreUtils.hpp:463
HYPRE_Int SuperLUDistSolve(HYPRE_Solver solver, HYPRE_ParCSRMatrix A, HYPRE_ParVector b, HYPRE_ParVector x)
The missing wrapper compatible with hypre solver solve signature.
HYPRE_Int SuperLUDistDestroy(HYPRE_Solver solver)
The missing wrapper compatible with hypre solver destroy signature.
void checkDeviceErrors(char const *msg, char const *file, int const line)
Definition: HypreUtils.hpp:150
HYPRE_Int chebyshevSolve(HYPRE_Solver solver, HYPRE_ParCSRMatrix A, HYPRE_ParVector b, HYPRE_ParVector x)
Solve with a Chebyshev smoother.
HYPRE_Int chebyshevCreate(HYPRE_Solver &solver, HYPRE_Int const order, HYPRE_Int const eigNumIter)
Create a Chebyshev smoother.
MGRInterpolationType
This enum class specifies the strategy for computing the level intepolation operator in MGR.
Definition: HypreUtils.hpp:534
@ classicalModifiedInterpolation
Classical modified interpolation.
@ approximateInverse
Approximate inverse.
constexpr LvArray::MemorySpace memorySpace
Memory space used by hypre matrix/vector objects.
Definition: HypreUtils.hpp:115
MGRGlobalSmootherType
This enum class specifies the global smoother type.
Definition: HypreUtils.hpp:602
@ none
no global smoothing is performed (default)
@ ilu0
incomplete LU factorization
parallelDevicePolicy<> execPolicy
Execution policy for operations on hypre data.
Definition: HypreUtils.hpp:113
GEOS_GLOBALINDEX_TYPE globalIndex
Global index type (for indexing objects across MPI partitions).
Definition: DataTypes.hpp:87
GEOS_LOCALINDEX_TYPE localIndex
Local index type (for indexing objects within an MPI partition).
Definition: DataTypes.hpp:84
Container for hypre preconditioner function pointers.
Definition: HypreUtils.hpp:50
HYPRE_Int(*)(HYPRE_Solver) DestroyFunc
Alias for destroy function type.
Definition: HypreUtils.hpp:58
HYPRE_Int(*)(HYPRE_Solver, HYPRE_ParCSRMatrix, HYPRE_ParVector, HYPRE_ParVector) SolveFunc
Alias for apply function type.
Definition: HypreUtils.hpp:55
HYPRE_Int(*)(HYPRE_Solver, HYPRE_ParCSRMatrix, HYPRE_ParVector, HYPRE_ParVector) SetupFunc
Alias for setup function type.
Definition: HypreUtils.hpp:52
HYPRE_Solver ptr
pointer to preconditioner
Definition: HypreUtils.hpp:60
SetupFunc setup
pointer to setup function
Definition: HypreUtils.hpp:61
DestroyFunc destroy
pointer to destroy function
Definition: HypreUtils.hpp:63
SolveFunc solve
pointer to apply function
Definition: HypreUtils.hpp:62
InterpType
AMG interpolation type (HYPRE only)
@ modifiedExtendedI
Modularized extended+i (GPU support)
@ modifiedExtendedE
Modularized extended+e (GPU support)
@ modifiedExtended
Modularized extended classical (GPU support)
@ directBAMG
Direct with separation of weights (GPU support)
@ extended
Extended classical (GPU support)
@ chebyshev
Chebyshev polynomial smoothing.
@ iluk
Incomplete LU with k-level of fill.
@ l1sgs
l1-Symmetric Gauss-Seidel smoothing
@ ilut
Incomplete LU with thresholding.
@ sgs
Symmetric Gauss-Seidel smoothing.
@ bgs
Gauss-Seidel smoothing (backward sweep)
@ fgs
Gauss-Seidel smoothing (forward sweep)
CoarseningType
AMG coarsening types (HYPRE only)
@ CLJP
A parallel coarsening algorithm using independent sets.
@ Falgout
Ruge-Stueben followed by CLJP.
@ RugeStueben
Classical Ruge-Stueben on each processor, followed by a third pass.
@ PMIS
Parallel coarsening as CLJP but with lower complexities (GPU support)
AggInterpType
AMG interpolation type for aggressive coarsening levels (HYPRE only)
@ modifiedExtendedI
Modularized Extended+i (GPU support)
@ extendedIStage2
Extended+i 2-stage (GPU support)
@ modifiedExtendedE
Modularized Extended+e (GPU support)
@ modifiedExtended
Modularized Extended (GPU support)
@ modifiedMultipass
Modularized Multipass (GPU support)
@ extendedStage2
Extended 2-stage (GPU support)
@ jacobi
Jacobi (GPU support in hypre)
@ chebyshev
Chebyshev polynomial (GPU support in hypre)
@ l1jacobi
l1-Jacobi (GPU support in hypre)
@ gsElimWInverse
Direct inverse with Gaussian Elimination.
@ direct
Direct solver as preconditioner.
@ gsElimWPivoting
Gaussian Elimination with pivoting direct solver.
@ bgs
Gauss-Seidel smoothing (backward sweep)
@ chebyshev
Chebyshev polynomial smoothing.
@ iluk
Incomplete LU with k-level of fill.
@ l1sgs
l1-Symmetric Gauss-Seidel smoothing
@ ilut
Incomplete LU with thresholding.
@ sgs
Symmetric Gauss-Seidel smoothing.
@ bgs
Gauss-Seidel smoothing (backward sweep)
@ fgs
Gauss-Seidel smoothing (forward sweep)