Add linear solvers

RELEASE_NOTES.md

@@ -77,6 +77,7 @@ Notable changes include:
     * Both ASPH and ASPHClassic now allow the user to override the final H evolution through optional functors added to the classes:
       - HidealFilter
       - RadialFunctor
+    * Added linear solvers from Eigen and Hypre.
     * FacetedSurfaceASPHHydro has been removed in favor of providing user filters to the ASPH methods (i.e., the RadialFunctor method).
     * Field resizing operations have been removed from the public interface.
     * Performance analysis tools are improved.

cmake/InstallTPLs.cmake

@@ -194,7 +194,7 @@ set_property(GLOBAL PROPERTY SPHERAL_FP_DIRS ${SPHERAL_FP_DIRS})
 
 message("-----------------------------------------------------------------------------")
 # Use find_package to get Sundials
-if (ENABLE_SUNDIALS)
+if (SPHERAL_ENABLE_SOLVERS)
   set(SUNDIALS_DIR "${sundials_DIR}")
   find_package(SUNDIALS REQUIRED NO_DEFAULT_PATH
     COMPONENTS kinsol nvecparallel nvecmpiplusx nvecserial 
@@ -205,6 +205,7 @@ if (ENABLE_SUNDIALS)
     list(APPEND SPHERAL_FP_DIRS ${sundials_DIR})
     message("Found SUNDIALS External Package")
   endif()
+  list(APPEND SPHERAL_EXTERN_LIBS hypre)
 endif()
 
 message("-----------------------------------------------------------------------------")

cmake/spheral/SpheralHandleTPL.cmake

@@ -86,6 +86,9 @@ function(Spheral_Handle_TPL lib_name TPL_CMAKE_DIR)
     INCLUDES ${${lib_name}_INCLUDE_DIR}
     LIBRARIES ${${lib_name}_LIBRARIES}
     EXPORTABLE ON)
+  if(${lib_name}_EXT_LIBRARIES)
+    target_link_libraries(${lib_name} INTERFACE ${${lib_name}_EXT_LIBRARIES})
+  endif()
   get_target_property(_is_imported ${lib_name} IMPORTED)
   if(NOT ${_is_imported})
     install(TARGETS ${lib_name}

cmake/tpl/hypre.cmake

@@ -0,0 +1 @@
+set(hypre_libs libHYPRE.a)

scripts/spack/packages/spheral/package.py

@@ -40,7 +40,7 @@ class Spheral(CachedCMakePackage, CudaPackage, ROCmPackage):
     variant('caliper', default=True, description='Enable Caliper timers.')
     variant('opensubdiv', default=True, description='Enable use of opensubdiv to do refinement.')
     variant('network', default=True, description='Disable to build Spheral from a local buildcache.')
-    variant('sundials', default=True, when="+mpi", description='Build Sundials package.')
+    variant('solvers', default=True, when="+mpi", description='Build Sundials and Hypre packages.')
     variant('leos', default=LEOSpresent, when="+mpi", description='Build LEOS package.')
 
     # -------------------------------------------------------------------------
@@ -85,7 +85,9 @@ class Spheral(CachedCMakePackage, CudaPackage, ROCmPackage):
     depends_on('polytope +python', type='build', when="+python")
     depends_on('polytope ~python', type='build', when="~python")
 
-    depends_on('sundials@7.0.0 ~shared cxxstd=17 cppflags="-fPIC"', type='build', when='+sundials')
+    with when("+solvers"):
+        depends_on('sundials@7.0.0 ~shared cxxstd=17 cppflags="-fPIC"', type='build')
+        depends_on('hypre@2.26.0 ~shared cppflags="-fPIC" cflags="-fPIC"', type='build')
 
     depends_on('leos@8.4.2', type='build', when='+leos')
 
@@ -265,9 +267,13 @@ def initconfig_package_entries(self):
         if "+python" in spec:
             entries.append(cmake_cache_path('python_DIR', spec['python'].prefix))
 
-        if spec.satisfies("+sundials"):
+        if spec.satisfies("+solvers"):
             entries.append(cmake_cache_path('sundials_DIR', spec['sundials'].prefix))
-            entries.append(cmake_cache_option('ENABLE_SUNDIALS', True))
+            entries.append(cmake_cache_path('hypre_DIR', spec['hypre'].prefix))
+            entries.append(cmake_cache_path('hypre_INCLUDES', spec['hypre'].prefix.include))
+            hypre_libs = spec["blas"].libs + spec["lapack"].libs
+            entries.append(cmake_cache_path('hypre_EXT_LIBRARIES', hypre_libs.joined(";")))
+            entries.append(cmake_cache_option('SPHERAL_ENABLE_SOLVERS', True))
 
         if spec.satisfies("+leos"):
             entries.append(cmake_cache_path('leos_DIR', spec['leos'].prefix))

src/CMakeLists.txt

@@ -69,7 +69,7 @@ if (SPHERAL_ENABLE_SVPH)
   list(APPEND _packages SVPH)
 endif()
 
-if (ENABLE_SUNDIALS)
+if (SPHERAL_ENABLE_SOLVERS)
   list(APPEND _packages Solvers)
 endif()
 

src/PYB11/CMakeLists.txt

@@ -73,7 +73,7 @@ if (SPHERAL_ENABLE_SVPH)
   list(APPEND _python_packages SVPH)
 endif()
 
-if (ENABLE_SUNDIALS)
+if (SPHERAL_ENABLE_SOLVERS)
   list(PREPEND _python_packages Solvers)
 endif()
 

src/PYB11/Solvers/EigenLinearSolver.py

@@ -0,0 +1,10 @@
+from PYB11Generator import *
+from LinearSolver import *
+
+@PYB11holder("std::shared_ptr")
+class EigenLinearSolver(LinearSolver):
+    def pyinit(self,
+               options = "std::shared_ptr<EigenOptions>"):
+        "Eigen solver"
+
+PYB11inject(LinearSolverAbstractMethods, EigenLinearSolver, virtual = True)

src/PYB11/Solvers/EigenOptions.py

@@ -0,0 +1,8 @@
+from PYB11Generator import *
+
+@PYB11holder("std::shared_ptr")
+class EigenOptions:
+    def pyinit(self):
+        "Holds the options for Eigen solvers"
+
+    qr = PYB11readwrite()

src/PYB11/Solvers/HypreLinearSolver.py

@@ -0,0 +1,10 @@
+from PYB11Generator import *
+from LinearSolver import *
+
+@PYB11holder("std::shared_ptr")
+class HypreLinearSolver(LinearSolver):
+    def pyinit(self,
+               options = "std::shared_ptr<HypreOptions>"):
+        "Hypre solver"
+
+PYB11inject(LinearSolverAbstractMethods, HypreLinearSolver, virtual = True)

src/PYB11/Solvers/HypreOptions.py

@@ -0,0 +1,62 @@
+from PYB11Generator import *
+
+@PYB11holder("std::shared_ptr")
+class HypreOptions:
+    def pyinit(self):
+        "Holds the options for Hypre preconditioners and solver"
+
+    quitIfDiverged = PYB11readwrite()
+
+    logLevel = PYB11readwrite()
+    printLevel = PYB11readwrite()
+
+    sumDuplicates = PYB11readwrite()
+    addToValues = PYB11readwrite()
+
+    saveLinearSystem = PYB11readwrite()
+    printIterations = PYB11readwrite()
+    maxNumberOfIterations = PYB11readwrite()
+    toleranceL2 = PYB11readwrite()
+    useRobustTolerance = PYB11readwrite()
+    absoluteTolerance = PYB11readwrite()
+
+    kDim = PYB11readwrite()
+    minIters = PYB11readwrite()
+
+    HyprePreconditionerType = PYB11enum(("NoPreconditioner",
+                                         "AMGPreconditioner",
+                                         "ILUPreconditioner"),
+                                         export_values = True,
+                                         doc = "Preconditioners available for Hypre")
+    preconditionerType = PYB11readwrite()
+
+    measure_type = PYB11readwrite()
+    pcTol = PYB11readwrite()
+    minItersAMG = PYB11readwrite()
+    maxItersAMG = PYB11readwrite()
+    cycleType = PYB11readwrite()
+    coarsenTypeAMG = PYB11readwrite()
+    strongThresholdAMG = PYB11readwrite()
+    maxRowSumAMG = PYB11readwrite()
+    interpTypeAMG = PYB11readwrite()
+    aggNumLevelsAMG = PYB11readwrite()
+    aggInterpTypeAMG = PYB11readwrite()
+    pMaxElmtsAMG = PYB11readwrite()
+    truncFactorAMG = PYB11readwrite()
+    printLevelAMG = PYB11readwrite()
+    logLevelAMG = PYB11readwrite()
+    relaxWeightAMG = PYB11readwrite()
+    relaxTypeAMG = PYB11readwrite()
+    relaxTypeCoarseAMG = PYB11readwrite()
+    cycleNumSweepsAMG = PYB11readwrite()
+    cycleNumSweepsCoarseAMG = PYB11readwrite()
+    maxLevelsAMG = PYB11readwrite()
+
+    useILUT = PYB11readwrite()
+    factorLevelILU = PYB11readwrite()
+    rowScaleILU = PYB11readwrite()
+    printLevelILU = PYB11readwrite()
+    dropToleranceILU = PYB11readwrite()
+
+    zeroOutNegativities = PYB11readwrite()
+

src/PYB11/Solvers/IncrementalStatistic.py

@@ -0,0 +1,46 @@
+from PYB11Generator import *
+
+@PYB11template("DataType")
+@PYB11holder("std::shared_ptr")
+class IncrementalStatistic:
+    def pyinit(self,
+               meanGuess = ("const %(DataType)s", 0),
+               name = ("const std::string", "\"IncrementalStatistic\""),
+               printAdd = ("bool", "false")):
+        "Allows tracking of incrementally added statistics"
+
+    def add(self,
+            data = "const %(DataType)s"):
+        return "void"
+
+    @PYB11const
+    def mean(self):
+        return "%(DataType)s"
+    @PYB11const
+    def variance(self):
+        return "%(DataType)s"
+    @PYB11const
+    def total(self):
+        return "%(DataType)s"
+    @PYB11const
+    def min(self):
+        return "%(DataType)s"
+    @PYB11const
+    def max(self):
+        return "%(DataType)s"
+    @PYB11const
+    def meanGuess(self):
+        return "%(DataType)s"
+    @PYB11const
+    def numPoints(self):
+        return "%(DataType)s"
+    @PYB11const
+    def shiftedSum(self):
+        return "%(DataType)s"
+    @PYB11const
+    def shiftedSum2(self):
+        return "%(DataType)s"
+    @PYB11cppname("print")
+    @PYB11const
+    def print1(self):
+        return "void"

src/PYB11/Solvers/LinearSolver.py

@@ -0,0 +1,164 @@
+from PYB11Generator import *
+
+@PYB11ignore
+class LinearSolverAbstractMethods:
+    @PYB11cppname("initialize")
+    def initializePtr(self,
+                    numLocVal = "const unsigned",
+                    firstGlobInd = "const unsigned",
+                    numValsPerRow = "const unsigned*"):
+        "Create matrices and vectors"
+        return "void"
+
+    def beginFill(self):
+        "Begin matrix fill"
+        return "void"
+
+    @PYB11cppname("setMatRow")
+    def setMatRowPtr(self,
+                     numVals = "const unsigned",
+                     globRowInd = "const unsigned",
+                     globColInd = "const unsigned*",
+                     colVal = "const double*"):
+        "Set values in matrix"
+        return "void"
+
+    @PYB11cppname("setMatRows")
+    def setMatRowsPtr(self,
+                      numRows = "const unsigned",
+                      numColsPerRow = "const unsigned*",
+                      globRowInd = "const unsigned*",
+                      globColInd = "const unsigned*",
+                      colVal = "const double*"):
+        "Set multiple rows in matrix"
+        return "void"
+
+    def assemble(self):
+        "Assemble matrix after fill"
+        return "void"
+
+    def finalize(self):
+        "Create solver/preconditioner"
+        return "void"
+
+    @PYB11cppname("set")
+    def setPtr(self,
+               c = "const LinearSolver::Component",
+               numVals = "const unsigned",
+               firstGlobInd = "const unsigned",
+               val = "const double*"):
+        "Set values from external vector"
+        return "void"
+
+    @PYB11cppname("get")
+    def getPtr(self,
+               c = "const LinearSolver::Component",
+               numVals = "const unsigned",
+               firstGlobInd = "const unsigned",
+               val = "double*"):
+        "Get values from internal vector"
+        return "void"
+
+    def solve(self):
+        "Solve the system of equations in place"
+        return "void"
+
+    def multiply(self):
+        "Multiply the matrix by a vector"
+        return "void"
+
+@PYB11holder("std::shared_ptr")
+class LinearSolver:
+    Component = PYB11enum(("LHS",
+                           "RHS"),
+                          export_values = True,
+                          doc = "Choose which vector to use when applying a function")
+
+    def pyinit(self):
+        "Pure virtual linear solver class"
+
+    @PYB11implementation("""[](LinearSolver& self,
+                               const unsigned numLocVal,
+                               const unsigned firstGlobInd,
+                               const std::vector<unsigned>& numValsPerRow) {
+                               self.initialize(numLocVal, firstGlobInd, &numValsPerRow[0]); }""")
+    @PYB11pyname("initialize")
+    def initializeVec(self,
+                      numLocVal = "const unsigned",
+                      firstGlobInd = "const unsigned",
+                      numValsPerRow = "const std::vector<unsigned>&"):
+        "Create matrices and vectors"
+        return "void"
+
+    @PYB11implementation("""[](LinearSolver& self,
+                               const unsigned numVals,
+                               const unsigned globRowInd,
+                               const std::vector<unsigned>& globColInd,
+                               const std::vector<double>& colVal) {
+                               self.setMatRow(numVals, globRowInd, &globColInd[0], &colVal[0]); }""")
+    @PYB11pyname("setMatRow")
+    def setMatRowVec(self,
+                     numVals = "const unsigned",
+                     globRowInd = "const unsigned",
+                     globColInd = "const vector<unsigned>&",
+                     colVal = "const std::vector<double>&"):
+        "Set values in matrix"
+        return "void"
+
+    @PYB11implementation("""[](LinearSolver& self,
+                               const unsigned numRows,
+                               const std::vector<unsigned>& numColsPerRow,
+                               const std::vector<unsigned>& globRowInd,
+                               const std::vector<unsigned>& globColInd,
+                               const std::vector<double>& colVal) {
+                               self.setMatRows(numRows, &numColsPerRow[0], &globRowInd[0], &globColInd[0], &colVal[0]); }""")
+    @PYB11pyname("setMatRows")
+    def setMatRowsVec(self,
+                      numRows = "const unsigned",
+                      numColsPerRow = "const vector<unsigned>&",
+                      globRowInd = "const vector<unsigned>&",
+                      globColInd = "const vector<unsigned>&",
+                      colVal = "const std::vector<double>&"):
+        "Set multiple rows in matrix"
+        return "void"
+
+    @PYB11implementation("""[](LinearSolver& self,
+                               const LinearSolver::Component c,
+                               const unsigned numVals,
+                               const unsigned firstGlobInd,
+                               const std::vector<double>& val) {
+                               self.set(c, numVals, firstGlobInd, &val[0]); }""")
+    @PYB11pyname("set")
+    def setVec(self,
+               c = "const LinearSolver::Component",
+               numVals = "const unsigned",
+               firstGlobInd = "const unsigned",
+               val = "const std::vector<double>&"):
+        "Set values in RHS"
+        return "void"
+
+    @PYB11implementation("""[](LinearSolver& self,
+                               const LinearSolver::Component c,
+                               const unsigned numVals,
+                               const unsigned firstGlobInd) {
+                               std::vector<double> result(numVals);
+                               self.get(c, numVals, firstGlobInd, &result[0]);
+                               return result; }""")
+    @PYB11pyname("get")
+    def getVec(self,
+               c = "const LinearSolver::Component",
+               numVals = "const unsigned",
+               firstGlobInd = "const unsigned"):
+        "Get values from RHS"
+        return "std::vector<double>"
+
+    description = PYB11property(returnType = "std::string",
+                                getter = "getDescription",
+                                setter = "setDescription")
+
+    @PYB11virtual
+    @PYB11const
+    def statistics(self):
+        return "std::vector<std::shared_ptr<IncrementalStatistic<double>>>"
+
+PYB11inject(LinearSolverAbstractMethods, LinearSolver, pure_virtual = True)