diff --git a/RELEASE_NOTES.md b/RELEASE_NOTES.md index 900c4faeb..5a16b5903 100644 --- a/RELEASE_NOTES.md +++ b/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. diff --git a/cmake/InstallTPLs.cmake b/cmake/InstallTPLs.cmake index a3f566334..852a5fa21 100644 --- a/cmake/InstallTPLs.cmake +++ b/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("-----------------------------------------------------------------------------") diff --git a/cmake/spheral/SpheralHandleTPL.cmake b/cmake/spheral/SpheralHandleTPL.cmake index 3e75a8978..d4f1d8f1f 100644 --- a/cmake/spheral/SpheralHandleTPL.cmake +++ b/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} diff --git a/cmake/tpl/hypre.cmake b/cmake/tpl/hypre.cmake new file mode 100644 index 000000000..4bd67476f --- /dev/null +++ b/cmake/tpl/hypre.cmake @@ -0,0 +1 @@ +set(hypre_libs libHYPRE.a) diff --git a/scripts/spack/packages/spheral/package.py b/scripts/spack/packages/spheral/package.py index 3afe13681..419151c33 100644 --- a/scripts/spack/packages/spheral/package.py +++ b/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)) diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 4b5997c53..fd5e4f64c 100644 --- a/src/CMakeLists.txt +++ b/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() diff --git a/src/PYB11/CMakeLists.txt b/src/PYB11/CMakeLists.txt index 0ebeaf762..3f9ddac75 100644 --- a/src/PYB11/CMakeLists.txt +++ b/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() diff --git a/src/PYB11/Solvers/EigenLinearSolver.py b/src/PYB11/Solvers/EigenLinearSolver.py new file mode 100644 index 000000000..10fdd5e0f --- /dev/null +++ b/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"): + "Eigen solver" + +PYB11inject(LinearSolverAbstractMethods, EigenLinearSolver, virtual = True) diff --git a/src/PYB11/Solvers/EigenOptions.py b/src/PYB11/Solvers/EigenOptions.py new file mode 100644 index 000000000..1153010d0 --- /dev/null +++ b/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() diff --git a/src/PYB11/Solvers/HypreLinearSolver.py b/src/PYB11/Solvers/HypreLinearSolver.py new file mode 100644 index 000000000..2d296746e --- /dev/null +++ b/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"): + "Hypre solver" + +PYB11inject(LinearSolverAbstractMethods, HypreLinearSolver, virtual = True) diff --git a/src/PYB11/Solvers/HypreOptions.py b/src/PYB11/Solvers/HypreOptions.py new file mode 100644 index 000000000..44f801eeb --- /dev/null +++ b/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() + diff --git a/src/PYB11/Solvers/IncrementalStatistic.py b/src/PYB11/Solvers/IncrementalStatistic.py new file mode 100644 index 000000000..242ef2e2f --- /dev/null +++ b/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" diff --git a/src/PYB11/Solvers/LinearSolver.py b/src/PYB11/Solvers/LinearSolver.py new file mode 100644 index 000000000..8f1d61056 --- /dev/null +++ b/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& numValsPerRow) { + self.initialize(numLocVal, firstGlobInd, &numValsPerRow[0]); }""") + @PYB11pyname("initialize") + def initializeVec(self, + numLocVal = "const unsigned", + firstGlobInd = "const unsigned", + numValsPerRow = "const std::vector&"): + "Create matrices and vectors" + return "void" + + @PYB11implementation("""[](LinearSolver& self, + const unsigned numVals, + const unsigned globRowInd, + const std::vector& globColInd, + const std::vector& colVal) { + self.setMatRow(numVals, globRowInd, &globColInd[0], &colVal[0]); }""") + @PYB11pyname("setMatRow") + def setMatRowVec(self, + numVals = "const unsigned", + globRowInd = "const unsigned", + globColInd = "const vector&", + colVal = "const std::vector&"): + "Set values in matrix" + return "void" + + @PYB11implementation("""[](LinearSolver& self, + const unsigned numRows, + const std::vector& numColsPerRow, + const std::vector& globRowInd, + const std::vector& globColInd, + const std::vector& colVal) { + self.setMatRows(numRows, &numColsPerRow[0], &globRowInd[0], &globColInd[0], &colVal[0]); }""") + @PYB11pyname("setMatRows") + def setMatRowsVec(self, + numRows = "const unsigned", + numColsPerRow = "const vector&", + globRowInd = "const vector&", + globColInd = "const vector&", + colVal = "const std::vector&"): + "Set multiple rows in matrix" + return "void" + + @PYB11implementation("""[](LinearSolver& self, + const LinearSolver::Component c, + const unsigned numVals, + const unsigned firstGlobInd, + const std::vector& 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&"): + "Set values in RHS" + return "void" + + @PYB11implementation("""[](LinearSolver& self, + const LinearSolver::Component c, + const unsigned numVals, + const unsigned firstGlobInd) { + std::vector 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" + + description = PYB11property(returnType = "std::string", + getter = "getDescription", + setter = "setDescription") + + @PYB11virtual + @PYB11const + def statistics(self): + return "std::vector>>" + +PYB11inject(LinearSolverAbstractMethods, LinearSolver, pure_virtual = True) diff --git a/src/PYB11/Solvers/Solvers_PYB11.py b/src/PYB11/Solvers/Solvers_PYB11.py index 5ba56323d..cbcb071da 100644 --- a/src/PYB11/Solvers/Solvers_PYB11.py +++ b/src/PYB11/Solvers/Solvers_PYB11.py @@ -5,13 +5,23 @@ """ from PYB11Generator import * +from SpheralCommon import * +from spheralDimensions import * +dims = spheralDimensions() -PYB11opaque = ["std::vector"] +# PYB11opaque = ["std::vector"] #------------------------------------------------------------------------------- # Includes #------------------------------------------------------------------------------- -PYB11includes = ['"Solvers/SolverFunction.hh"', +PYB11includes = ['"Solvers/containsConstantNodes.hh"', + '"Solvers/EigenLinearSolver.hh"', + '"Solvers/EigenOptions.hh"', + '"Solvers/HypreLinearSolver.hh"', + '"Solvers/HypreOptions.hh"', + '"Solvers/IncrementalStatistic.hh"', + '"Solvers/LinearSolver.hh"', + '"Solvers/SolverFunction.hh"', '"Solvers/KINSOL.hh"'] #------------------------------------------------------------------------------- @@ -19,5 +29,39 @@ #------------------------------------------------------------------------------- PYB11namespaces = ["Spheral"] -from SolverFunction import * +from EigenLinearSolver import * +from EigenOptions import * +from HypreLinearSolver import * +from HypreOptions import * +from IncrementalStatistic import * from KINSOL import * +from LinearSolver import * +from SolverFunction import * + +#------------------------------------------------------------------------------- +# Instantiations +#------------------------------------------------------------------------------- +for ndim in dims: + # Dimension-dependent + exec(''' +@PYB11pycppname("containsConstantNodes") +def containsConstantNodes1%(ndim)id(boundary = "const Boundary>*"): + "Does this boundary contain constant boundary nodes?" + return "bool" +@PYB11pycppname("containsConstantNodes") +def containsConstantNodes2%(ndim)id(boundaries = "const std::vector>*>&"): + "Do these boundaries contain constant boundary nodes?" + return "bool" +''' % {"ndim" : ndim}) + +#------------------------------------------------------------------------------- +# Scalar instantiations +#------------------------------------------------------------------------------- +scalar_types = (("int", "Ordinal"), + ("double", "Scalar")) + +for (value, label) in scalar_types: + exec(''' +%(label)sIncrementalStatistic = PYB11TemplateClass(IncrementalStatistic, template_parameters=("%(value)s")) +''' % {"value" : value, + "label" : label}) diff --git a/src/Solvers/CMakeLists.txt b/src/Solvers/CMakeLists.txt index acbcfc17f..cb13c072a 100644 --- a/src/Solvers/CMakeLists.txt +++ b/src/Solvers/CMakeLists.txt @@ -2,14 +2,30 @@ include_directories(.) set(Solvers_inst ) set(Solvers_sources + EigenLinearSolver.cc + HypreFunctions.cc + HypreLinearSolver.cc KINSOL.cc + LinearSolver.cc ) instantiate(Solvers_inst Solvers_sources) set(Solvers_headers - SolverFunction.hh + containsConstantNodes.hh + containsConstantNodesInline.hh + EigenLinearSolver.hh + EigenOptions.hh + HypreFunctions.hh + HypreLinearSolver.hh + HypreLinearSolverInline.hh + HypreOptions.hh + IncrementalStatistic.hh + IncrementalStatisticInline.hh KINSOL.hh + LinearSolver.hh + LinearSolverInline.hh + SolverFunction.hh ) spheral_add_obj_library(Solvers SPHERAL_OBJ_LIBS) diff --git a/src/Solvers/EigenLinearSolver.cc b/src/Solvers/EigenLinearSolver.cc new file mode 100644 index 000000000..56c65d573 --- /dev/null +++ b/src/Solvers/EigenLinearSolver.cc @@ -0,0 +1,176 @@ +//---------------------------------Spheral++----------------------------------// +// EigenLinearSolver +// +// Represents a solver for a linear system of equations +//----------------------------------------------------------------------------// + +#include "EigenLinearSolver.hh" + +#include +#include // for std::accumulate +#include "Utilities/DBC.hh" +#include "Distributed/Process.hh" + +namespace Spheral { + +//------------------------------------------------------------------------------ +// Constructor +//------------------------------------------------------------------------------ +EigenLinearSolver:: +EigenLinearSolver(std::shared_ptr options) : + mGraphChangedSinceFill(true), + mMatrixChangedSinceFactorization(true), + mOptions(options) { + this->setDescription("EigenLinearSolver"); +} + +//------------------------------------------------------------------------------ +// Initialize the matrices and vectors +//------------------------------------------------------------------------------ +void +EigenLinearSolver:: +initialize(const unsigned numLocVal, + const unsigned firstGlobInd, + const unsigned* numValsPerRow) { + VERIFY(Process::getTotalNumberOfProcesses() == 1); + + mMatrix.setZero(); + mMatrix.resize(numLocVal, numLocVal); + mRhs.resize(numLocVal); + mLhs.resize(numLocVal); + + mGraphChangedSinceFill = true; +} + +//------------------------------------------------------------------------------ +// Set matrix values +//------------------------------------------------------------------------------ +void +EigenLinearSolver:: +setMatRow(const unsigned numVals, + const unsigned globRowInd, + const unsigned* globColInd, + const double* colVal) { + VERIFY(globRowInd < mMatrix.rows()); + for (auto j = 0u; j < numVals; ++j) { + mMatrix.coeffRef(globRowInd, globColInd[j]) = colVal[j]; + } + mMatrixChangedSinceFactorization = true; +} + +void +EigenLinearSolver:: +setMatRows(const unsigned numRows, + const unsigned* numColsPerRow, + const unsigned* globRowInd, + const unsigned* globColInd, + const double* colVal) { + VERIFY(numRows > 0); + auto index = 0u; + for (auto i = 0u; i < numRows; ++i) { + setMatRow(numColsPerRow[i], globRowInd[i], globColInd + index, colVal + index); + index += numColsPerRow[i]; + } + mMatrixChangedSinceFactorization = true; +} + +//------------------------------------------------------------------------------ +// Assemble matrix after fill +//------------------------------------------------------------------------------ +void +EigenLinearSolver:: +assemble() { + mMatrix.makeCompressed(); + + mGraphChangedSinceFill = false; + mMatrixChangedSinceFactorization = true; +} + +//------------------------------------------------------------------------------ +// Create solver and preconditioner +//------------------------------------------------------------------------------ +void +EigenLinearSolver:: +finalize() { + VERIFY(!mGraphChangedSinceFill); + + if (mOptions->qr) { + mQRSolver.analyzePattern(mMatrix); + mQRSolver.factorize(mMatrix); + } + else { + mLUSolver.analyzePattern(mMatrix); + mLUSolver.factorize(mMatrix); + } + + mMatrixChangedSinceFactorization = false; +} + +//------------------------------------------------------------------------------ +// Set/get values in LHS and RHS +//------------------------------------------------------------------------------ +void +EigenLinearSolver:: +set(const Component c, + const unsigned numVals, + const unsigned firstGlobInd, + const double* val) { + switch (c) { + case RHS: + for (auto i = 0u; i < numVals; ++i) { + mRhs(firstGlobInd + i) = val[i]; + } + break; + case LHS: + for (auto i = 0u; i < numVals; ++i) { + mLhs(firstGlobInd + i) = val[i]; + } + break; + } +} + +void +EigenLinearSolver:: +get(const Component c, + const unsigned numVals, + const unsigned firstGlobInd, + double* val) { + switch (c) { + case RHS: + for (auto i = 0u; i < numVals; ++i) { + val[i] = mRhs(firstGlobInd + i); + } + break; + case LHS: + for (auto i = 0u; i < numVals; ++i) { + val[i] = mLhs(firstGlobInd + i); + } + break; + } +} + +//------------------------------------------------------------------------------ +// Solve the system of equations in place +//------------------------------------------------------------------------------ +void +EigenLinearSolver:: +solve() { + VERIFY(!mMatrixChangedSinceFactorization); + if (mOptions->qr) { + mLhs = mQRSolver.solve(mRhs); + } + else { + mLhs = mLUSolver.solve(mRhs); + } +} + +//------------------------------------------------------------------------------ +// Multiply by the matrix +//------------------------------------------------------------------------------ +void +EigenLinearSolver:: +multiply() { + mLhs = mMatrix * mRhs; +} + +} // end namespace Spheral diff --git a/src/Solvers/EigenLinearSolver.hh b/src/Solvers/EigenLinearSolver.hh new file mode 100644 index 000000000..92fd54f37 --- /dev/null +++ b/src/Solvers/EigenLinearSolver.hh @@ -0,0 +1,96 @@ +//---------------------------------Spheral++----------------------------------// +// EigenLinearSolver +// +// Solver linear system using Eigen +// Only works on one processor +//----------------------------------------------------------------------------// +#ifndef __Spheral_EigenLinearSolver_hh__ +#define __Spheral_EigenLinearSolver_hh__ + +#include +#include + +#include "Eigen/OrderingMethods" +#include "Eigen/Sparse" +#include "Eigen/SparseLU" +#include "Eigen/SparseQR" + +#include "EigenOptions.hh" +#include "LinearSolver.hh" + +namespace Spheral { + +class EigenLinearSolver : public LinearSolver { +public: + typedef Eigen::SparseMatrix MatrixType; + typedef Eigen::Triplet DataType; + typedef Eigen::VectorXd VectorType; + typedef Eigen::SparseLU> LUSolverType; + typedef Eigen::SparseQR> QRSolverType; + + using LinearSolver::Component; + + // Constructor + EigenLinearSolver(std::shared_ptr options); + + // Create matrices and vectors + virtual void initialize(const unsigned numLocVal, + const unsigned firstGlobInd, + const unsigned* numValsPerRow) override; + + // Begin the matrix fill + virtual void beginFill() override {} + + // Set values in matrix + virtual void setMatRow(const unsigned numVals, // Number of values in this row + const unsigned globRowInd, // Global index of this row + const unsigned* globColInd, // [numVals] Global indices for columns + const double* colVal) override; // [numVals] Values + virtual void setMatRows(const unsigned numRows, // Number of rows + const unsigned* numColsPerRow, // [numRows] Number of columns for each of these rows + const unsigned* globRowInd, // [numRows] Global indices for these rows + const unsigned* globColInd, // [numColsPerRow[numRows]] Global column indies for each value + const double* colVal) override; // [numColsPerRow[numRows]] Values + + // Assemble matrix after fill + virtual void assemble() override; + + // Create solver/preconditioner + virtual void finalize() override; + + // Set/get values + virtual void set(const Component c, + const unsigned numVals, // Number of vector values to set + const unsigned firstGlobInd, // First global index for values + const double* val) override; // [numVals] Values + virtual void get(const Component c, + const unsigned numVals, // Number of vector values to set + const unsigned firstGlobInd, // First global index for values + double* val) override; // [numVals] Values + + // Solve the system of equations in place + virtual void solve() override; + + // Multiply the matrix by a vector + virtual void multiply() override; + +private: + // Have we initialized things correctly? + bool mGraphChangedSinceFill; + bool mMatrixChangedSinceFactorization; + + // Input data + std::shared_ptr mOptions; + + // Eigen data + MatrixType mMatrix; + LUSolverType mLUSolver; + QRSolverType mQRSolver; + VectorType mRhs; + VectorType mLhs; + +}; // end class EigenLinearSolver + +} // end namespace Spheral + +#endif diff --git a/src/Solvers/EigenOptions.hh b/src/Solvers/EigenOptions.hh new file mode 100644 index 000000000..2da74277f --- /dev/null +++ b/src/Solvers/EigenOptions.hh @@ -0,0 +1,20 @@ +//---------------------------------Spheral++----------------------------------// +// EigenOptions +// +// Holds the options for Eigen solvers and preconditioners +//----------------------------------------------------------------------------// +#ifndef __Spheral_EigenOptions_hh__ +#define __Spheral_EigenOptions_hh__ + +namespace Spheral { + +struct EigenOptions { + EigenOptions() { } + + bool qr = false; + +}; // end struct EigenOptions + +} // end namespace Spheral + +#endif diff --git a/src/Solvers/HypreFunctions.cc b/src/Solvers/HypreFunctions.cc new file mode 100644 index 000000000..cbe77c4a7 --- /dev/null +++ b/src/Solvers/HypreFunctions.cc @@ -0,0 +1,693 @@ +//---------------------------------Spheral++----------------------------------// +// HypreFunctions +// +// Convenience functions for Hypre +//----------------------------------------------------------------------------// + +#include "HypreFunctions.hh" + +#include +#include // for std::iota + +#include "_hypre_parcsr_ls.h" +#include "HYPRE.h" +#include "HYPRE_IJ_mv.h" + +#include "Distributed/Communicator.hh" +#include "Solvers/HypreOptions.hh" +#include "Solvers/IncrementalStatistic.hh" +#include "Utilities/DBC.hh" +#include "Utilities/Timer.hh" + +static_assert(sizeof(int) == sizeof(HYPRE_Int), + "if this fails, need to refactor to convert to HYPRE_Int"); +static_assert(sizeof(int) == sizeof(HYPRE_BigInt), + "if this fails, need to refactor to convert to HYPRE_BigInt"); + +namespace Spheral { + +//------------------------------------------------------------------------------ +// Create Hypre vector that should self-destruct when it goes out of scope. +//------------------------------------------------------------------------------ +std::shared_ptr +HypreFunctions:: +createVector(const unsigned numLocVal, + const unsigned firstGlobInd) { + TIME_FUNCTION; + const auto lastGlobInd = static_cast(firstGlobInd + numLocVal) - 1; + + int hypreStatus; + VectorType* tempVector; + + hypreStatus = HYPRE_IJVectorCreate(Communicator::communicator(), + firstGlobInd, + lastGlobInd, + &tempVector); + VERIFY2(hypreStatus == 0, + "HYPRE_IJVectorCreate, error code " + std::to_string(hypreStatus)); + + hypreStatus = HYPRE_IJVectorSetObjectType(tempVector, HYPRE_PARCSR); + VERIFY2(hypreStatus == 0, + "HYPRE_IJVectorSetObjectType, error code " + std::to_string(hypreStatus)); + + hypreStatus = HYPRE_IJVectorInitialize(tempVector); + VERIFY2(hypreStatus == 0, + "HYPRE_IJVectorInitialize, error code " + std::to_string(hypreStatus)); + + return std::shared_ptr(tempVector, HYPRE_IJVectorDestroy); +} + +//------------------------------------------------------------------------------ +// Create Hypre matrix that should self-destruct when it goes out of scope. +//------------------------------------------------------------------------------ +std::shared_ptr +HypreFunctions:: +createMatrix(const unsigned numLocVal, + const unsigned firstGlobInd, + const unsigned* numValsPerRow) { + TIME_FUNCTION; + const auto lastGlobInd = static_cast(firstGlobInd + numLocVal) - 1; + + int hypreStatus; + MatrixType* tempMatrix; + hypreStatus = HYPRE_IJMatrixCreate(Communicator::communicator(), + firstGlobInd, + lastGlobInd, + firstGlobInd, + lastGlobInd, + &tempMatrix); + VERIFY2(hypreStatus == 0, + "HYPRE_IJMatrixCreate, error code " + std::to_string(hypreStatus)); + + hypreStatus = HYPRE_IJMatrixSetObjectType(tempMatrix, HYPRE_PARCSR); + VERIFY2(hypreStatus == 0, + "HYPRE_IJMatrixSetObjectType, error code " + std::to_string(hypreStatus)); + + hypreStatus = HYPRE_IJMatrixSetRowSizes(tempMatrix, (const int*)numValsPerRow); + VERIFY2(hypreStatus == 0, + "HYPRE_IJMatrixSetRowSizes, error code " + std::to_string(hypreStatus)); + + hypreStatus = HYPRE_IJMatrixInitialize(tempMatrix); + VERIFY2(hypreStatus == 0, + "HYPRE_IJMatrixInitialize, error code " + std::to_string(hypreStatus)); + + return std::shared_ptr(tempMatrix, HYPRE_IJMatrixDestroy); +} + +//------------------------------------------------------------------------------ +// Create Hypre matrix that should self-destruct when it goes out of scope. +//------------------------------------------------------------------------------ +void +HypreFunctions:: +initializeMatrix(std::shared_ptr matrix) { + TIME_FUNCTION; + int hypreStatus; + hypreStatus = HYPRE_IJMatrixInitialize(matrix.get()); + VERIFY2(hypreStatus == 0, + "HYPRE_IJMatrixInitialize, error code " + std::to_string(hypreStatus)); +} + +//------------------------------------------------------------------------------ +// Put values into the matrix +//------------------------------------------------------------------------------ +void +HypreFunctions:: +setMatrixValues(const unsigned numRows, + const unsigned* numColsPerRow, + const unsigned* globRowInd, + const unsigned* globColInd, + const double* colVal, + std::shared_ptr opt, + std::shared_ptr matrix) { + TIME_FUNCTION; + if (numRows == 0) { + return; + } + BEGIN_CONTRACT_SCOPE + { + auto index = 0u; + for (auto i = 0u; i < numRows; ++i) { + for (auto j = 0u; j < numColsPerRow[i]; ++j, ++index) { + CHECK2(std::isfinite(colVal[index]), + "NaN found in input to Hypre matrix at index " + std::to_string(globRowInd[i]) + ", " + std::to_string(globColInd[index])); + } + } + } + END_CONTRACT_SCOPE + + int hypreStatus; + if (opt->addToValues) { + hypreStatus = HYPRE_IJMatrixAddToValues(matrix.get(), + numRows, + (int*)numColsPerRow, + (const int*)globRowInd, + (const int*)globColInd, + colVal); + } + else { + hypreStatus = HYPRE_IJMatrixSetValues(matrix.get(), + numRows, + (int*)numColsPerRow, + (const int*)globRowInd, + (const int*)globColInd, + colVal); + } + VERIFY2(hypreStatus == 0, + "HYPRE_IJMatrixSetValues, error code " + std::to_string(hypreStatus)); +} + +//------------------------------------------------------------------------------ +// Put values into hypre vector +//------------------------------------------------------------------------------ +void +HypreFunctions:: +setVectorValues(const unsigned numVals, + const unsigned firstGlobInd, + const double* val, + std::shared_ptr vec) { + TIME_FUNCTION; + if (numVals == 0) { + return; + } + BEGIN_CONTRACT_SCOPE + { + for (auto i = 0u; i < numVals; ++i) { + CHECK2(std::isfinite(val[i]), "NaN found in input to Hypre vector at index " + std::to_string(firstGlobInd + i)); + } + } + END_CONTRACT_SCOPE + + std::vector globalIndices(numVals); + std::iota(globalIndices.begin(), globalIndices.end(), firstGlobInd); + + int hypreStatus; + hypreStatus = HYPRE_IJVectorSetValues(vec.get(), + numVals, + &globalIndices[0], + val); + VERIFY2(hypreStatus == 0, + "HYPRE_IJVectorSetValues, error code " + std::to_string(hypreStatus)); + +} + +//------------------------------------------------------------------------------ +// Get values from hypre vector +//------------------------------------------------------------------------------ +void +HypreFunctions:: +getVectorValues(const unsigned numVals, + const unsigned firstGlobInd, + double* val, + std::shared_ptr vec) { + TIME_FUNCTION; + if (numVals == 0) { + return; + } + + std::vector globalIndices(numVals); + std::iota(globalIndices.begin(), globalIndices.end(), firstGlobInd); + + int hypreStatus; + hypreStatus = HYPRE_IJVectorGetValues(vec.get(), + numVals, + &globalIndices[0], + val); + VERIFY2(hypreStatus == 0, + "HYPRE_IJVectorGetValues, error code " + std::to_string(hypreStatus)); + + BEGIN_CONTRACT_SCOPE + { + for (auto i = 0u; i < numVals; ++i) { + CHECK2(std::isfinite(val[i]), "NaN found in Hypre output at index " + std::to_string(firstGlobInd + i)); + } + } + END_CONTRACT_SCOPE +} + +void +HypreFunctions:: +assembleMatrix(std::shared_ptr matrix) { + TIME_FUNCTION; + int hypreStatus; + hypreStatus = HYPRE_IJMatrixAssemble(matrix.get()); + VERIFY2(hypreStatus == 0, + "HYPRE_IJMatrixAssemble, error code " + std::to_string(hypreStatus)); +} + +//------------------------------------------------------------------------------ +// Create Hypre solver that should self-destruct when it goes out of scope. +// Right now, the constants and type (GMRES) are hardcoded. At the least, the +// kDim should be changed to a variable. +//------------------------------------------------------------------------------ +std::shared_ptr +HypreFunctions:: +createSolver(std::shared_ptr opt) { + TIME_FUNCTION; + int hypreStatus; + SolverType* tempSolver; + hypreStatus = HYPRE_ParCSRGMRESCreate(Communicator::communicator(), &tempSolver); + VERIFY2(hypreStatus == 0, + "HYPRE_ParCSRGMRESCreate, error code " + std::to_string(hypreStatus)); + + hypreStatus = HYPRE_ParCSRGMRESSetKDim(tempSolver, opt->kDim); + VERIFY2(hypreStatus == 0, + "HYPRE_PARCSRGMRESSetKDim, error code " + std::to_string(hypreStatus)); + + hypreStatus = HYPRE_ParCSRGMRESSetMinIter(tempSolver, opt->minIters); + VERIFY2(hypreStatus == 0, + "HYPRE_ParCSRGMRESSetMinIter, error code " + std::to_string(hypreStatus)); + + hypreStatus = HYPRE_ParCSRGMRESSetLogging(tempSolver, opt->logLevel); + VERIFY2(hypreStatus == 0, + "HYPRE_ParCSRGMRESSetLogging, error code " + std::to_string(hypreStatus)); + + hypreStatus = HYPRE_ParCSRGMRESSetPrintLevel(tempSolver, opt->printLevel); + VERIFY2(hypreStatus == 0, + "HYPRE_ParCSRMGRESSetPrintLevel, error code " + std::to_string(hypreStatus)); + + return std::shared_ptr(tempSolver, HYPRE_ParCSRGMRESDestroy); +} + +//------------------------------------------------------------------------------ +// Create Hypre preconditioner that should self-destruct when it goes out of +// scope. +// The constants and type (AMG) are hardcoded to start out with. +//------------------------------------------------------------------------------ +std::shared_ptr +HypreFunctions:: +createPreconditioner(std::shared_ptr opt, + std::shared_ptr solver) { + TIME_FUNCTION; + CHECK(solver); + + int hypreStatus; + SolverType* tempPreconditioner; + + switch (opt->preconditionerType) { + case HypreOptions::HyprePreconditionerType::NoPreconditioner: + return std::shared_ptr(); + case HypreOptions::HyprePreconditionerType::AMGPreconditioner: + hypreStatus = HYPRE_BoomerAMGCreate(&tempPreconditioner); + VERIFY2(hypreStatus == 0, + "HYPRE_BoomerAMGCreate, error code " + std::to_string(hypreStatus)); + + hypreStatus = HYPRE_BoomerAMGSetCoarsenType(tempPreconditioner, opt->coarsenTypeAMG); + VERIFY2(hypreStatus == 0, + "HYPRE_BoomerAMGSetCoarsenType, error code " + std::to_string(hypreStatus)); + + hypreStatus = HYPRE_BoomerAMGSetMeasureType(tempPreconditioner, opt->measure_type); + VERIFY2(hypreStatus == 0, + "HYPRE_BoomerAMGSetMeasureType, error code " + std::to_string(hypreStatus)); + + hypreStatus = HYPRE_BoomerAMGSetTol(tempPreconditioner, opt->pcTol); + VERIFY2(hypreStatus == 0, + "HYPRE_BoomerAMGSetTol, error code " + std::to_string(hypreStatus)); + + hypreStatus + = HYPRE_BoomerAMGSetStrongThreshold(tempPreconditioner, opt->strongThresholdAMG); + VERIFY2(hypreStatus == 0, + "HYPRE_BoomerAMGSetStrongThreshold, error code " + std::to_string(hypreStatus)); + + hypreStatus = HYPRE_BoomerAMGSetMaxRowSum(tempPreconditioner, opt->maxRowSumAMG); + VERIFY2(hypreStatus == 0, + "HYPRE_BoomerAMGSetMaxRowSum, error code " + std::to_string(hypreStatus)); + + if (opt->interpTypeAMG >= 0) { + hypreStatus = HYPRE_BoomerAMGSetInterpType(tempPreconditioner, opt->interpTypeAMG); + VERIFY2(hypreStatus == 0, + "HYPRE_BoomerAMGSetInterpType, error code " + std::to_string(hypreStatus)); + } + + if (opt->aggNumLevelsAMG >= 0) { + hypreStatus + = HYPRE_BoomerAMGSetAggNumLevels(tempPreconditioner, opt->aggNumLevelsAMG); + VERIFY2(hypreStatus == 0, + "HYPRE_BoomerAMGSetAggNumLevels, error code " + std::to_string(hypreStatus)); + } + if (opt->aggInterpTypeAMG >= 0) { + hypreStatus + = HYPRE_BoomerAMGSetAggInterpType(tempPreconditioner, opt->aggInterpTypeAMG); + VERIFY2(hypreStatus == 0, + "HYPRE_BoomerAMGSetAggInterpType, error code " + std::to_string(hypreStatus)); + } + if (opt->pMaxElmtsAMG >= 0){ + hypreStatus = HYPRE_BoomerAMGSetPMaxElmts(tempPreconditioner, opt->pMaxElmtsAMG); + VERIFY2(hypreStatus == 0, + "HYPRE_BoomerAMGSetPMaxElmts, error code " + std::to_string(hypreStatus)); + } + + hypreStatus = HYPRE_BoomerAMGSetTruncFactor(tempPreconditioner, opt->truncFactorAMG); + VERIFY2(hypreStatus == 0, + "HYPRE_BoomerAMGSetTruncFactor, error code " + std::to_string(hypreStatus)); + + if (opt->printLevelAMG == -1) { + opt->printLevelAMG = opt->printLevel; + } + hypreStatus = HYPRE_BoomerAMGSetPrintLevel(tempPreconditioner, opt->printLevelAMG); + VERIFY2(hypreStatus == 0, + "HYPRE_BoomerAMGSetPrintLevel, error code " + std::to_string(hypreStatus)); + + hypreStatus = HYPRE_BoomerAMGSetLogging(tempPreconditioner, opt->logLevelAMG); + VERIFY2(hypreStatus == 0, + "HYPRE_BoomerAMGSetLogging, error code " + std::to_string(hypreStatus)); + + hypreStatus = HYPRE_BoomerAMGSetMinIter(tempPreconditioner, opt->minItersAMG); + VERIFY2(hypreStatus == 0, + "HYPRE_BoomerAMGSetMinIter, error code " + std::to_string(hypreStatus)); + + hypreStatus = HYPRE_BoomerAMGSetMaxIter(tempPreconditioner, opt->maxItersAMG); + VERIFY2(hypreStatus == 0, + "HYPRE_BoomerAMGSetMaxIter, error code " + std::to_string(hypreStatus)); + + hypreStatus = HYPRE_BoomerAMGSetCycleType(tempPreconditioner, opt->cycleType); + VERIFY2(hypreStatus == 0, + "HYPRE_BoomerAMGSetCycleType, error code " + std::to_string(hypreStatus)); + + hypreStatus = HYPRE_BoomerAMGSetRelaxWt(tempPreconditioner, opt->relaxWeightAMG); + VERIFY2(hypreStatus == 0, + "HYPRE_BoomerAMGSetRelaxWt, error code " + std::to_string(hypreStatus)); + + if (opt->relaxTypeCoarseAMG == -1) { + opt->relaxTypeCoarseAMG = opt->relaxTypeAMG; + } + + hypreStatus + = HYPRE_BoomerAMGSetCycleRelaxType(tempPreconditioner, + opt->relaxTypeAMG, 1); + VERIFY2(hypreStatus == 0, + "HYPRE_BoomerAMGSetCycleRelaxType, error code " + std::to_string(hypreStatus)); + + hypreStatus + = HYPRE_BoomerAMGSetCycleRelaxType(tempPreconditioner, + opt->relaxTypeAMG, 2); + VERIFY2(hypreStatus == 0, + "HYPRE_BoomerAMGSetCycleRelaxType, error code " + std::to_string(hypreStatus)); + + hypreStatus + = HYPRE_BoomerAMGSetCycleRelaxType(tempPreconditioner, + opt->relaxTypeCoarseAMG, 3); + VERIFY2(hypreStatus == 0, + "HYPRE_BoomerAMGSetCycleRelaxType, error code " + std::to_string(hypreStatus)); + + hypreStatus + = HYPRE_BoomerAMGSetCycleNumSweeps(tempPreconditioner, + opt->cycleNumSweepsAMG, 1); + VERIFY2(hypreStatus == 0, + "HYPRE_BoomerAMGSetCycleNumSweeps, error code " + std::to_string(hypreStatus)); + + hypreStatus + = HYPRE_BoomerAMGSetCycleNumSweeps(tempPreconditioner, + opt->cycleNumSweepsAMG, 2); + VERIFY2(hypreStatus == 0, + "HYPRE_BoomerAMGSetCycleNumSweeps, error code " + std::to_string(hypreStatus)); + + if (opt->relaxTypeCoarseAMG == 19 + || opt->relaxTypeCoarseAMG == 29 + || opt->relaxTypeCoarseAMG == 9) { + opt->cycleNumSweepsCoarseAMG = 1; + } + else if (opt->cycleNumSweepsCoarseAMG == -1) { + opt->cycleNumSweepsCoarseAMG = opt->cycleNumSweepsAMG; + } + hypreStatus = HYPRE_BoomerAMGSetCycleNumSweeps(tempPreconditioner, + opt->cycleNumSweepsCoarseAMG, 3); + VERIFY2(hypreStatus == 0, + "HYPRE_BoomerAMGSetCycleNumSweeps, error code " + std::to_string(hypreStatus)); + + hypreStatus = HYPRE_BoomerAMGSetMaxLevels(tempPreconditioner, opt->maxLevelsAMG); + VERIFY2(hypreStatus == 0, + "HYPRE_BoomerAMGSetMaxLevels, error code " + std::to_string(hypreStatus)); + + hypreStatus = HYPRE_ParCSRGMRESSetPrecond(solver.get(), + HYPRE_BoomerAMGSolve, + HYPRE_BoomerAMGSetup, + tempPreconditioner); + VERIFY2(hypreStatus == 0, + "HYPRE_ParCSRGMRESSetPrecond, error code " + std::to_string(hypreStatus)); + + // Set up preconditioner + hypreStatus = HYPRE_BoomerAMGSetSetupType(tempPreconditioner, 1); + VERIFY2(hypreStatus == 0, + "HYPRE_BoomerAMGSetSetupType, error code " + std::to_string(hypreStatus)); + + return std::shared_ptr(tempPreconditioner, HYPRE_BoomerAMGDestroy); + case HypreOptions::HyprePreconditionerType::ILUPreconditioner: + hypreStatus = HYPRE_EuclidCreate(Communicator::communicator(), &tempPreconditioner); + VERIFY2(hypreStatus == 0, + "HYPRE_BoomerAMGDestroy, error code " + std::to_string(hypreStatus)); + + hypreStatus = HYPRE_EuclidSetLevel(tempPreconditioner, opt->factorLevelILU); + VERIFY2(hypreStatus == 0, + "HYPRE_EuclidSetLevel, error code " + std::to_string(hypreStatus)); + + if (opt->printLevelILU == -1) { + opt->printLevelILU = opt->printLevel; + } + hypreStatus = HYPRE_EuclidSetStats(tempPreconditioner, opt->printLevelILU); + VERIFY2(hypreStatus == 0, + "HYPRE_EuclidSetStats, error code " + std::to_string(hypreStatus)); + + if (opt->useILUT) { + hypreStatus = HYPRE_EuclidSetILUT(tempPreconditioner, opt->dropToleranceILU); + } + else { + hypreStatus = HYPRE_EuclidSetSparseA(tempPreconditioner, opt->dropToleranceILU); + VERIFY2(hypreStatus == 0, + "HYPRE_EuclidSetILUT, error code " + std::to_string(hypreStatus)); + } + + hypreStatus = HYPRE_EuclidSetRowScale(tempPreconditioner, opt->rowScaleILU); + VERIFY2(hypreStatus == 0, + "HYPRE_EuclidSetRowScale, error code " + std::to_string(hypreStatus)); + + hypreStatus = HYPRE_ParCSRGMRESSetPrecond(solver.get(), + HYPRE_EuclidSolve, + HYPRE_EuclidSetup, + tempPreconditioner); + VERIFY2(hypreStatus == 0, + "HYPRE_ParCSRGMRESSetPrecond, error code " + std::to_string(hypreStatus)); + + return std::shared_ptr(tempPreconditioner, HYPRE_EuclidDestroy); + default: + VERIFY2(false, "incorrect preconditioner type"); + return std::shared_ptr(); // So compiler doesn't complain + } +} + +void +HypreFunctions:: +setupSolver(std::shared_ptr opt, + std::shared_ptr solver, + std::shared_ptr matrix, + std::shared_ptr lhs, + std::shared_ptr rhs) { + TIME_FUNCTION; + // Get storage in CRS format + HYPRE_ParCSRMatrix hypreParMatrix; + HYPRE_ParVector hypreParLHS; + HYPRE_ParVector hypreParRHS; + + int hypreStatus; + hypreStatus = HYPRE_IJMatrixGetObject(matrix.get(), + (void **)&hypreParMatrix); + VERIFY2(hypreStatus == 0, + "HYPRE_IJMatrixGetObject, error code " + std::to_string(hypreStatus)); + + hypreStatus = HYPRE_IJVectorGetObject(lhs.get(), + (void **)&hypreParLHS); + VERIFY2(hypreStatus == 0, + "HYPRE_IJVectorGetObject, error code " + std::to_string(hypreStatus)); + + hypreStatus = HYPRE_IJVectorGetObject(rhs.get(), + (void **)&hypreParRHS); + VERIFY2(hypreStatus == 0, + "HYPRE_IJVectorGetObject, error code " + std::to_string(hypreStatus)); + + hypreStatus = HYPRE_ParCSRGMRESSetLogging(solver.get(), opt->logLevel); + VERIFY2(hypreStatus == 0, + "HYPRE_IJMatrixPrint, error code " + std::to_string(hypreStatus)); + + hypreStatus = HYPRE_ParCSRGMRESSetPrintLevel(solver.get(), opt->printLevel); + VERIFY2(hypreStatus == 0, + "HYPRE_ParCSRGMRESSetPrintLevel, error code " + std::to_string(hypreStatus)); + + hypreStatus + = HYPRE_ParCSRGMRESSetMaxIter(solver.get(), opt->maxNumberOfIterations); + VERIFY2(hypreStatus == 0, + "HYPRE_ParCSRGMRESSetMaxIter, error code " + std::to_string(hypreStatus)); + + hypreStatus = HYPRE_ParCSRGMRESSetTol(solver.get(), opt->toleranceL2); + VERIFY2(hypreStatus == 0, + "HYPRE_ParCSRGMRESSetTol, error code " + std::to_string(hypreStatus)); + + if (opt->useRobustTolerance) { + hypreStatus = HYPRE_GMRESSetSkipRealResidualCheck(solver.get(), 1); + VERIFY2(hypreStatus == 0, + "HYPRE_GMRESSetSkipRealResidualCheck, error code " + std::to_string(hypreStatus)); + } + + hypreStatus + = HYPRE_ParCSRGMRESSetAbsoluteTol(solver.get(), opt->absoluteTolerance); + VERIFY2(hypreStatus == 0, + "HYPRE_ParCSRGMRESSetAbsoluteTol, error code " + std::to_string(hypreStatus)); + + hypreStatus = HYPRE_ParCSRGMRESSetup(solver.get(), + hypreParMatrix, + hypreParRHS, + hypreParLHS); + VERIFY2(hypreStatus == 0, + "ParCSRGMRESSetup, error code " + std::to_string(hypreStatus)); +} + +//------------------------------------------------------------------------------ +// Solve the linear system +//------------------------------------------------------------------------------ +static int solveCall = 0; +std::pair +HypreFunctions:: +solveSystem(std::shared_ptr opt, + std::shared_ptr solver, + std::shared_ptr matrix, + std::shared_ptr lhs, + std::shared_ptr rhs, + std::string description) { + TIME_FUNCTION; + // Get storage in CRS format + HYPRE_ParCSRMatrix hypreParMatrix; + HYPRE_ParVector hypreParLHS; + HYPRE_ParVector hypreParRHS; + + int hypreStatus; + hypreStatus = HYPRE_IJMatrixGetObject(matrix.get(), + (void **)&hypreParMatrix); + VERIFY2(hypreStatus == 0, + "HYPRE_IJMatrixGetObject, error code " + std::to_string(hypreStatus)); + + hypreStatus = HYPRE_IJVectorGetObject(lhs.get(), + (void **)&hypreParLHS); + VERIFY2(hypreStatus == 0, + "HYPRE_IJVectorGetObject, error code " + std::to_string(hypreStatus)); + + hypreStatus = HYPRE_IJVectorGetObject(rhs.get(), + (void **)&hypreParRHS); + VERIFY2(hypreStatus == 0, + "HYPRE_IJVectorGetObject, error code " + std::to_string(hypreStatus)); + + if (opt->saveLinearSystem) { + const auto matName = "HypreMatrix_" + description + "_" + std::to_string(solveCall); + const auto rhsName = "HypreRHS_" + description + "_" + std::to_string(solveCall); + const auto initName = "HypreInit_" + description + "_" + std::to_string(solveCall); + HYPRE_IJMatrixPrint(matrix.get(), matName.c_str()); + HYPRE_IJVectorPrint(rhs.get(), rhsName.c_str()); + HYPRE_IJVectorPrint(lhs.get(), initName.c_str()); + } + + int solveStatus = HYPRE_ParCSRGMRESSolve(solver.get(), + hypreParMatrix, + hypreParRHS, + hypreParLHS); + int numIterations = 0; + if (solveStatus) { + VERIFY2(!HYPRE_CheckError(solveStatus, HYPRE_ERROR_GENERIC), + "INF or NaN in Matrix, RHS, or initial guess (" + description + ")"); + VERIFY2(!HYPRE_CheckError(solveStatus, HYPRE_ERROR_MEMORY), + "HYPRE ParCSRGMRESSetupcannot allocate memory (" + description + ")"); + VERIFY2(HYPRE_CheckError(solveStatus, HYPRE_ERROR_CONV), + "GMRES failed with (unknown) solver status (" + description + "):" << solveStatus); + HYPRE_ClearError(HYPRE_ERROR_CONV); + numIterations = opt->maxNumberOfIterations; + } + else { + hypreStatus + = HYPRE_ParCSRGMRESGetNumIterations(solver.get(), &numIterations); + VERIFY2(hypreStatus == 0, + "HYPRE_ParCSRGMRESGetNumIterations, error code " + std::to_string(hypreStatus)); + } + + double finalResNorm; + hypreStatus = HYPRE_ParCSRGMRESGetFinalRelativeResidualNorm(solver.get(), + &finalResNorm); + VERIFY2(hypreStatus == 0, + "HYPRE_ParCSRGMRESGetFinalRelativeResidualNorm, error code " + std::to_string(hypreStatus)); + + if (opt->saveLinearSystem) { + const auto lhsName = "HypreLHS_" + description + "_" + std::to_string(solveCall); + HYPRE_IJVectorPrint(lhs.get(), lhsName.c_str()); + ++solveCall; + } + + if (numIterations >= opt->maxNumberOfIterations) { + std::ostringstream oss; + oss << std::scientific << std::setprecision(2); + oss << "Hypre solver (" << description << ") did not converge, norm: " << finalResNorm; + const auto message = oss.str(); + if (opt->quitIfDiverged) { + VERIFY2(false, message); + } + else { + if (opt->warnIfDiverged && Process::getRank() == 0) { + std::cerr << message << std::endl; + } + } + } + + return std::make_pair(numIterations, finalResNorm); +} + +//------------------------------------------------------------------------------ +// Multiply the linear system +//------------------------------------------------------------------------------ +void +HypreFunctions:: +multiplySystem(std::shared_ptr opt, + std::shared_ptr matrix, + std::shared_ptr lhs, + std::shared_ptr rhs, + std::string description) { + TIME_FUNCTION; + // Get storage in CRS format + HYPRE_ParCSRMatrix hypreParMatrix; + HYPRE_ParVector hypreParLHS; + HYPRE_ParVector hypreParRHS; + + int hypreStatus; + hypreStatus = HYPRE_IJMatrixGetObject(matrix.get(), + (void **)&hypreParMatrix); + VERIFY2(hypreStatus == 0, + "HYPRE_IJMatrixGetObject, error code " + std::to_string(hypreStatus)); + + hypreStatus = HYPRE_IJVectorGetObject(lhs.get(), + (void **)&hypreParLHS); + VERIFY2(hypreStatus == 0, + "HYPRE_IJVectorGetObject, error code " + std::to_string(hypreStatus)); + + hypreStatus = HYPRE_IJVectorGetObject(rhs.get(), + (void **)&hypreParRHS); + VERIFY2(hypreStatus == 0, + "HYPRE_IJVectorGetObject, error code " + std::to_string(hypreStatus)); + + if (opt->saveLinearSystem) { + const auto matName = "HypreMatrix_" + description + "_" + std::to_string(solveCall); + const auto initName = "HypreLHS_" + description + "_" + std::to_string(solveCall); + HYPRE_IJMatrixPrint(matrix.get(), matName.c_str()); + HYPRE_IJVectorPrint(lhs.get(), initName.c_str()); + } + + // Perform matrix multiplication, y = \alpha * A x + \beta y + hypreStatus = HYPRE_ParCSRMatrixMatvec(1.0, // \alpha + hypreParMatrix, // A + hypreParLHS, // x + 0.0, // \beta + hypreParRHS); // y + + if (opt->saveLinearSystem) { + const auto rhsName = "HypreRHS_" + description + "_" + std::to_string(solveCall); + HYPRE_IJVectorPrint(rhs.get(), rhsName.c_str()); + ++solveCall; + } + + VERIFY2(hypreStatus == 0, + "HYPRE_ParCSRMatrixMatvec, error code " + std::to_string(hypreStatus)); +} + +} // end namespace Spheral diff --git a/src/Solvers/HypreFunctions.hh b/src/Solvers/HypreFunctions.hh new file mode 100644 index 000000000..a02fa0590 --- /dev/null +++ b/src/Solvers/HypreFunctions.hh @@ -0,0 +1,102 @@ +//---------------------------------Spheral++----------------------------------// +// HypreFunctions +// +// Convenience functions for Hypre +//----------------------------------------------------------------------------// +#ifndef __Spheral_HypreFunctions_hh__ +#define __Spheral_HypreFunctions_hh__ + +#include +#include + +struct hypre_IJMatrix_struct; +struct hypre_IJVector_struct; +struct hypre_Solver_struct; + +namespace Spheral { + +class HypreOptions; + +// These functions do not check for safety, so make sure objects are initialized +// Here is an example of how these might be called in order: +// rhs = createVector(...); +// lhs = createVector(...); +// mat = createMatrix(...); +// setMatrixValues(...); +// assembleMatrix(...); +// sol = createSolver(...); +// prec = createPreconditioner(...); +// setupSolver(...); +// setVectorValues(...); // rhs +// setVectorValues(...); // lhs +// solveSystem(...); +// getVectorValues(...); // lhs +struct HypreFunctions { + typedef struct hypre_IJMatrix_struct MatrixType; + typedef struct hypre_IJVector_struct VectorType; + typedef struct hypre_Solver_struct SolverType; + + // Get Hypre vectors and matrices that self-destruct + static std::shared_ptr createVector(const unsigned numLocVal, + const unsigned firstGlobInd); + static std::shared_ptr createMatrix(const unsigned numLocVal, + const unsigned firstGlobInd, + const unsigned* numValsPerRow); + + // Reinitialize matrix before setting new values + static void initializeMatrix(std::shared_ptr matrix); + + // Set values in the matrix + static void setMatrixValues(const unsigned numRows, + const unsigned* numColsPerRow, + const unsigned* globRowInd, + const unsigned* globColInd, + const double* colVal, + std::shared_ptr opt, + std::shared_ptr matrix); + + // Set or get the values of a Hypre vector + static void setVectorValues(const unsigned numVals, + const unsigned firstGlobInd, + const double* val, + std::shared_ptr hypreVector); + static void getVectorValues(const unsigned numVals, + const unsigned firstGlobInd, + double* val, + std::shared_ptr hypreVector); + + // Assemble the matrix + static void assembleMatrix(std::shared_ptr matrix); + + // Get solver and preconditioner that self-destructs + static std::shared_ptr createSolver(std::shared_ptr opt); + static std::shared_ptr createPreconditioner(std::shared_ptr opt, + std::shared_ptr solver); + + // Set remaining solver options and set up preconditioner + static void setupSolver(std::shared_ptr opt, + std::shared_ptr solver, + std::shared_ptr matrix, + std::shared_ptr lhs, + std::shared_ptr rhs); + + // Solve the linear system + // Return iterations and residual norm + static std::pair solveSystem(std::shared_ptr opt, + std::shared_ptr solver, + std::shared_ptr matrix, + std::shared_ptr lhs, + std::shared_ptr rhs, + std::string description); + + // Multiply the linear system + static void multiplySystem(std::shared_ptr opt, + std::shared_ptr matrix, + std::shared_ptr lhs, + std::shared_ptr rhs, + std::string description); +}; // end struct HypreFunctions + +} // end namespace Spheral + +#endif diff --git a/src/Solvers/HypreLinearSolver.cc b/src/Solvers/HypreLinearSolver.cc new file mode 100644 index 000000000..7bacbae84 --- /dev/null +++ b/src/Solvers/HypreLinearSolver.cc @@ -0,0 +1,200 @@ +//---------------------------------Spheral++----------------------------------// +// HypreLinearSolver +// +// Represents a solver for a linear system of equations +//----------------------------------------------------------------------------// + +#include "HypreLinearSolver.hh" + +#include "Solvers/HypreFunctions.hh" +#include "Solvers/HypreOptions.hh" +#include "Solvers/IncrementalStatistic.hh" +#include "Utilities/DBC.hh" +#include "Utilities/Timer.hh" + +namespace Spheral { + +//------------------------------------------------------------------------------ +// Constructor +//------------------------------------------------------------------------------ +HypreLinearSolver:: +HypreLinearSolver(std::shared_ptr options) : + mInitialized(false), + mFilling(false), + mAssembled(false), + mFinalized(false), + mHypreOptions(options), + mIterationStatistics( + std::make_shared>(options->meanIterationsGuess, + "HypreLinearSolver iterations", + options->printIterations)), // print + mFinalResidualStatistics( + std::make_shared>(options->toleranceL2, + "HypreLinearSolver residual", + false)) { // print + this->setDescription("HypreLinearSolver"); +} + +//------------------------------------------------------------------------------ +// Initialize the matrices and vectors +//------------------------------------------------------------------------------ +void +HypreLinearSolver:: +initialize(const unsigned numLocVal, + const unsigned firstGlobInd, + const unsigned* numValsPerRow) { + TIME_SCOPE("HypreLinearSolver::initialize"); + mInitialized = false; + mAssembled = false; + mFinalized = false; + + mHypreLHS = HypreFunctions::createVector(numLocVal, firstGlobInd); + mHypreRHS = HypreFunctions::createVector(numLocVal, firstGlobInd); + mHypreMatrix = HypreFunctions::createMatrix(numLocVal, firstGlobInd, numValsPerRow); + mHypreSolver = HypreFunctions::createSolver(mHypreOptions); + mHyprePreconditioner = HypreFunctions::createPreconditioner(mHypreOptions, mHypreSolver); + + mInitialized = true; +} + +void +HypreLinearSolver:: +beginFill() { + TIME_SCOPE("HypreLinearSolver::beginFill"); + VERIFY(mInitialized); + HypreFunctions::initializeMatrix(mHypreMatrix); + mFilling = true; +} + +//------------------------------------------------------------------------------ +// Set matrix values +//------------------------------------------------------------------------------ +void +HypreLinearSolver:: +setMatRow(const unsigned numVals, + const unsigned globRowInd, + const unsigned* globColInd, + const double* colVal) { + TIME_SCOPE("HypreLinearSolver::setMatRow"); + VERIFY(mFilling); + HypreFunctions::setMatrixValues(1, &numVals, &globRowInd, + globColInd, colVal, + mHypreOptions, mHypreMatrix); +} + +void +HypreLinearSolver:: +setMatRows(const unsigned numRows, + const unsigned* numColsPerRow, + const unsigned* globRowInd, + const unsigned* globColInd, + const double* colVal) { + TIME_SCOPE("HypreLinearSolver::setMatRows"); + VERIFY(mFilling); + HypreFunctions::setMatrixValues(numRows, numColsPerRow, globRowInd, + globColInd, colVal, + mHypreOptions, mHypreMatrix); +} + +//------------------------------------------------------------------------------ +// Assemble matrix after fill +//------------------------------------------------------------------------------ +void +HypreLinearSolver:: +assemble() { + TIME_SCOPE("HypreLinearSolver::assemble"); + VERIFY(mFilling); + HypreFunctions::assembleMatrix(mHypreMatrix); + mFilling = false; + mAssembled = true; +} + +//------------------------------------------------------------------------------ +// Create solver and preconditioner +//------------------------------------------------------------------------------ +void +HypreLinearSolver:: +finalize() { + TIME_SCOPE("HypreLinearSolver::finalize"); + VERIFY(mAssembled); + HypreFunctions::setupSolver(mHypreOptions, mHypreSolver, mHypreMatrix, mHypreLHS, mHypreRHS); + mFinalized = true; +} + +//------------------------------------------------------------------------------ +// Set/get values in LHS and RHS +//------------------------------------------------------------------------------ +void +HypreLinearSolver:: +set(const Component c, + const unsigned numVals, + const unsigned firstGlobInd, + const double* val) { + TIME_SCOPE("HypreLinearSolver::set"); + VERIFY(mInitialized); + switch (c) { + case RHS: + HypreFunctions::setVectorValues(numVals, firstGlobInd, val, mHypreRHS); + break; + case LHS: + HypreFunctions::setVectorValues(numVals, firstGlobInd, val, mHypreLHS); + break; + } +} + +void +HypreLinearSolver:: +get(const Component c, + const unsigned numVals, + const unsigned firstGlobInd, + double* val) { + TIME_SCOPE("HypreLinearSolver::get"); + VERIFY(mInitialized); + switch (c) { + case RHS: + HypreFunctions::getVectorValues(numVals, firstGlobInd, val, mHypreRHS); + break; + case LHS: + HypreFunctions::getVectorValues(numVals, firstGlobInd, val, mHypreLHS); + break; + } +} + +//------------------------------------------------------------------------------ +// Solve the system of equations in place +//------------------------------------------------------------------------------ +void +HypreLinearSolver:: +solve() { + TIME_SCOPE("HypreLinearSolver::solve"); + VERIFY(mFinalized); + auto [it, norm] = HypreFunctions::solveSystem(mHypreOptions, mHypreSolver, mHypreMatrix, mHypreLHS, mHypreRHS, this->getDescription()); + mIterationStatistics->add(it); + mFinalResidualStatistics->add(norm); +} + +//------------------------------------------------------------------------------ +// Multiply by the matrix +//------------------------------------------------------------------------------ +void +HypreLinearSolver:: +multiply() { + TIME_SCOPE("HypreLinearSolver::multiply"); + VERIFY(mAssembled); + HypreFunctions::multiplySystem(mHypreOptions, mHypreMatrix, mHypreLHS, mHypreRHS, this->getDescription()); +} + +//------------------------------------------------------------------------------ +// Set the description, including for the stats +//------------------------------------------------------------------------------ +inline +void +HypreLinearSolver:: +setDescription(std::string desc) { + this->mDescription = desc; + mIterationStatistics->setName(desc + " iter"); + mFinalResidualStatistics->setName(desc + " norm"); +} + + +} // end namespace Spheral diff --git a/src/Solvers/HypreLinearSolver.hh b/src/Solvers/HypreLinearSolver.hh new file mode 100644 index 000000000..05583c054 --- /dev/null +++ b/src/Solvers/HypreLinearSolver.hh @@ -0,0 +1,105 @@ +//---------------------------------Spheral++----------------------------------// +// HypreLinearSolver +// +// Represents a solver for a linear system of equations +//----------------------------------------------------------------------------// +#ifndef __Spheral_HypreLinearSolver_hh__ +#define __Spheral_HypreLinearSolver_hh__ + +#include +#include + +#include "LinearSolver.hh" + +struct hypre_IJMatrix_struct; +struct hypre_IJVector_struct; +struct hypre_Solver_struct; + +namespace Spheral { + +class HypreOptions; +template class IncrementalStatistic; + +class HypreLinearSolver : public LinearSolver { +public: + typedef struct hypre_IJMatrix_struct MatrixType; + typedef struct hypre_IJVector_struct VectorType; + typedef struct hypre_Solver_struct SolverType; + + using LinearSolver::Component; + + // Constructor + HypreLinearSolver(std::shared_ptr options); + + // Create matrices and vectors + virtual void initialize(const unsigned numLocVal, + const unsigned firstGlobInd, + const unsigned* numValsPerRow) override; + + // Begin the matrix fill + virtual void beginFill() override; + + // Set values in matrix + virtual void setMatRow(const unsigned numVals, // Number of values in this row + const unsigned globRowInd, // Global index of this row + const unsigned* globColInd, // [numVals] Global indices for columns + const double* colVal) override; // [numVals] Values + virtual void setMatRows(const unsigned numRows, // Number of rows + const unsigned* numColsPerRow, // [numRows] Number of columns for each of these rows + const unsigned* globRowInd, // [numRows] Global indices for these rows + const unsigned* globColInd, // [numColsPerRow[numRows]] Global column indies for each value + const double* colVal) override; // [numColsPerRow[numRows]] Values + + // Assemble matrix after fill + virtual void assemble() override; + + // Create solver/preconditioner + virtual void finalize() override; + + // Set/get values + virtual void set(const Component c, + const unsigned numVals, // Number of vector values to set + const unsigned firstGlobInd, // First global index for values + const double* val) override; // [numVals] Values + virtual void get(const Component c, + const unsigned numVals, // Number of vector values to set + const unsigned firstGlobInd, // First global index for values + double* val) override; // [numVals] Values + + // Solve the system of equations in place + virtual void solve() override; + + // Multiply the matrix by a vector + virtual void multiply() override; + + // Set the description + virtual void setDescription(std::string desc) override; + + // Return statistics + virtual std::vector>> statistics() const override; + +private: + // Which stages have we completed? + bool mInitialized; + bool mFilling; + bool mAssembled; + bool mFinalized; + + // Data + std::shared_ptr mHypreOptions; + std::shared_ptr mHypreMatrix; + std::shared_ptr mHypreLHS; + std::shared_ptr mHypreRHS; + std::shared_ptr mHypreSolver; + std::shared_ptr mHyprePreconditioner; + + // Iteration statistics + std::shared_ptr> mIterationStatistics; + std::shared_ptr> mFinalResidualStatistics; +}; // end class HypreLinearSolver + +} // end namespace Spheral + +#include "HypreLinearSolverInline.hh" + +#endif diff --git a/src/Solvers/HypreLinearSolverInline.hh b/src/Solvers/HypreLinearSolverInline.hh new file mode 100644 index 000000000..b22a11319 --- /dev/null +++ b/src/Solvers/HypreLinearSolverInline.hh @@ -0,0 +1,14 @@ +namespace Spheral { + +//------------------------------------------------------------------------------ +// Return statistics for number of iterations and tolerance reached +//------------------------------------------------------------------------------ +inline +std::vector>> +HypreLinearSolver:: +statistics() const { + return {mIterationStatistics, mFinalResidualStatistics}; +} + +} // end namespace Spheral + diff --git a/src/Solvers/HypreOptions.hh b/src/Solvers/HypreOptions.hh new file mode 100644 index 000000000..f5aa02eab --- /dev/null +++ b/src/Solvers/HypreOptions.hh @@ -0,0 +1,83 @@ +//---------------------------------Spheral++----------------------------------// +// HypreOptions +// +// Holds the options for Hypre solvers and preconditioners +//----------------------------------------------------------------------------// +#ifndef __Spheral_HypreOptions_hh__ +#define __Spheral_HypreOptions_hh__ + +namespace Spheral { + +struct HypreOptions { + HypreOptions() { } + + // Quit program if solution did not converge + bool quitIfDiverged = false; + bool warnIfDiverged = true; + + // Sum all row values that share an index + bool sumDuplicates = true; + bool addToValues = false; + + // Shared options + int logLevel = 0; + int printLevel = 0; + + // Options used in GMRES solve + bool saveLinearSystem = false; + bool printIterations = false; + int meanIterationsGuess = 5; + int maxNumberOfIterations = 100; + double toleranceL2 = 1.e-14; + int useRobustTolerance = 1; + double absoluteTolerance = 0.0; + + // Options used in GMRES initialization + int kDim = 10; + int minIters = 1; + + // General preconditioner options + enum HyprePreconditionerType { + NoPreconditioner, + AMGPreconditioner, + ILUPreconditioner + }; + HyprePreconditionerType preconditionerType = HyprePreconditionerType::AMGPreconditioner; + + // Options used in BoomerAMG initialization + int measure_type = 0; + double pcTol = 0.; + int minItersAMG = 1; + int maxItersAMG = 1; + int cycleType = 1; + int coarsenTypeAMG = 6; + double strongThresholdAMG = 0.25; + double maxRowSumAMG = 0.5; + int interpTypeAMG = -1; + int aggNumLevelsAMG = -1; + int aggInterpTypeAMG = -1; + int pMaxElmtsAMG = -1; + double truncFactorAMG = 0.0; + int printLevelAMG = -1; + int logLevelAMG = 0; + double relaxWeightAMG = 1.0; + int relaxTypeAMG = 8; + int relaxTypeCoarseAMG = 19; + int cycleNumSweepsAMG = 1; + int cycleNumSweepsCoarseAMG = -1; + int maxLevelsAMG = 25; + + // Options used in Euclid (ILU) initialization + bool useILUT = false; + int factorLevelILU = 1; + int rowScaleILU = 0; + int printLevelILU = -1; + double dropToleranceILU = 0; + + // Postprocessing + bool zeroOutNegativities = false; +}; // end struct HypreOptions + +} // end namespace Spheral + +#endif diff --git a/src/Solvers/IncrementalStatistic.hh b/src/Solvers/IncrementalStatistic.hh new file mode 100644 index 000000000..9f85454d5 --- /dev/null +++ b/src/Solvers/IncrementalStatistic.hh @@ -0,0 +1,75 @@ +//---------------------------------Spheral++----------------------------------// +// IncrementalStatistic +// +// Allows tracking of incrementally added stats, such as iteration counts +// For now, all variables use the same DataType +// The mean guess prevents overflow for data with a large mean +//----------------------------------------------------------------------------// +#ifndef __Spheral_IncrementalStatistic_hh__ +#define __Spheral_IncrementalStatistic_hh__ + +#include + +namespace Spheral { + +template +class IncrementalStatistic { +public: + // Constructor + IncrementalStatistic(const DataType meanGuess = 0, + const std::string name = "IncrementalStatistic", + bool printAdd = false); + + // Add a point of data + void add(const DataType data); + + // Calculated data + DataType mean() const; + DataType variance() const; + DataType total() const; + + // Stored data + std::string name() const; + DataType meanGuess() const; + DataType min() const; + DataType max() const; + DataType numPoints() const; + DataType shiftedSum() const; + DataType shiftedSum2() const; + + // Print to console + void print() const; + + // Print summary to console + void printSummary() const; + + // Set name + void setName(std::string newName); + +private: + // Convenience function + void getStats(DataType& min, + DataType& max, + DataType& numPoints, + DataType& mean, + DataType& variance, + DataType& total) const; + + // Input data + const DataType mMeanGuess; + std::string mName; + bool mPrint; + + // Incrementally updated data + DataType mMin; + DataType mMax; + DataType mNumPoints; + DataType mShiftedSum; + DataType mShiftedSum2; +}; + +} // end namespace Spheral + +#include "IncrementalStatisticInline.hh" + +#endif diff --git a/src/Solvers/IncrementalStatisticInline.hh b/src/Solvers/IncrementalStatisticInline.hh new file mode 100644 index 000000000..c6c52cbae --- /dev/null +++ b/src/Solvers/IncrementalStatisticInline.hh @@ -0,0 +1,295 @@ +#include +#include +#include +#include + +#include "Distributed/Communicator.hh" +#include "Utilities/DataTypeTraits.hh" + +namespace { // anonymous + +template +inline +DataType +calculateMean(DataType numPoints, + DataType shiftedSum, + DataType meanGuess) { + return shiftedSum / numPoints + meanGuess; +} + +template +inline +DataType +calculateVariance(DataType numPoints, + DataType shiftedSum, + DataType shiftedSum2) { + return (shiftedSum2 - shiftedSum * shiftedSum / numPoints) / numPoints; +} + +template +inline +DataType +calculateTotal(DataType numPoints, + DataType shiftedSum, + DataType meanGuess) { + return shiftedSum + numPoints * meanGuess; +} + +} // namespace anonymous + +namespace Spheral { + +//------------------------------------------------------------------------------ +// Constructor +//------------------------------------------------------------------------------ +template +inline +IncrementalStatistic:: +IncrementalStatistic(const DataType meanGuess, + std::string name, + bool print): + mMeanGuess(meanGuess), + mName(name), + mPrint(print), + mMin(std::numeric_limits::max()), + mMax(std::numeric_limits::min()), + mNumPoints(0), + mShiftedSum(0), + mShiftedSum2(0) { +} + +//------------------------------------------------------------------------------ +// Add a data point +//------------------------------------------------------------------------------ +template +inline +void +IncrementalStatistic:: +add(const DataType data) { + if (data < mMin) mMin = data; + if (data > mMax) mMax = data; + mNumPoints += 1; + const DataType delta = data - mMeanGuess; + mShiftedSum += delta; + mShiftedSum2 += delta * delta; + + if (mPrint) { + std::cout << mName; + std::cout << "\tadd: " << data; + std::cout << std::endl; + } +} + +//------------------------------------------------------------------------------ +// Calculate mean +//------------------------------------------------------------------------------ +template +inline +DataType +IncrementalStatistic:: +mean() const { + return calculateMean(mNumPoints, mShiftedSum, mMeanGuess); +} + +//------------------------------------------------------------------------------ +// Calculate variance +//------------------------------------------------------------------------------ +template +inline +DataType +IncrementalStatistic:: +variance() const { + return calculateVariance(mNumPoints, mShiftedSum, mShiftedSum2); +} + +//------------------------------------------------------------------------------ +// Calculate total +//------------------------------------------------------------------------------ +template +inline +DataType +IncrementalStatistic:: +total() const { + return calculateTotal(mNumPoints, mShiftedSum, mMeanGuess); +} + + +//------------------------------------------------------------------------------ +// Set name +//------------------------------------------------------------------------------ +template +inline +void +IncrementalStatistic:: +setName(std::string name) { + mName = name; + return; +} + +//------------------------------------------------------------------------------ +// Get mean guess +//------------------------------------------------------------------------------ +template +inline +DataType +IncrementalStatistic:: +min() const { + return mMin; +} + +//------------------------------------------------------------------------------ +// Get mean guess +//------------------------------------------------------------------------------ +template +inline +DataType +IncrementalStatistic:: +max() const { + return mMax; +} + +//------------------------------------------------------------------------------ +// Get mean guess +//------------------------------------------------------------------------------ +template +inline +DataType +IncrementalStatistic:: +meanGuess() const { + return mMeanGuess; +} + +//------------------------------------------------------------------------------ +// Get numPoints +//------------------------------------------------------------------------------ +template +inline +DataType +IncrementalStatistic:: +numPoints() const { + return mNumPoints; +} + +//------------------------------------------------------------------------------ +// Get shiftedSum +//------------------------------------------------------------------------------ +template +inline +DataType +IncrementalStatistic:: +shiftedSum() const { + return mShiftedSum; +} + +//------------------------------------------------------------------------------ +// Get shiftedSum2 +//------------------------------------------------------------------------------ +template +inline +DataType +IncrementalStatistic:: +shiftedSum2() const { + return mShiftedSum2; +} + +//------------------------------------------------------------------------------ +// Print results +//------------------------------------------------------------------------------ +template +inline +void +IncrementalStatistic:: +print() const { + typedef std::pair PairType; + + auto outputPair = [](PairType val) { + std::cout << std::setw(20) << val.first; + std::cout << std::setw(14) << std::setprecision(6) << val.second; + }; + + DataType min, max, numPoints, mean, variance, total; + getStats(min, max, numPoints, mean, variance, total); + + if (Process::getRank() == 0) { + std::cout << mName << std::endl; + if (numPoints == 0) { + PairType val = {"points", numPoints}; + outputPair(val); + std::cout << "(no data recorded)" << std::endl; + } + else { + // Output in two columns + std::vector vals = + {{"points", numPoints}, {"total", total}, + {"mean", mean}, {"minimum", min}, + {"variance", variance}, {"maximum", max}}; + auto index = 0; + for (auto val : vals) { + outputPair(val); + if (index % 2 == 1) { + std::cout << std::endl; + } + index += 1; + } + } + } +} + +template +inline +void +IncrementalStatistic:: +printSummary() const { + DataType min, max, numPoints, mean, variance, total; + getStats(min, max, numPoints, mean, variance, total); + + if (Process::getRank() == 0) { + std::cout << std::left << std::setw(22) << mName; + if (numPoints == 0) { + std::cout << " (no data recorded)" << std::endl; + } + else { + std::cout << std::scientific << std::setprecision(4); + std::cout << " mean=" << mean + << " min=" << min + << " max=" << max + << " tot=" << total << std::endl; + } + } +} + +// Get raw stats communicated over all ranks +template +inline +void +IncrementalStatistic:: +getStats(DataType& min, + DataType& max, + DataType& numPoints, + DataType& mean, + DataType& variance, + DataType& total) const { + // Communicate raw data + DataType shiftedSum, shiftedSum2; + const auto mpiDataType = DataTypeTraits::MpiDataType(); + MPI_Allreduce(&mMin, &min, 1, mpiDataType, MPI_MIN, Communicator::communicator()); + MPI_Allreduce(&mMax, &max, 1, mpiDataType, MPI_MAX, Communicator::communicator()); + MPI_Allreduce(&mNumPoints, &numPoints, 1, mpiDataType, MPI_SUM, Communicator::communicator()); + MPI_Allreduce(&mShiftedSum, &shiftedSum, 1, mpiDataType, MPI_SUM, Communicator::communicator()); + MPI_Allreduce(&mShiftedSum2, &shiftedSum2, 1, mpiDataType, MPI_SUM, Communicator::communicator()); + + // Calculate mean and variance with total num points, since they are divided through + mean = calculateMean(numPoints, shiftedSum, mMeanGuess); + variance = calculateVariance(numPoints, shiftedSum, shiftedSum2); + total = calculateTotal(numPoints, shiftedSum, mMeanGuess); + + // Divide total and numPoints by procs, since these are usually added similarly on every proc + const auto numProcs = Process::getTotalNumberOfProcesses(); + numPoints /= numProcs; + total /= numProcs; +} + + + + +} // end namespace Spheral diff --git a/src/Solvers/LinearSolver.cc b/src/Solvers/LinearSolver.cc new file mode 100644 index 000000000..401603c05 --- /dev/null +++ b/src/Solvers/LinearSolver.cc @@ -0,0 +1,19 @@ +//---------------------------------Spheral++----------------------------------// +// LinearSolver +// +// Represents a solver for a linear system of equations +//----------------------------------------------------------------------------// + +#include "LinearSolver.hh" + +namespace Spheral { + +//------------------------------------------------------------------------------ +// Constructor +//------------------------------------------------------------------------------ +LinearSolver:: +LinearSolver(): + mDescription("LinearSolver") { +} + +} // end namespace Spheral diff --git a/src/Solvers/LinearSolver.hh b/src/Solvers/LinearSolver.hh new file mode 100644 index 000000000..c43fd4ba9 --- /dev/null +++ b/src/Solvers/LinearSolver.hh @@ -0,0 +1,100 @@ +//---------------------------------Spheral++----------------------------------// +// LinearSolver +// +// Represents a solver for a linear system of equations +// A given type only needs to fill in the pure virtual methods, and the rest +// is done through the setMap and setData functions. +//----------------------------------------------------------------------------// +#ifndef __Spheral_LinearSolver_hh__ +#define __Spheral_LinearSolver_hh__ + +#include +#include +#include + +namespace Spheral { + +template class IncrementalStatistic; + +/* To set up the class for a solve (each time the matrix changes) + initialize(...); + for (auto i = 0; i < numLocalRows; ++i) { + setMatRow(...); + } + assemble(...); + finalize(...); + To solve or multiply (as many times as needed after initialization) + setLHS(...); // Set guess + setRHS(...); // Set source + solve(); // Solve in place + getLHS(...); // Get solution + +*/ +class LinearSolver { +public: + // Specify which part of the system to apply set/get to + enum Component { + LHS, + RHS, + }; + + // Constructor + LinearSolver(); + + // Create matrices and vectors + virtual void initialize(const unsigned numLocVal, + const unsigned firstGlobInd, + const unsigned* numValsPerRow) = 0; + + // Begin matrix fill + virtual void beginFill() = 0; + + // Set values in matrix + virtual void setMatRow(const unsigned numVals, // Number of values in this row + const unsigned globRowInd, // Global index of this row + const unsigned* globColInd, // [numVals] Global indices for columns + const double* colVal) = 0; // [numVals] Values + virtual void setMatRows(const unsigned numRows, // Number of rows + const unsigned* numColsPerRow, // [numRows] Number of columns for each of these rows + const unsigned* globRowInd, // [numRows] Global indices for these rows + const unsigned* globColInd, // [numColsPerRow[numRows]] Global column indies for each value + const double* colVal) = 0; // [numColsPerRow[numRows]] Values + + // Assemble matrix after fill + virtual void assemble() = 0; + + // Create solver/preconditioner + virtual void finalize() = 0; + + // Set/get values in vectors + virtual void set(const Component c, // Which vector to apply this to + const unsigned numVals, // Number of vector values to set + const unsigned firstGlobInd, // First global index for values + const double* val) = 0; // [numVals] Values + virtual void get(const Component c, // Which vector to apply this to + const unsigned numVals, // Number of vector values to set + const unsigned firstGlobInd, // First global index for values + double* val) = 0; // [numVals] Values + + // Solve the system of equations in place + virtual void solve() = 0; + + // Multiply the matrix by a vector + virtual void multiply() = 0; + + // Set/get description + virtual void setDescription(std::string desc); + virtual std::string getDescription() const; + + // Return statistics if available, empty vector if not + virtual std::vector>> statistics() const; + +protected: + std::string mDescription; +}; // end class LinearSolver + +} // end namespace Spheral + +#include "LinearSolverInline.hh" + +#endif diff --git a/src/Solvers/LinearSolverInline.hh b/src/Solvers/LinearSolverInline.hh new file mode 100644 index 000000000..447a25651 --- /dev/null +++ b/src/Solvers/LinearSolverInline.hh @@ -0,0 +1,30 @@ +namespace Spheral { + +//------------------------------------------------------------------------------ +// Return empty statistics +//------------------------------------------------------------------------------ +inline +std::vector>> +LinearSolver:: +statistics() const { + return std::vector>>(); +} + +//------------------------------------------------------------------------------ +// Set/get description +//------------------------------------------------------------------------------ +inline +void +LinearSolver:: +setDescription(std::string desc) { + mDescription = desc; +} + +inline +std::string +LinearSolver:: +getDescription() const { + return mDescription; +} + +} // end namespace Spheral diff --git a/src/Solvers/containsConstantNodes.hh b/src/Solvers/containsConstantNodes.hh new file mode 100644 index 000000000..54e004d85 --- /dev/null +++ b/src/Solvers/containsConstantNodes.hh @@ -0,0 +1,25 @@ +//---------------------------------Spheral++----------------------------------// +// containsConstantNodes +// +// Checks whether boundaries contain constant nodes +//----------------------------------------------------------------------------// +#ifndef __Spheral_containsConstantNodes_hh__ +#define __Spheral_containsConstantNodes_hh__ + +#include + +namespace Spheral { + +template class Boundary; + +template +inline bool containsConstantNodes(const Boundary* boundary); + +template +inline bool containsConstantNodes(const std::vector*>& boundaries); + +} // end namespace Spheral + +#include "containsConstantNodesInline.hh" + +#endif diff --git a/src/Solvers/containsConstantNodesInline.hh b/src/Solvers/containsConstantNodesInline.hh new file mode 100644 index 000000000..b14ca0c22 --- /dev/null +++ b/src/Solvers/containsConstantNodesInline.hh @@ -0,0 +1,28 @@ +#include "Boundary/ConstantBoundary.hh" +#include "Boundary/InflowOutflowBoundary.hh" + +namespace Spheral { + +template +inline +bool +containsConstantNodes(const Boundary* boundary) { + // This is bad. It would probably be better to add a bool to the boundaries + // that is true for boundaries that contain constant nodes. + return (dynamic_cast*>(boundary) != nullptr + || dynamic_cast*>(boundary) != nullptr); +} + +template +inline +bool +containsConstantNodes(const std::vector*>& boundaries) { + for (auto& boundary : boundaries) { + if (containsConstantNodes(boundary)) { + return true; + } + } + return false; +} + +} // end namespace Spheral diff --git a/tests/integration.ats b/tests/integration.ats index cebe3c304..774223986 100644 --- a/tests/integration.ats +++ b/tests/integration.ats @@ -93,6 +93,7 @@ source("unit/Mesh/testLineMesh.py") # Solvers unit tests source("unit/Solvers/testKinsol.py") +source("unit/Solvers/testLinearSolver.py") # MPI python interface unit tests source("unit/Distributed/testMPI4PYasPYMPI.py") diff --git a/tests/unit/Solvers/testLinearSolver.py b/tests/unit/Solvers/testLinearSolver.py new file mode 100644 index 000000000..de3d8254c --- /dev/null +++ b/tests/unit/Solvers/testLinearSolver.py @@ -0,0 +1,406 @@ +#ATS:test(SELF, "--dimension 1 --nx 64 --preconditioner 'none' --numpyTest True", np=1, label="Linear solver, 1d GMRES") +#ATS:test(SELF, "--dimension 2 --nx 16 --preconditioner 'ilu' --reflecting True --numpyTest True", np=4, label="Linear solver, 2d ILU") +#ATS:test(SELF, "--dimension 3 --nx 8 --preconditioner 'amg' --periodic True --numpyTest True --nPerh 4.01", np=8, label="Linear solver, 3d AMG") +# #ATS:test(SELF, "--dimension 3 --nx 50 --preconditioner 'amg' --reflecting True --nPerh 4.01 --numTests 10 --numSolves 5 --caliperConfig 'spot,runtime-report'", np=32, label="Linear solver, timing check") + +#------------------------------------------------------------------------------- +# Set up a system of equations and solve them with the linear solver operator +#------------------------------------------------------------------------------- +from Spheral import * +from SpheralTestUtilities import * +import numpy as np +import time, mpi +title("Linear solver test") + +#------------------------------------------------------------------------------- +# Generic problem parameters +#------------------------------------------------------------------------------- +commandLine( + # Spatial stuff + dimension = 1, + nx = 32, + x0 = 0.0, + x1 = 1.0, + + # Interpolation kernel + nPerh = 2.01, + + # Solver options + preconditioner = "amg", + + # Test reinitialization + numTests = 4, # Multiple tests, same graph but different data + numSolves = 2, # Solve the same thing multiple times for timing + + # Include boundary conditions + reflecting = False, + periodic = False, + + # Test options + tolerance = 1e-4, + numpyTest = False, + seed = None, + saveSystem = False, # Save hypre matrices + exactGuess = False, +) +exec("from Spheral%id import *" % dimension, globals()) + +if seed is None: + seed = np.random.randint(0, 2**8) +output("seed") +np.random.seed(seed) + +#------------------------------------------------------------------------------- +# Create nodes +#------------------------------------------------------------------------------ +units = MKS() +output("units") + +eos = GammaLawGas(5.0/3.0, 1.0, units) +output("eos") + +WT = TableKernel(WendlandC4Kernel(), 1000) +kernelExtent = WT.kernelExtent +output("WT") + +delta = (x1 - x0) / nx +hmid = delta * nPerh * kernelExtent +hmin = hmid * 1.e-3 +hmax = hmid * 1.e3 +nodes = makeFluidNodeList("nodes", eos, + hmin = hmin, + hmax = hmax, + nPerh = nPerh, + kernelExtent = kernelExtent) +output("nodes") +output("nodes.hmin") +output("nodes.hmax") +output("nodes.nodesPerSmoothingScale") + +#------------------------------------------------------------------------------- +# Seed the nodes +#------------------------------------------------------------------------------- +rho0 = 1.0 +if dimension == 1: + from DistributeNodes import distributeNodesInRange1d + distributeNodesInRange1d([(nodes, nx, rho0, (x0, x1))], + nPerh = nPerh) +elif dimension == 2: + from GenerateNodeDistribution2d import * + generator = GenerateNodeDistribution2d(distributionType="lattice", + nRadial = nx, nTheta = nx, + xmin = (x0, x0), + xmax = (x1, x1), + rho = rho0, + nNodePerh = nPerh) + if mpi.procs > 1: + from VoronoiDistributeNodes import distributeNodes2d + else: + from DistributeNodes import distributeNodes2d + distributeNodes2d((nodes, generator)) +else: + from GenerateNodeDistribution3d import * + generator = GenerateNodeDistribution3d(distributionType="lattice", + n1 = nx, n2 = nx, n3 = nx, + xmin = (x0, x0, x0), + xmax = (x1, x1, x1), + rho=rho0, + nNodePerh = nPerh) + if mpi.procs > 1: + from VoronoiDistributeNodes import distributeNodes3d + else: + from DistributeNodes import distributeNodes3d + distributeNodes3d((nodes, generator)) +output("nodes.numNodes") +numLocal = nodes.numInternalNodes +output("numLocal") + +#------------------------------------------------------------------------------- +# Create DataBase +#------------------------------------------------------------------------------- +dataBase = DataBase() +dataBase.appendNodeList(nodes) +output("dataBase") + +#------------------------------------------------------------------------------- +# Set initial conditions +#------------------------------------------------------------------------------ +inp = dataBase.newFluidScalarFieldList(0.0, "input") +out = dataBase.newFluidScalarFieldList(0.0, "output") +inp0 = dataBase.newFluidScalarFieldList(0.0, "initial input") +out0 = dataBase.newFluidScalarFieldList(0.0, "initial output") +for i in range(numLocal): + inp[0][i] = np.random.uniform(0.0, 2.0) + inp0[0][i] = inp[0][i] + out[0][i] = np.random.uniform(-2.0, 2.0) + out0[0][i] = out[0][i] +pos = dataBase.fluidPosition + +#------------------------------------------------------------------------------- +# Construct boundary conditions. +#------------------------------------------------------------------------------- +limits = [x0, x1] +bounds = [] +assert(not (reflecting and periodic)) +if reflecting or periodic: + for d in range(dimension): + planes = [] + for n in range(2): + point = Vector.zero + point[d] = limits[n] + normal = Vector.zero + normal[d] = 1.0 if n == 0 else -1.0 + planes.append(Plane(point, normal)) + if reflecting: + bounds.append(ReflectingBoundary(planes[n])) + if periodic: + bounds.append(PeriodicBoundary(planes[0], planes[1])) +if mpi.procs > 1: + bounds.append(TreeDistributedBoundary.instance()) +output("bounds") + +#------------------------------------------------------------------------------- +# Iterate h +#------------------------------------------------------------------------------- +method = SPHSmoothingScale(IdealH, WT) +iterateIdealH(dataBase, + [method], + bounds, + 100, # max h iterations + 1.e-4) # h tolerance + +#------------------------------------------------------------------------------- +# Finish setting up connectivity +#------------------------------------------------------------------------------- +nodes.numGhostNodes = 0 +nodes.neighbor().updateNodes() +for bc in bounds: + bc.setAllGhostNodes(dataBase) + bc.finalizeGhostBoundary() + nodes.neighbor().updateNodes() +dataBase.updateConnectivityMap() +for bc in bounds: + for fl in [inp, out, inp0, out0]: + bc.applyFieldListGhostBoundary(fl) +for bc in bounds: + bc.finalizeGhostBoundary() + +#------------------------------------------------------------------------------- +# Create map using generated points +#------------------------------------------------------------------------------ +flatConnectivity = FlatConnectivity() +flatConnectivity.computeIndices(dataBase) +flatConnectivity.computeGlobalIndices(dataBase, bounds) +flatConnectivity.computeBoundaryInformation(dataBase, bounds) +output("flatConnectivity") +numGlobal = flatConnectivity.numGlobalNodes() +output("numGlobal") + +#------------------------------------------------------------------------------- +# Test functions +#------------------------------------------------------------------------------ +maxErr = 0.0 +def norm(s1, s2, eps = 1.0e-16): + v1 = np.array(s1) + v2 = np.array(s2) + return np.amax(2.0 * np.abs(v1 - v2) / (np.abs(v1) + np.abs(v2) + eps)) +def check(s1, s2, label, tol = tolerance, eps = 1.0e-16): + err = norm(s1, s2, eps) + global maxErr + if err > maxErr: + maxErr = err + if err > tol: + message = "{}\n\tcalculated={} does not agree with\n\texpected={}\n\terr={}".format(label, s1, s2, err) + if ignoreErr: + print(message) + else: + raise ValueError(message) + +#------------------------------------------------------------------------------- +# Generate matrix +#------------------------------------------------------------------------------ +globRowInd = [] +numColsPerRow = [] +globColInd = [] +colVal = [[] for t in range(numTests)] +multDirect = [np.zeros(numLocal) for t in range(numTests)] +if numpyTest: + spRow = [] + spCol = [] + spDat = [[] for t in range(numTests)] +for i in range(numLocal): + # Also check unique neighbors as a bonus since this function is designed for matrices + globali = flatConnectivity.localToGlobal(i) + numNeighbors = flatConnectivity.numNonConstNeighbors(i) + numConst = flatConnectivity.numConstNeighbors(i) + numUnique, localCols, globalCols, constCols, indexMap = flatConnectivity.uniqueNeighborIndices(i) + check(localCols, flatConnectivity.nonConstNeighborIndices(i), "local cols") + check(globalCols, flatConnectivity.globalNeighborIndices(i), "global cols") + if numConst > 0: + check(constCols, flatConnectivity.constNeighborIndices(i), "const cols") + + # Set row information + numColsPerRow.append(numUnique) + globRowInd.append(globali) + + # Initialize the column data + firstInd = len(globColInd) + globColInd.extend([-1 for j in range(numUnique)]) + for t in range(numTests): + colVal[t].extend([0.0 for j in range(numUnique)]) + + for j in range(numNeighbors): + ind = localCols[j] + globInd = globalCols[j] + uniqueInd = indexMap[j] + flatInd = firstInd + uniqueInd # Index for globColInd and colVal + globColInd[flatInd] = globInd + if numpyTest: + spRow.append(globali) + spCol.append(globInd) + for t in range(numTests): + val = -1.0 * np.random.uniform(1.5, 2.5) if ind == i else np.random.uniform(0.5, 1.5) / numNeighbors + colVal[t][flatInd] += val + multDirect[t][i] += val * inp0[0][ind] + if numpyTest: + spDat[t].append(val) +if numpyTest: + import scipy.sparse as sps + globDat = [mpi.allreduce(spDat[t]) for t in range(numTests)] + globRow = mpi.allreduce(spRow) + globCol = mpi.allreduce(spCol) + npMat = [sps.coo_matrix((globDat[t], (globRow, globCol)), shape=(numGlobal, numGlobal)).tocsr() for t in range(numTests)] + lhsVec = [inp0[0][i] for i in range(numLocal)] + lhsVec = np.array(mpi.allreduce(lhsVec)) + multNumpy = [npMat[t] @ lhsVec for t in range(numTests)] + +#------------------------------------------------------------------------------- +# Initialize solver +#------------------------------------------------------------------------------ +options = HypreOptions() +options.maxNumberOfIterations = 1000 +options.saveLinearSystem = saveSystem +if preconditioner == "amg": + options.preconditionerType = HypreOptions.AMGPreconditioner +elif preconditioner == "ilu": + options.preconditionerType = HypreOptions.ILUPreconditioner + options.factorLevelILU = 5 +else: + options.preconditionerType = HypreOptions.NoPreconditioner +solver = HypreLinearSolver(options); +output("solver") + +start = time.perf_counter() +firstGlobInd = flatConnectivity.firstGlobalIndex() +mpi.barrier() +solver.initialize(numLocal, firstGlobInd, vector_of_unsigned(numColsPerRow)) +hypre_init_time = mpi.allreduce(time.perf_counter() - start, mpi.MAX) +output("hypre_init_time") + +#------------------------------------------------------------------------------- +# Run tests +#------------------------------------------------------------------------------ +for t in range(numTests): + mpi.barrier() + start = time.perf_counter() + solver.beginFill() + solver.setMatRows(numLocal, + vector_of_unsigned(numColsPerRow), + vector_of_unsigned(globRowInd), + vector_of_unsigned(globColInd), + vector_of_double(colVal[t])) + hypre_fill_time = mpi.allreduce(time.perf_counter() - start, mpi.MAX) + output("hypre_fill_time") + + mpi.barrier() + start = time.perf_counter() + solver.assemble() + hypre_assemble_time = mpi.allreduce(time.perf_counter() - start, mpi.MAX) + output("hypre_assemble_time") + + mpi.barrier() + start = time.perf_counter() + solver.finalize() + hypre_finalize_time = mpi.allreduce(time.perf_counter() - start, mpi.MAX) + output("hypre_finalize_time") + + #------------------------------------------------------------------------------- + # Multiply through to get RHS + #------------------------------------------------------------------------------ + # output("inp0[0][0]") + + lhsVec = vector_of_double([inp0[0][i] for i in range(numLocal)]) + mpi.barrier() + solver.set(LinearSolver.LHS, numLocal, firstGlobInd, lhsVec) + mpi.barrier() + solver.multiply() + + rhsVec = solver.get(LinearSolver.RHS, numLocal, firstGlobInd) + for i in range(numLocal): + out[0][i] = rhsVec[i] + + for i in range(numLocal): + if abs(out[0][i] - out0[0][i]) < 1e-16: + raise ValueError("output has not been changed") + + for bc in bounds: + for fl in [inp, out]: + bc.applyFieldListGhostBoundary(fl) + for bc in bounds: + bc.finalizeGhostBoundary() + + #------------------------------------------------------------------------------- + # Check that the multiplication result is as expected + #------------------------------------------------------------------------------ + if numpyTest: + firstGlobal = flatConnectivity.localToGlobal(0) + for i in range(numLocal): + globali = flatConnectivity.localToGlobal(i) + if abs(multNumpy[t][globali] - multDirect[t][i]) > tolerance: + print("warning: numpy and direct mult differ: np={}, dir={}".format(multNumpy[t][globali], multDirect[t][globali])) + if abs(multNumpy[t][globali] - out[0][i]) > tolerance: + raise ValueError("numpy comparison failed for test {}, point {}: exp={}, calc={}".format(t, i, multDirect[t][i], out[0][i])) + print("numpy comparison successful: {} -> {}, {}".format(lhsVec[0], multDirect[t][0], multNumpy[t][firstGlobal])) + + for i in range(numLocal): + if abs(multDirect[t][i] - out[0][i]) > tolerance: + raise ValueError("multiplication incorrect for test {}, point {}: exp={}, calc={}".format(t, i, multDirect[t][i], out[0][i])) + + #------------------------------------------------------------------------------- + # Solve the linear system to get back the LHS + #------------------------------------------------------------------------------ + if exactGuess: + for i in range(numLocal): + inp[0][i] = inp0[0][i] + else: + for i in range(numLocal): + inp[0][i] = np.random.uniform(0.0, 1.0) + + for i in range(numSolves): + lhsVec = vector_of_double([inp[0][i] for i in range(numLocal)]) + rhsVec = vector_of_double([out[0][i] for i in range(numLocal)]) + + start = time.perf_counter() + mpi.barrier() + solver.set(LinearSolver.LHS, numLocal, firstGlobInd, lhsVec) + solver.set(LinearSolver.RHS, numLocal, firstGlobInd, rhsVec) + mpi.barrier() + solver.solve() + hypre_solve_time = mpi.allreduce(time.perf_counter() - start, mpi.MAX) + output("hypre_solve_time") + + mpi.barrier() + lhsVec = solver.get(LinearSolver.LHS, numLocal, firstGlobInd) + for i in range(numLocal): + inp[0][i] = lhsVec[i] + + # output("inp[0][0]") + + #------------------------------------------------------------------------------- + # Check the results + #------------------------------------------------------------------------------ + for i in range(numLocal): + if abs(inp[0][i] - inp0[0][i]) > tolerance: + raise ValueError("error too high for test {} point {}: begin={}, end={}".format(t + 1, i, inp0[0][i], inp[0][i])) + print("test {} passed".format(t+1)) +