Add LinearSolver interface

CHANGELOG.md

@@ -68,6 +68,9 @@
 - Added IDA option to suppress algebraic variables in local error tests.
 - Removed `COO_Matrix` class.
 - Added portable `Vector` class and policy-based memory utilities.
+- Add new `LinearSolver` interface for linear solvers.
+- Added new `Rosenbrock` integrator.
+- Clarified naming conventions for macros.
 
 ## v0.1
 

CMakeLists.txt

@@ -33,6 +33,9 @@ option(GridKit_ENABLE_SUNDIALS "Enable SUNDIALS support" OFF)
 # Enzyme support is disabled by default
 option(GridKit_ENABLE_ENZYME "Enable automatic differentiation with Enzyme" OFF)
 
+# ReSolve support is disabled by default
+option(GridKit_ENABLE_RESOLVE "Enable ReSolve support" OFF)
+
 # Enable the address sanitizer
 option(GridKit_ENABLE_ASAN "Enable the address sanitizer" OFF)
 
@@ -52,6 +55,7 @@ list(
   IPOPT
   SUNDIALS
   ENZYME
+  RESOLVE
   ASAN
   UBSAN
   OPENMP)
@@ -147,6 +151,20 @@ if(GRIDKIT_ENABLE_THREADS)
   find_package(Threads REQUIRED)
 endif()
 
+if(GRIDKIT_ENABLE_RESOLVE)
+  find_package(
+    ReSolve
+    CONFIG
+    PATHS
+    ${ReSolve_DIR}
+    ${ReSolve_DIR}/share/resolve/cmake
+    ENV
+    LD_LIBRARY_PATH
+    ENV
+    DYLD_LIBRARY_PATH
+    REQUIRED)
+endif()
+
 # Macro that adds libraries
 include(GridkitAddLibrary)
 

CONTRIBUTING.md

@@ -285,6 +285,10 @@ name. Use all caps (screaming snake case).
    constexpr double EXP = 2.7183      // Yes
 ```
 
+### Macros
+
+Macros should use all caps, same as constants.
+
 ### Enums (enumerated types)
 
 Always define `enum`s inside `GridKit` namespace. The `enum` name should

GridKit/LinearAlgebra/CMakeLists.txt

@@ -1,3 +1,4 @@
 add_subdirectory(SparseMatrix)
 add_subdirectory(DenseMatrix)
 add_subdirectory(Vector)
+add_subdirectory(Solver)

GridKit/LinearAlgebra/DenseMatrix/DenseMatrix.hpp

@@ -6,18 +6,18 @@
 #include <limits>
 #include <vector>
 
-/**
- * @brief Class to provide dense matrices.
- *
- * @deprecated This class is deprecated. Use Vector(n, k) instead.
- *
- * This is intended for small matrices that store model Jacobians to be subsequently copied
- * into large sparse matrices.
- */
 namespace GridKit
 {
   namespace LinearAlgebra
   {
+    /**
+     * @brief Class to provide dense matrices.
+     *
+     * @deprecated This class is deprecated. Use Vector(n, k) instead.
+     *
+     * This is intended for small matrices that store model Jacobians to be subsequently copied
+     * into large sparse matrices.
+     */
     template <typename RealT, typename IdxT>
     class DenseMatrix
     {

GridKit/LinearAlgebra/Solver/CMakeLists.txt

@@ -0,0 +1,7 @@
+if(GRIDKIT_ENABLE_RESOLVE)
+  gridkit_add_library(
+    resolve_system_solver
+    SOURCES ResolveSystemSolver.cpp
+    HEADERS ResolveSystemSolver.hpp
+    LINK_LIBRARIES PRIVATE ReSolve::ReSolve)
+endif()

GridKit/LinearAlgebra/Solver/LinearSolver.hpp

@@ -0,0 +1,55 @@
+#pragma once
+
+#include <GridKit/LinearAlgebra/SparseMatrix/CsrMatrix.hpp>
+#include <GridKit/ScalarTraits.hpp>
+
+#include "GridKit/LinearAlgebra/Vector/Vector.hpp"
+
+namespace GridKit
+{
+  namespace LinearAlgebra
+  {
+
+    /**
+     * @brief An interface for linear solvers to be used in GridKit, such as in \ref Integrator::Rosenbrock.
+     *
+     */
+    template <class ScalarT, typename IdxT>
+    class LinearSolver
+    {
+    public:
+      using RealT = GridKit::ScalarTraits<ScalarT>::RealT;
+
+      /**
+       * @brief Configure the solver by analyzing the given matrix. The sparsity pattern of the given matrix must be set, and future calls
+       * to other methods will use this sparsity pattern, as well as the data pointer of the given matrix.
+       *
+       * @param matrix The matrix to configure the solver with.
+       * @return int An error code, or 0 if none.
+       */
+      virtual int configureSolver(GridKit::LinearAlgebra::CsrMatrix<RealT, IdxT>& matrix)                                       = 0;
+      /**
+       * @brief Setup the solver with matrix data.
+       *
+       * @param reuse_factors If true, factors from a previous call to `setupSolver()` are reused.
+       * @pre \ref configureSolver() must be called first and every time the sparsity pattern of the matrix changes.
+       * @pre If `reuse_factors` is true, then `setupSolver` must have been called once before with the same sparsity pattern.
+       * @return int An error code, or 0 if none.
+       */
+      virtual int setupSolver(bool reuse_factors = false)                                                                       = 0;
+      /**
+       * @brief Perform a linear solve using the configured matrix.
+       *
+       * @param rhs The right-hand-side vector (input). No data is changed.
+       * @param lhs The left-hand-side vector (output).
+       * @pre \ref setupSolver must be called first and every time the data of the matrix changes.
+       * @return int An error code, or 0 if none.
+       *
+       * @todo The data of the `rhs` parameter doesn't change, but we can't mark as cosntant because ReSolve currently requires
+       * a non-`const` pointer for vectors, even if it will only read from that vector. In the future this might change, so this
+       * can be updated to be `const`.
+       */
+      virtual int solve(GridKit::LinearAlgebra::Vector<ScalarT, IdxT>& rhs, GridKit::LinearAlgebra::Vector<ScalarT, IdxT>& lhs) = 0;
+    };
+  } // namespace LinearAlgebra
+} // namespace GridKit

GridKit/LinearAlgebra/Solver/ResolveSystemSolver.cpp

@@ -0,0 +1,105 @@
+#include "ResolveSystemSolver.hpp"
+
+#include "GridKit/MemoryUtilities/MemoryUtils.hpp"
+#include <resolve/vector/Vector.hpp>
+
+namespace GridKit
+{
+  namespace LinearAlgebra
+  {
+    /**
+     * @brief Helper function to convert GridKit `MemorySpace` into ReSolve `MemorySpace`.
+     *
+     */
+    ReSolve::memory::MemorySpace memorySpaceAsResolve(memory::MemorySpace memspace)
+    {
+      switch (memspace)
+      {
+      case memory::HOST:
+        return ReSolve::memory::HOST;
+      case memory::DEVICE:
+        return ReSolve::memory::DEVICE;
+      default:
+        throw "Memory space not supported";
+      }
+    }
+
+    template <class ScalarT, typename IdxT>
+    ResolveSystemSolver<ScalarT, IdxT>::ResolveSystemSolver(ReSolve::SystemSolver& lin_solver, GridKit::memory::MemorySpace memspace)
+      : lin_solver_(lin_solver), memspace_(memorySpaceAsResolve(memspace))
+    {
+    }
+
+    /**
+     * @brief Creates a new ReSolve matrix with the dimensions of the given matrix, then sets its data pointers to those of the given matrix.
+     * Sets up the linear solver to use this new ReSolve matrix.
+     *
+     * @todo Right now preconditioning doesn't work. There should be a ReSolve PR soon for this.
+     *
+     */
+    template <class ScalarT, typename IdxT>
+    int ResolveSystemSolver<ScalarT, IdxT>::configureSolver(GridKit::LinearAlgebra::CsrMatrix<RealT, IdxT>& matrix)
+    {
+      matrix_ = std::make_unique<ReSolve::matrix::Csr>(matrix.getNumRows(), matrix.getNumColumns(), matrix.getNnz());
+
+      if (int err = matrix_->setDataPointers(matrix.getRowData(), matrix.getColData(), matrix.getValues(), memspace_))
+        return err;
+
+      if (int err = lin_solver_.setMatrix(matrix_.get()))
+        return err;
+
+      if (int err = lin_solver_.analyze())
+        return err;
+
+      // TODO: This function needs to be called to properly use a preconditioner in ReSolve (if there is any), but currently will error
+      // if there is no preconditioner configured. Once we can detect if a preconditioner is configured, we can restore this functionality.
+      // Also, we will always want to use *right* preconditioning.
+      // BUBBLE_FAIL(lin_solver_.preconditionerSetup("right"));
+
+      return 0;
+    }
+
+    template <class ScalarT, typename IdxT>
+    int ResolveSystemSolver<ScalarT, IdxT>::setupSolver(bool reuse_factors)
+    {
+      if (reuse_factors)
+      {
+        if (int err = lin_solver_.refactorize())
+          return err;
+      }
+      else
+      {
+        if (int err = lin_solver_.factorize())
+          return err;
+      }
+
+      return 0;
+    }
+
+    /**
+     * @brief Perform the solver by creating two new ReSolve vectors and setting their data to point to their GridKit counterparts.
+     * Then ask ReSolve to perform the solve. Since the matrix should be pointing at the correct GridKit buffers and the vectors
+     * as well, this solve should correctly fill `lhs`'s data buffer.
+     *
+     */
+    template <class ScalarT, typename IdxT>
+    int ResolveSystemSolver<ScalarT, IdxT>::solve(GridKit::LinearAlgebra::Vector<ScalarT, IdxT>& rhs, GridKit::LinearAlgebra::Vector<ScalarT, IdxT>& lhs)
+    {
+      ReSolve::vector::Vector resolve_rhs(rhs.getSize());
+      ReSolve::vector::Vector resolve_lhs(lhs.getSize());
+
+      if (int err = resolve_rhs.setData(rhs.getData(), memspace_))
+        return err;
+
+      if (int err = resolve_lhs.setData(lhs.getData(), memspace_))
+        return err;
+
+      if (int err = lin_solver_.solve(&resolve_rhs, &resolve_lhs))
+        return err;
+
+      return 0;
+    }
+
+    template class ResolveSystemSolver<double, int>;
+  } // namespace LinearAlgebra
+} // namespace GridKit

GridKit/LinearAlgebra/Solver/ResolveSystemSolver.hpp

@@ -0,0 +1,45 @@
+#pragma once
+
+#include <memory>
+
+#include <GridKit/LinearAlgebra/Solver/LinearSolver.hpp>
+
+#include <resolve/MemoryUtils.hpp>
+#include <resolve/SystemSolver.hpp>
+#include <resolve/matrix/Csr.hpp>
+
+namespace GridKit
+{
+  namespace LinearAlgebra
+  {
+    template <class ScalarT, typename IdxT>
+    class ResolveSystemSolver : public LinearSolver<ScalarT, IdxT>
+    {
+    public:
+      using RealT = LinearSolver<ScalarT, IdxT>::RealT;
+
+      ResolveSystemSolver(ReSolve::SystemSolver& lin_solver, GridKit::memory::MemorySpace memspace = GridKit::memory::HOST);
+
+      int configureSolver(GridKit::LinearAlgebra::CsrMatrix<RealT, IdxT>& matrix) final;
+      int setupSolver(bool reuse_factors = false) final;
+      int solve(GridKit::LinearAlgebra::Vector<ScalarT, IdxT>& rhs, GridKit::LinearAlgebra::Vector<ScalarT, IdxT>& lhs) final;
+
+    private:
+      /**
+       * @brief The ReSolve matrix used for the linear solve. Its data pointers should just point to a \ref GridKit::LinearAlgebra::CsrMatrix.
+       *
+       */
+      std::unique_ptr<ReSolve::matrix::Csr> matrix_;
+      /**
+       * @brief The ReSolve linear solver to use.
+       *
+       */
+      ReSolve::SystemSolver&                lin_solver_;
+      /**
+       * @brief The memspace the data is stored in.
+       *
+       */
+      ReSolve::memory::MemorySpace          memspace_;
+    };
+  } // namespace LinearAlgebra
+} // namespace GridKit

GridKit/Solver/Dynamic/CMakeLists.txt

@@ -5,13 +5,11 @@
 #    - Cameron Rutherford <cameron.rutherford@pnnl.gov>
 #]]
 
-set(_install_headers DynamicSolver.hpp Ida.hpp)
-
 if(GRIDKIT_ENABLE_SUNDIALS_SPARSE)
   gridkit_add_library(
     solvers_dyn
     SOURCES Ida.cpp
-    HEADERS ${_install_headers}
+    HEADERS Ida.hpp DynamicSolver.hpp
     LINK_LIBRARIES
       PUBLIC
       SUNDIALS::nvecserial
@@ -29,7 +27,7 @@ else()
   gridkit_add_library(
     solvers_dyn
     SOURCES Ida.cpp
-    HEADERS ${_install_headers}
+    HEADERS Ida.hpp DynamicSolver.hpp
     LINK_LIBRARIES
       PUBLIC
       SUNDIALS::nvecserial
@@ -42,3 +40,26 @@ else()
       PUBLIC
       GridKit::sparse_matrix)
 endif()
+
+if(GRIDKIT_ENABLE_RESOLVE)
+  gridkit_add_library(
+    rosenbrock
+    SOURCES Rosenbrock.cpp RosenbrockTableaus.cpp
+    HEADERS Rosenbrock.hpp
+            Native/ErrorNorm.hpp
+            Native/StepControl.hpp
+            Native/StepController.hpp
+    LINK_LIBRARIES
+      PRIVATE
+      ReSolve::ReSolve
+      PUBLIC
+      GridKit::dense_vector
+      PUBLIC
+      GridKit::definitions
+      PUBLIC
+      GridKit::utilities_logger
+      PUBLIC
+      GridKit::sparse_matrix)
+endif()
+
+add_subdirectory(Native)

GridKit/Solver/Dynamic/Native/AdaptiveStep.cpp

@@ -0,0 +1,33 @@
+#include "AdaptiveStep.hpp"
+
+#include <algorithm>
+#include <cmath>
+
+#include <GridKit/Constants.hpp>
+
+namespace AnalysisManager
+{
+  namespace NativeDynamicSolver
+  {
+    /**
+     * @brief Standard textbook adaptive controller. Accept if `err <= 1` and use
+     *
+     * \f[h_{new} = h * \min \left\{fac_{max}, \max\left\{fac_{min}, fac_{scale} \cdot e ^{-1/p}\right\}\right\}.\f]
+     *
+     */
+    template <typename RealT>
+    StepControl<RealT> AdaptiveStep<RealT>::nextStep(RealT err, StepControl<RealT> prev_step, uint8_t method_order)
+    {
+      StepControl<RealT> next_step = prev_step;
+
+      double h_mult = std::min(params_.fac_max_, std::max(params_.fac_scale_ * std::pow(err, GridKit::MINUS_ONE<RealT> / method_order), params_.fac_min_));
+
+      next_step.accept_     = err <= 1;
+      next_step.step_size_ *= h_mult;
+
+      return next_step;
+    }
+
+    template class AdaptiveStep<double>;
+  } // namespace NativeDynamicSolver
+} // namespace AnalysisManager