ABBA GMRES Solver

resolve/CMakeLists.txt

@@ -22,6 +22,7 @@ set(ReSolve_SRC
     SystemSolver.cpp
     Preconditioner.cpp
     PreconditionerLU.cpp
+    PreconditionerABBA.cpp
 )
 
 set(ReSolve_KLU_SRC LinSolverDirectKLU.cpp)
@@ -53,6 +54,7 @@ set(ReSolve_HEADER_INSTALL
     MemoryUtils.hpp
     Preconditioner.hpp
     PreconditionerLU.hpp
+    PreconditionerABBA.hpp
 )
 
 set(ReSolve_KLU_HEADER_INSTALL LinSolverDirectKLU.hpp)

resolve/LinSolverIterativeFGMRES.cpp

@@ -128,8 +128,18 @@ namespace ReSolve
     vec_Z_->setToZero(memspace_);
     vec_V_->setToZero(memspace_);
 
+    // V[0] = b-A*x_0
     rhs->copyDataTo(vec_V_->getData(memspace_), 0, memspace_);
     matrix_handler_->matvec(A_, x, vec_V_, &MINUS_ONE, &ONE, memspace_);
+
+    // Calculate r_0 = ||B(b - Ax_0)|| for BAGMRES
+    if (preconditioner_->getSide() == "left")
+    {
+      vector_type temp_vec(n_);
+      temp_vec.copyDataFrom(vec_V_->getData(0, memspace_), memspace_, memspace_);
+      matrix_handler_->matvec(preconditioner_->getPrecMatrix(), &temp_vec, vec_V_, &ONE, &ZERO, memspace_);
+    }
+
     rnorm = 0.0;
     bnorm = vector_handler_->dot(rhs, rhs, memspace_);
     rnorm = vector_handler_->dot(vec_V_, vec_V_, memspace_);
@@ -197,14 +207,29 @@ namespace ReSolve
         {
           vec_z.setData(vec_Z_->getData(0, memspace_), memspace_);
         }
-        this->precV(&vec_v, &vec_z);
-        mem_.deviceSynchronize();
+        // this->precV(&vec_v, &vec_z);
+        // mem_.deviceSynchronize();
 
-        // V_{i+1}=A*Z_i
+        // V_{i+1}=A*Z_i or B * Z_i for BAGMRES
 
-        vec_v.setData(vec_V_->getData(i + 1, memspace_), memspace_);
+        // vec_v.setData(vec_V_->getData(i + 1, memspace_), memspace_);
 
-        matrix_handler_->matvec(A_, &vec_z, &vec_v, &ONE, &ZERO, memspace_);
+        // Right preconditioned
+        if (preconditioner_->getSide() == "right")
+        {
+          this->precV(&vec_v, &vec_z);
+          mem_.deviceSynchronize();
+          vec_v.setData(vec_V_->getData(i + 1, memspace_), memspace_);
+          matrix_handler_->matvec(A_, &vec_z, &vec_v, &ONE, &ZERO, memspace_);
+        }
+        // Left Preconditioned
+        else
+        {
+          matrix_handler_->matvec(A_, &vec_v, &vec_z, &ONE, &ZERO, memspace_);
+          mem_.deviceSynchronize();
+          vec_v.setData(vec_V_->getData(i + 1, memspace_), memspace_);
+          this->precV(&vec_z, &vec_v);
+        }
 
         // orthogonalize V[i+1], form a column of h_H_
 
@@ -284,12 +309,18 @@ namespace ReSolve
           vec_v.setData(vec_V_->getData(j, memspace_), memspace_);
           vector_handler_->axpy(&h_rs_[j], &vec_v, &vec_z, memspace_);
         }
-        // now multiply d_Z by precon
-
-        vec_v.setData(vec_V_->getData(memspace_), memspace_);
-        this->precV(&vec_z, &vec_v);
-        // and add to x
-        vector_handler_->axpy(&ONE, &vec_v, x, memspace_);
+        // now multiply d_Z by precon for left preconditioned
+        if (preconditioner_->getSide() == "right")
+        {
+          this->precV(&vec_z, &vec_v);
+          // and add to x
+          vector_handler_->axpy(&ONE, &vec_v, x, memspace_);
+        }
+        // Directly use vec_z for right preconditioned
+        else
+        {
+          vector_handler_->axpy(&ONE, &vec_z, x, memspace_);
+        }
       }
 
       /* test solution */

resolve/Preconditioner.cpp

@@ -22,4 +22,19 @@ namespace ReSolve
     return 1;
   }
 
+  std::string Preconditioner::getSide()
+  {
+    return side_;
+  }
+
+  /**
+   * @brief Used to get the preconditioning matrix for Preconditioner Matvec
+   *
+   * Should not be called unless using PreconditionerMatvec
+   */
+  matrix::Sparse* Preconditioner::getPrecMatrix()
+  {
+    return nullptr;
+  }
+
 } // namespace ReSolve

resolve/Preconditioner.hpp

@@ -4,6 +4,9 @@
  * @brief  Declaration of preconditioner base class.
  *
  */
+
+#include <string>
+
 #pragma once
 
 namespace ReSolve
@@ -32,8 +35,13 @@ namespace ReSolve
     Preconditioner();
     virtual ~Preconditioner();
 
-    virtual int setup(matrix_type* A)                   = 0;
-    virtual int apply(vector_type* rhs, vector_type* x) = 0;
-    virtual int reset(matrix_type* /* A */);
+    virtual int          setup(matrix_type* A)                   = 0;
+    virtual int          apply(vector_type* rhs, vector_type* x) = 0;
+    virtual int          reset(matrix_type* /* A */);
+    virtual std::string  getSide();       // Gets the preconditioning side
+    virtual matrix_type* getPrecMatrix(); // Get the preconditioner matrix
+
+  private:
+    std::string side_ = "right"; // Right preconditioning by default
   };
 } // namespace ReSolve

resolve/PreconditionerABBA.cpp

@@ -0,0 +1,109 @@
+/**
+ * @file   PreconditionerABBA.cpp
+ * @author Jeffery Zhang (jefferyz@vt.edu)
+ * @brief  Declaration of PreconditionerMatrix class.
+ */
+
+#include "PreconditionerABBA.hpp"
+
+namespace ReSolve
+{
+  using out = io::Logger;
+
+  /**
+   * @brief Constructor for PreconditionerMatrix.
+   *
+   * @param[in] A - Pointer to the forward operator
+   * @param[in] matrix_handler - Pointer to the matrix handler
+   */
+  PreconditionerABBA::PreconditionerABBA(matrix::Sparse* A, MatrixHandler* matrix_handler)
+  {
+    A_              = A;
+    matrix_handler_ = matrix_handler;
+    setMemorySpace();
+  }
+
+  /**
+   * @brief Destructor for PreconditionerLU
+   */
+  PreconditionerABBA::~PreconditionerABBA()
+  {
+  }
+
+  /**
+   * @brief Set the preconditioning matrix
+   *
+   * @param[in] B - Pointer to the preconditioning matrix
+   */
+  int PreconditionerABBA::setup(matrix::Sparse* B)
+  {
+    B_ = B;
+    return 0;
+  }
+
+  /**
+   * @brief Get the preconditioner matrix
+   *
+   * Necessary for some calculations
+   *
+   */
+  matrix::Sparse* PreconditionerABBA::getPrecMatrix()
+  {
+    return B_;
+  }
+
+  /**
+   * @brief Setter for the side
+   *
+   * @param[in] The new side
+   */
+  int PreconditionerABBA::setSide(std::string side)
+  {
+    if (side == "left" || side == "right")
+    {
+      side_ = side;
+      return 0;
+    }
+    out::error() << "Choose either left (BA) or right (AB)\n";
+    return 1;
+  }
+
+  /**
+   * @brief getter for the side
+   */
+  std::string PreconditionerABBA::getSide()
+  {
+    return side_;
+  }
+
+  /**
+   * @brief Applies the preconditioner depending on the side
+   *
+   * @param[in] rhs - Right-hand-side vector
+   * @param[in] x   - Solution vector
+   *
+   * @return int 0 if successful, 1 if fails
+   */
+  int PreconditionerABBA::apply(vector_type* rhs, vector_type* x)
+  {
+    using namespace constants;
+    matrix_handler_->matvec(B_, rhs, x, &ONE, &ZERO, memspace_);
+    return 0;
+  }
+
+  void PreconditionerABBA::setMemorySpace()
+  {
+    bool is_matrix_handler_cuda = matrix_handler_->getIsCudaEnabled();
+    bool is_matrix_handler_hip  = matrix_handler_->getIsHipEnabled();
+
+    if (is_matrix_handler_cuda || is_matrix_handler_hip)
+    {
+      memspace_ = memory::DEVICE;
+    }
+    else
+    {
+      memspace_ = memory::HOST;
+    }
+  }
+
+} // namespace ReSolve

resolve/PreconditionerABBA.hpp

@@ -0,0 +1,62 @@
+/**
+ * @file   PreconditionerABBA.hpp
+ * @author Jeffery Zhang (jefferyz@vt.edu)
+ * @brief  Declaration of left and right preconditioner class
+ */
+
+#include <string>
+
+#include "Common.hpp"
+#include <resolve/MemoryUtils.hpp>
+#include <resolve/Preconditioner.hpp>
+#include <resolve/matrix/MatrixHandler.hpp>
+
+namespace ReSolve
+{
+
+  namespace matrix
+  {
+    class Sparse;
+  } // namespace matrix
+
+  namespace vector
+  {
+    class Vector;
+  } // namespace vector
+
+  // Forward declaration of MatrixHandler class
+  class MatrixHandler;
+
+  /**
+   * @brief Class allows for a user specified matrix to be used for preconditioning
+   * Allows user to switch between right preconditioning (AB) and
+   * left preconditioning (BA)
+   *
+   * @author Jeffery Zhang (jefferyz@vt.edu)
+   *
+   */
+  class PreconditionerABBA : public Preconditioner
+  {
+  public:
+    using vector_type = vector::Vector;
+    using matrix_type = matrix::Sparse;
+
+    PreconditionerABBA(matrix_type* A, MatrixHandler* matrix_handler);
+    ~PreconditionerABBA();
+
+    int          apply(vector_type* rhs, vector_type* x) override; // Applies preconditioning
+    int          setup(matrix_type*) override;
+    std::string  getSide() override;
+    matrix_type* getPrecMatrix() override;  // Used to get the preconditioning matrix for calculation of initial residual for BAGMRES
+    int          setSide(std::string side); // Changes the preconditioning side for BAGMRES
+
+  private:
+    void setMemorySpace();
+
+    matrix_type*        A_{nullptr};
+    matrix_type*        B_{nullptr};
+    std::string         side_ = "right"; // Defaults to ABGMRES
+    MatrixHandler*      matrix_handler_{nullptr};
+    memory::MemorySpace memspace_;
+  };
+} // namespace ReSolve

tests/functionality/CMakeLists.txt

@@ -17,6 +17,10 @@ target_link_libraries(rand_gmres_test.exe PRIVATE ReSolve)
 add_executable(sys_rand_gmres_test.exe testSysGmres.cpp)
 target_link_libraries(sys_rand_gmres_test.exe PRIVATE ReSolve)
 
+# Test for ABBA GMRES
+add_executable(trips_py_test.exe testTripsProblem.cpp)
+target_link_libraries(trips_py_test.exe PRIVATE ReSolve)
+
 if(RESOLVE_USE_KLU)
   # Build KLU+KLU test
   add_executable(klu_klu_test.exe testKlu.cpp)
@@ -48,6 +52,7 @@ endif()
 set(installable_tests version.exe)
 list(APPEND installable_tests rand_gmres_test.exe)
 list(APPEND installable_tests sys_rand_gmres_test.exe)
+list(APPEND installable_tests trips_py_test.exe)
 
 if(RESOLVE_USE_KLU)
   list(APPEND installable_tests klu_klu_test.exe)