From 50a10483e4045b5cdf01b0763872b8d0bf248512 Mon Sep 17 00:00:00 2001
From: kswirydo <kasia.swirydowicz@gmail.com>
Date: Mon, 11 Dec 2023 21:53:15 -0500
Subject: [PATCH] fixing failing test (#109)

* fixing failing test klu_rf_fgmres_test.exe

---------

Co-authored-by: kswirydo <kasia.swirydowicz@pnnl.gov>
Co-authored-by: pelesh <peless@ornl.gov>
---
 resolve/LinSolverIterativeFGMRES.cpp      | 15 +++++++++++++++
 tests/functionality/testKLU_Rf_FGMRES.cpp | 12 +++++++-----
 2 files changed, 22 insertions(+), 5 deletions(-)

diff --git a/resolve/LinSolverIterativeFGMRES.cpp b/resolve/LinSolverIterativeFGMRES.cpp
index 3585ea78..12e97a16 100644
--- a/resolve/LinSolverIterativeFGMRES.cpp
+++ b/resolve/LinSolverIterativeFGMRES.cpp
@@ -1,6 +1,7 @@
 #include <iostream>
 #include <cassert>
 #include <cmath>
+#include <iomanip>
 
 #include <resolve/utilities/logger/Logger.hpp>
 #include <resolve/matrix/MatrixHandler.hpp>
@@ -112,6 +113,8 @@ namespace ReSolve
   {
     using namespace constants;
 
+    //io::Logger::setVerbosity(io::Logger::EVERYTHING);
+    
     int outer_flag = 1;
     int notconv = 1; 
     int i = 0;
@@ -140,6 +143,9 @@ namespace ReSolve
     //rnorm = ||V_1||
     rnorm = sqrt(rnorm);
     bnorm = sqrt(bnorm);
+    io::Logger::misc() << "it 0: norm of residual "
+                       << std::scientific << std::setprecision(16) 
+                       << rnorm << " Norm of rhs: " << bnorm << "\n";
     initial_residual_norm_ = rnorm;
     while(outer_flag) {
       // check if maybe residual is already small enough?
@@ -227,12 +233,18 @@ namespace ReSolve
 
         // residual norm estimate
         rnorm = fabs(h_rs_[i + 1]);
+        io::Logger::misc() << "it: "<<it<< " --> norm of the residual "
+                           << std::scientific << std::setprecision(16)
+                           << rnorm << "\n";
         // check convergence
         if(i + 1 >= restart_ || rnorm <= tolrel || it >= maxit_) {
           notconv = 0;
         }
       } // inner while
 
+      io::Logger::misc() << "End of cycle, ESTIMATED norm of residual "
+                         << std::scientific << std::setprecision(16)
+                         << rnorm << "\n";
       // solve tri system
       h_rs_[i] = h_rs_[i] / h_H_[i * (restart_ + 1) + i];
       for(int ii = 2; ii <= i + 1; ii++) {
@@ -282,6 +294,9 @@ namespace ReSolve
       if(!outer_flag) {
         final_residual_norm_ = rnorm;
         total_iters_ = it;
+        io::Logger::misc() << "End of cycle, COMPUTED norm of residual "
+                           << std::scientific << std::setprecision(16)
+                           << rnorm << "\n";
       }
     } // outer while
     return 0;
diff --git a/tests/functionality/testKLU_Rf_FGMRES.cpp b/tests/functionality/testKLU_Rf_FGMRES.cpp
index 664ea446..488abfd4 100644
--- a/tests/functionality/testKLU_Rf_FGMRES.cpp
+++ b/tests/functionality/testKLU_Rf_FGMRES.cpp
@@ -39,7 +39,7 @@ int main(int argc, char *argv[])
   ReSolve::LinSolverDirectKLU* KLU = new ReSolve::LinSolverDirectKLU;
 
   ReSolve::LinSolverDirectCuSolverRf* Rf = new ReSolve::LinSolverDirectCuSolverRf;
-  ReSolve::GramSchmidt* GS = new ReSolve::GramSchmidt(vector_handler, ReSolve::GramSchmidt::mgs_pm);
+  ReSolve::GramSchmidt* GS = new ReSolve::GramSchmidt(vector_handler, ReSolve::GramSchmidt::cgs2);
   ReSolve::LinSolverIterativeFGMRES* FGMRES = new ReSolve::LinSolverIterativeFGMRES(matrix_handler, vector_handler, GS);
   // Input to this code is location of `data` directory where matrix files are stored
   const std::string data_path = (argc == 2) ? argv[1] : "./";
@@ -167,8 +167,10 @@ int main(int argc, char *argv[])
   index_type* Q = KLU->getQOrdering();
   error_sum += Rf->setup(A, L, U, P, Q); 
 
-  FGMRES->setMaxit(200); 
+  FGMRES->setMaxit(400); 
   FGMRES->setRestart(100); 
+  FGMRES->setFlexible(1); 
+  FGMRES->setTol(1e-17); 
 
   GS->setup(A->getNumRows(), FGMRES->getRestart()); 
   status =  FGMRES->setup(A); 
@@ -196,7 +198,7 @@ int main(int argc, char *argv[])
 
   A->updateFromCoo(A_coo, ReSolve::memory::DEVICE);
   vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE);
-  Rf->setNumericalProperties(1e-12, 1e-1);
+  Rf->setNumericalProperties(1e-14, 1e-1);
 
   status = Rf->refactorize();
   error_sum += status;
@@ -242,9 +244,9 @@ int main(int argc, char *argv[])
   std::cout<<"\t ||x-x_true||_2              : "<<normDiffMatrix2<<" (solution error)"<<std::endl;
   std::cout<<"\t ||x-x_true||_2/||x_true||_2 : "<<normDiffMatrix2/normXtrue<<" (scaled solution error)"<<std::endl;
   std::cout<<"\t ||b-A*x_exact||_2           : "<<exactSol_normRmatrix2<<" (control; residual norm with exact solution)"<<std::endl;
-  std::cout<<"\t IR iterations               : "<<FGMRES->getNumIter()<<" (max 200, restart 100)"<<std::endl;
+  std::cout<<"\t IR iterations               : "<<FGMRES->getNumIter()<<" (max 400, restart 100)"<<std::endl;
   std::cout<<"\t IR starting res. norm       : "<<FGMRES->getInitResidualNorm()<<" "<<std::endl;
-  std::cout<<"\t IR final res. norm          : "<<FGMRES->getFinalResidualNorm()<<" (tol 1e-14)"<<std::endl<<std::endl;
+  std::cout<<"\t IR final res. norm          : "<<FGMRES->getFinalResidualNorm()<<" "<<std::endl<<std::endl;
 
   if ((normRmatrix1/normB1 > 1e-12 ) || (normRmatrix2/normB2 > 1e-15)) {
     std::cout << "Result inaccurate!\n";