Feature/fgmres update

resolve/GramSchmidt.cpp

@@ -172,7 +172,7 @@ namespace ReSolve
       }
       else
       {
-        assert(0 && "Gram-Schmidt failed, vector with ZERO norm\n");
+        std::cout << "Encountered loss of Orthogonality, exiting gracefully at step " << i + 1 << std::endl;
         return 1;
       }
       return 0;
@@ -222,7 +222,7 @@ namespace ReSolve
       }
       else
       {
-        assert(0 && "Gram-Schmidt failed, vector with ZERO norm\n");
+        std::cout << "Encountered loss of Orthogonality, exiting gracefully at step " << i + 1 << std::endl;
         return 1;
       }
       return 0;
@@ -279,7 +279,7 @@ namespace ReSolve
       }
       else
       {
-        assert(0 && "Iterative refinement failed, Krylov vector with ZERO norm\n");
+        std::cout << "Encountered loss of Orthogonality, exiting gracefully at step " << i + 1 << std::endl;
         return 1;
       }
       h_rv = nullptr;
@@ -364,7 +364,7 @@ namespace ReSolve
       }
       else
       {
-        assert(0 && "Iterative refinement failed, Krylov vector with ZERO norm\n");
+        std::cout << "Encountered loss of Orthogonality, exiting gracefully at step " << i + 1 << std::endl;
         return 1;
       }
       h_rv = nullptr;
@@ -395,7 +395,7 @@ namespace ReSolve
       }
       else
       {
-        assert(0 && "Gram-Schmidt failed, vector with ZERO norm\n");
+        std::cout << "Encountered loss of Orthogonality, exiting gracefully at step " << i + 1 << std::endl;
         return 1;
       }
       return 0;

resolve/LinSolverIterativeFGMRES.cpp

@@ -108,37 +108,63 @@ namespace ReSolve
 
     // io::Logger::setVerbosity(io::Logger::EVERYTHING);
 
-    int outer_flag = 1;
-    int notconv    = 1;
-    int i          = 0;
-    int it         = 0;
-    int j          = 0;
-    int k          = 0;
-    int k1         = 0;
+    int outer_flag  = 1;
+    int update_flag = 1;
+    int notconv     = 1;
+    int i           = 0;
+    int it          = 0;
+    int j           = 0;
+    int k           = 0;
+    int k1          = 0;
 
     real_type   t     = 0.0;
     real_type   rnorm = 0.0;
     real_type   bnorm = 0.0;
     real_type   tolrel;
+    real_type   rhsnorm = 0.0;
     vector_type vec_v(n_);
     vector_type vec_z(n_);
-    // V[0] = b-A*x_0
-    // debug
+
+    // Compute the residual
+    vec_R_->setToZero(memspace_);
+    vec_R_->copyDataFrom(rhs, memspace_, memspace_);
+    matrix_handler_->matvec(A_, x, vec_R_, &MINUS_ONE, &ONE, memspace_);
+
+    // Normalizing the RHS
+    bnorm = vector_handler_->dot(vec_R_, vec_R_, memspace_);
+    bnorm = std::sqrt(bnorm);
+    t     = 1 / bnorm;
+    vector_handler_->scal(&t, vec_R_, memspace_);
+
+    // Arnoldi Basis
     vec_Z_->setToZero(memspace_);
     vec_V_->setToZero(memspace_);
 
-    rhs->copyDataTo(vec_V_->getData(memspace_), 0, memspace_);
-    matrix_handler_->matvec(A_, x, vec_V_, &MINUS_ONE, &ONE, memspace_);
-    rnorm = 0.0;
-    bnorm = vector_handler_->dot(rhs, rhs, memspace_);
+    // Initializing Residual and Update
+    vec_Y_->setToZero(memspace_);
+
+    // Computing the first Arnodi basis vector
+    vec_R_->copyDataTo(vec_V_->getData(memspace_), 0, memspace_);
     rnorm = vector_handler_->dot(vec_V_, vec_V_, memspace_);
-    // rnorm = ||V_1||
     rnorm = std::sqrt(rnorm);
-    bnorm = std::sqrt(bnorm);
+
+    // Checking if bnorm > norm of RHS
+    rhsnorm = vector_handler_->dot(rhs, rhs, memspace_);
+    rhsnorm = std::sqrt(rhsnorm);
+    if (bnorm > rhsnorm)
+    {
+      std::cout << "Initial guess is invalid. Quitting iterative refinement" << std::endl;
+      initial_residual_norm_ = bnorm;
+      final_residual_norm_   = bnorm; // Set final norm to initial norm as no work was done.
+      total_iters_           = 0;
+      return 0;
+    }
+
     io::Logger::misc() << "it 0: norm of residual "
                        << std::scientific << std::setprecision(16)
                        << rnorm << " Norm of rhs: " << bnorm << "\n";
-    initial_residual_norm_ = rnorm;
+    initial_residual_norm_ = bnorm;
+
     while (outer_flag)
     {
       if (it == 0)
@@ -169,13 +195,11 @@ namespace ReSolve
         outer_flag             = 0;
         final_residual_norm_   = rnorm;
         initial_residual_norm_ = rnorm;
+        update_flag            = 1;
         total_iters_           = 0;
         break;
       }
 
-      // normalize first vector
-      t = 1.0 / rnorm;
-      vector_handler_->scal(&t, vec_V_, memspace_);
       // initialize norm history
       h_rs_[0] = rnorm;
       i        = -1;
@@ -200,14 +224,18 @@ namespace ReSolve
         mem_.deviceSynchronize();
 
         // V_{i+1}=A*Z_i
-
         vec_v.setData(vec_V_->getData(i + 1, memspace_), memspace_);
-
         matrix_handler_->matvec(A_, &vec_z, &vec_v, &ONE, &ZERO, memspace_);
 
         // orthogonalize V[i+1], form a column of h_H_
 
-        GS_->orthogonalize(n_, vec_V_, h_H_, i);
+        int gs_status = GS_->orthogonalize(n_, vec_V_, h_H_, i);
+
+        if (gs_status != 0) // checking for successful breakdown
+        {
+          notconv = 0; // exiting outer loop after one inner loop iteration
+        }
+
         if (i != 0)
         {
           for (index_type k = 1; k <= i; k++)
@@ -252,7 +280,18 @@ namespace ReSolve
                          << std::scientific << std::setprecision(16)
                          << rnorm << "\n";
       // solve tri system
-      h_rs_[i] = h_rs_[i] / h_H_[i * (restart_ + 1) + i];
+      real_type H_ii_diag = h_H_[i * (restart_ + 1) + i];
+
+      // FIX: Check for division by zero before starting back-substitution
+      if (std::abs(H_ii_diag) < MACHINE_EPSILON)
+      {
+        h_rs_[i] = ZERO; // Last element of solution is set to 0
+      }
+      else
+      {
+        h_rs_[i] = h_rs_[i] / H_ii_diag;
+      }
+
       for (int ii = 2; ii <= i + 1; ii++)
       {
         k  = i - ii + 1;
@@ -262,7 +301,17 @@ namespace ReSolve
         {
           t -= h_H_[j * (restart_ + 1) + k] * h_rs_[j];
         }
-        h_rs_[k] = t / h_H_[k * (restart_ + 1) + k];
+
+        // Checking for division by zero inside the loop as well
+        real_type H_diag_k = h_H_[k * (restart_ + 1) + k];
+        if (std::abs(H_diag_k) < MACHINE_EPSILON)
+        {
+          h_rs_[k] = ZERO;
+        }
+        else
+        {
+          h_rs_[k] = t / H_diag_k;
+        }
       }
 
       // get solution
@@ -271,7 +320,7 @@ namespace ReSolve
         for (j = 0; j <= i; j++)
         {
           vec_z.setData(vec_Z_->getData(j, memspace_), memspace_);
-          vector_handler_->axpy(&h_rs_[j], &vec_z, x, memspace_);
+          vector_handler_->axpy(&h_rs_[j], &vec_z, vec_Y_, memspace_);
         }
       }
       else
@@ -287,20 +336,20 @@ namespace ReSolve
 
         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_);
+        // and add to update
+        vector_handler_->axpy(&ONE, &vec_v, vec_Y_, memspace_);
       }
 
-      /* test solution */
+      /* test update */
 
       if (rnorm <= tolrel || it >= maxit_)
       {
         // rnorm_aux = rnorm;
         outer_flag = 0;
       }
 
-      rhs->copyDataTo(vec_V_->getData(memspace_), 0, memspace_);
-      matrix_handler_->matvec(A_, x, vec_V_, &MINUS_ONE, &ONE, memspace_);
+      vec_R_->copyDataTo(vec_V_->getData(memspace_), 0, memspace_);
+      matrix_handler_->matvec(A_, vec_Y_, vec_V_, &MINUS_ONE, &ONE, memspace_);
       rnorm = vector_handler_->dot(vec_V_, vec_V_, memspace_);
       // rnorm = ||V_1||
       rnorm = std::sqrt(rnorm);
@@ -314,6 +363,37 @@ namespace ReSolve
                            << rnorm << "\n";
       }
     } // outer while
+
+    if (update_flag)
+    {
+      // renormalizing
+      t = bnorm;
+      vector_handler_->scal(&t, vec_Y_, memspace_);
+      vector_handler_->axpy(&ONE, x, vec_Y_, memspace_);
+      vec_R_->copyDataFrom(rhs, memspace_, memspace_);
+      matrix_handler_->matvec(A_, vec_Y_, vec_R_, &MINUS_ONE, &ONE, memspace_);
+      final_residual_norm_ = vector_handler_->dot(vec_R_, vec_R_, memspace_);
+      final_residual_norm_ = std::sqrt(final_residual_norm_);
+
+      // Compare this with initial and update x accordingly
+      if (final_residual_norm_ < initial_residual_norm_)
+      {
+        std::cout << "Update to intial guess is successful, final residual (solution plus update) "
+                  << std::scientific << std::setprecision(16)
+                  << final_residual_norm_ / rhsnorm << "\n";
+
+        x->copyDataFrom(vec_Y_, memspace_, memspace_);
+      }
+      else
+      {
+        std::cout << "Update to intial guess is not successful, final residual greater than initial residual "
+                  << std::scientific << std::setprecision(16)
+                  << final_residual_norm_ / rhsnorm << "\n";
+
+        final_residual_norm_ = initial_residual_norm_;
+      }
+    }
+
     return 0;
   }
 
@@ -559,6 +639,11 @@ namespace ReSolve
   {
     vec_V_ = new vector_type(n_, restart_ + 1);
     vec_V_->allocate(memspace_);
+    vec_Y_ = new vector_type(n_);
+    vec_Y_->allocate(memspace_);
+    vec_R_ = new vector_type(n_);
+    vec_R_->allocate(memspace_);
+
     if (flexible_)
     {
       vec_Z_ = new vector_type(n_, restart_ + 1);
@@ -585,13 +670,17 @@ namespace ReSolve
     delete[] h_rs_;
     delete vec_V_;
     delete vec_Z_;
+    delete vec_R_;
+    delete vec_Y_;
 
     h_H_   = nullptr;
     h_c_   = nullptr;
     h_s_   = nullptr;
     h_rs_  = nullptr;
     vec_V_ = nullptr;
     vec_Z_ = nullptr;
+    vec_Y_ = nullptr;
+    vec_R_ = nullptr;
 
     return 0;
   }

resolve/LinSolverIterativeFGMRES.hpp

@@ -97,6 +97,8 @@ namespace ReSolve
 
     vector_type* vec_V_{nullptr};
     vector_type* vec_Z_{nullptr};
+    vector_type* vec_R_{nullptr};
+    vector_type* vec_Y_{nullptr};
 
     real_type* h_H_{nullptr};
     real_type* h_c_{nullptr};