diff --git a/resolve/LinSolverIterativeFGMRES.cpp b/resolve/LinSolverIterativeFGMRES.cpp
index 877d32cc..8a7d28f5 100644
--- a/resolve/LinSolverIterativeFGMRES.cpp
+++ b/resolve/LinSolverIterativeFGMRES.cpp
@@ -15,11 +15,12 @@ namespace ReSolve
     memspace_ = memspace;
     this->matrix_handler_ = nullptr;
     this->vector_handler_ = nullptr;
-    tol_ = 1e-14; //default
-    maxit_= 100; //default
-    restart_ = 10;
-    conv_cond_ = 0;//default
-    flexible_ = 1;
+    // Defaults:
+    // tol_ = 1e-14;
+    // maxit_= 100;
+    // restart_ = 10;
+    // conv_cond_ = 0;
+    // flexible_ = true;
 
     d_V_ = nullptr;
     d_Z_ = nullptr;
@@ -35,11 +36,12 @@ namespace ReSolve
     this->vector_handler_ = vector_handler;
     this->GS_ = gs;
 
-    tol_ = 1e-14; //default
-    maxit_= 100; //default
-    restart_ = 10;
-    conv_cond_ = 0;//default
-    flexible_ = 1;
+    // Defaults:
+    // tol_ = 1e-14;
+    // maxit_= 100;
+    // restart_ = 10;
+    // conv_cond_ = 0;
+    // flexible_ = true;
 
     d_V_ = nullptr;
     d_Z_ = nullptr;
@@ -63,7 +65,7 @@ namespace ReSolve
     maxit_= maxit; 
     restart_ = restart;
     conv_cond_ = conv_cond;
-    flexible_ = 1;
+    flexible_ = true;
 
     d_V_ = nullptr;
     d_Z_ = nullptr;
diff --git a/resolve/SystemSolver.cpp b/resolve/SystemSolver.cpp
index 340114f2..e5389d17 100644
--- a/resolve/SystemSolver.cpp
+++ b/resolve/SystemSolver.cpp
@@ -46,38 +46,44 @@ namespace ReSolve
   }
 
 #ifdef RESOLVE_USE_CUDA
-  SystemSolver::SystemSolver(LinAlgWorkspaceCUDA* workspace, std::string ir) : workspaceCuda_(workspace), irMethod_(ir)
+  SystemSolver::SystemSolver(LinAlgWorkspaceCUDA*  workspaceCuda, 
+                   std::string factor,
+                   std::string refactor,
+                   std::string solve,
+                   std::string ir) 
+    : workspaceCuda_(workspaceCuda),
+      factorizationMethod_(factor),
+      refactorizationMethod_(refactor),
+      solveMethod_(solve),
+      irMethod_(ir)
   {
     // Instantiate handlers
     matrixHandler_ = new MatrixHandler(workspaceCuda_);
     vectorHandler_ = new VectorHandler(workspaceCuda_);
 
-    //set defaults:
     memspace_ = "cuda";
-    factorizationMethod_ = "klu";
-    refactorizationMethod_ = "glu";
-    solveMethod_ = "glu";
-    // irMethod_ = "none";
-    gsMethod_ = "cgs2";
     
     initialize();
   }
 #endif
 
 #ifdef RESOLVE_USE_HIP
-  SystemSolver::SystemSolver(LinAlgWorkspaceHIP* workspace, std::string ir) : workspaceHip_(workspace), irMethod_(ir)
+  SystemSolver::SystemSolver(LinAlgWorkspaceHIP*  workspaceHip, 
+                   std::string factor,
+                   std::string refactor,
+                   std::string solve,
+                   std::string ir) 
+    : workspaceHip_(workspaceHip),
+      factorizationMethod_(factor),
+      refactorizationMethod_(refactor),
+      solveMethod_(solve),
+      irMethod_(ir)
   {
     // Instantiate handlers
     matrixHandler_ = new MatrixHandler(workspaceHip_);
     vectorHandler_ = new VectorHandler(workspaceHip_);
 
-    //set defaults:
     memspace_ = "hip";
-    factorizationMethod_ = "klu";
-    refactorizationMethod_ = "rocsolverrf";
-    solveMethod_ = "rocsolverrf";
-    // irMethod_ = "none";
-    gsMethod_ = "cgs2";
     
     initialize();
   }
@@ -165,7 +171,9 @@ namespace ReSolve
         gs_ = new GramSchmidt(vectorHandler_, GramSchmidt::cgs1);
       } else {
         out::warning() << "Gram-Schmidt variant " << gsMethod_ << " not recognized.\n";
-        gs_ = nullptr;
+        out::warning() << "Using default cgs2 Gram-Schmidt variant.\n";
+        gs_ = new GramSchmidt(vectorHandler_, GramSchmidt::cgs2);
+        gsMethod_ = "cgs2";
       }
 
       iterativeSolver_ = new LinSolverIterativeFGMRES(matrixHandler_,
@@ -354,10 +362,47 @@ namespace ReSolve
     // initialize();
   }
 
-  void SystemSolver::setRefinementMethod(std::string method)
+  void SystemSolver::setRefinementMethod(std::string method, std::string gsMethod)
   {
-    irMethod_ = method;
-    // initialize();
+    if (iterativeSolver_ != nullptr)
+      delete iterativeSolver_;
+
+    if(gs_ != nullptr)
+      delete gs_;
+    
+    if(method == "none")
+      return;
+
+    gsMethod_ = gsMethod;
+
+#if defined(RESOLVE_USE_HIP) || defined(RESOLVE_USE_CUDA)
+    if (method == "fgmres") {
+      if (gsMethod == "cgs2") {
+        gs_ = new GramSchmidt(vectorHandler_, GramSchmidt::cgs2);
+      } else if (gsMethod == "mgs") {
+        gs_ = new GramSchmidt(vectorHandler_, GramSchmidt::mgs);
+      } else if (gsMethod == "mgs_two_synch") {
+        gs_ = new GramSchmidt(vectorHandler_, GramSchmidt::mgs_two_synch);
+      } else if (gsMethod == "mgs_pm") {
+        gs_ = new GramSchmidt(vectorHandler_, GramSchmidt::mgs_pm);
+      } else if (gsMethod == "cgs1") {
+        gs_ = new GramSchmidt(vectorHandler_, GramSchmidt::cgs1);
+      } else {
+        out::warning() << "Gram-Schmidt variant " << gsMethod_ << " not recognized.\n";
+        out::warning() << "Using default cgs2 Gram-Schmidt variant.\n";
+        gs_ = new GramSchmidt(vectorHandler_, GramSchmidt::cgs2);
+        gsMethod_ = "cgs2";
+      }
+
+      iterativeSolver_ = new LinSolverIterativeFGMRES(matrixHandler_,
+                                                      vectorHandler_,
+                                                      gs_,
+                                                      memspace_);
+      irMethod_ = method;
+    } else {
+      out::error() << "Iterative refinement method " << method << " not recognized.\n";
+    }
+#endif
   }
 
   real_type SystemSolver::getResidualNorm(vector_type* rhs, vector_type* x)
diff --git a/resolve/SystemSolver.hpp b/resolve/SystemSolver.hpp
index 4d6a0f61..de6b3d44 100644
--- a/resolve/SystemSolver.hpp
+++ b/resolve/SystemSolver.hpp
@@ -27,9 +27,16 @@ namespace ReSolve
       using matrix_type = matrix::Sparse;
 
       SystemSolver();
-      SystemSolver(LinAlgWorkspaceCUDA* workspaceCuda, std::string ir = "none");
-      SystemSolver(LinAlgWorkspaceHIP*  workspaceHip, std::string ir = "none");
-      SystemSolver(std::string factorizationMethod, std::string refactorizationMethod, std::string solveMethod, std::string IRMethod);
+      SystemSolver(LinAlgWorkspaceCUDA* workspaceCuda, 
+                   std::string factor = "klu",
+                   std::string refactor = "glu",
+                   std::string solve = "glu",
+                   std::string ir = "none");
+      SystemSolver(LinAlgWorkspaceHIP*  workspaceHip, 
+                   std::string factor = "klu",
+                   std::string refactor = "rocsolverrf",
+                   std::string solve = "rocsolverrf",
+                   std::string ir = "none");
 
       ~SystemSolver();
 
@@ -58,12 +65,13 @@ namespace ReSolve
       const std::string getRefactorizationMethod() const;
       const std::string getSolveMethod() const;
       const std::string getRefinementMethod() const;
+      const std::string getOrthogonalizationMethod() const;
 
       // Set solver parameters
       void setFactorizationMethod(std::string method);
       void setRefactorizationMethod(std::string method);
       void setSolveMethod(std::string method);
-      void setRefinementMethod(std::string method);
+      void setRefinementMethod(std::string method, std::string gs = "cgs2");
 
     private:
       LinSolverDirect* factorizationSolver_{nullptr};
diff --git a/tests/functionality/testSysHipRefine.cpp b/tests/functionality/testSysHipRefine.cpp
index 26e0ba89..2c92251d 100644
--- a/tests/functionality/testSysHipRefine.cpp
+++ b/tests/functionality/testSysHipRefine.cpp
@@ -22,7 +22,7 @@ using namespace ReSolve::constants;
 int main(int argc, char *argv[])
 {
   // Use ReSolve data types.
-  using real_type  = ReSolve::real_type;
+  using real_type   = ReSolve::real_type;
   using vector_type = ReSolve::vector::Vector;
 
   //we want error sum to be 0 at the end
@@ -31,14 +31,15 @@ int main(int argc, char *argv[])
   int error_sum = 0;
   int status = 0;
 
-  ReSolve::LinAlgWorkspaceHIP* workspace_HIP = new ReSolve::LinAlgWorkspaceHIP();
-  workspace_HIP->initializeHandles();
-  ReSolve::MatrixHandler* matrix_handler =  new ReSolve::MatrixHandler(workspace_HIP);
-  ReSolve::VectorHandler* vector_handler =  new ReSolve::VectorHandler(workspace_HIP);
+  ReSolve::LinAlgWorkspaceHIP workspace_HIP;
+  workspace_HIP.initializeHandles();
+  ReSolve::MatrixHandler matrix_handler(&workspace_HIP);
+  ReSolve::VectorHandler vector_handler(&workspace_HIP);
 
-  ReSolve::SystemSolver* solver = new ReSolve::SystemSolver(workspace_HIP, "fgmres");
-  solver->getIterativeSolver().setRestart(100);
-  solver->getIterativeSolver().setMaxit(200);
+  ReSolve::SystemSolver solver(&workspace_HIP);
+  solver.setRefinementMethod("fgmres", "cgs2");
+  solver.getIterativeSolver().setRestart(100);
+  solver.getIterativeSolver().setMaxit(200);
 
   // Input to this code is location of `data` directory where matrix files are stored
   const std::string data_path = (argc == 2) ? argv[1] : "./";
@@ -81,21 +82,21 @@ int main(int argc, char *argv[])
   rhs1_file.close();
 
   // Convert first matrix to CSR format
-  matrix_handler->coo2csr(A_coo, A, "cpu");
+  matrix_handler.coo2csr(A_coo, A, "cpu");
   vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::HOST);
   vec_rhs->setDataUpdated(ReSolve::memory::HOST);
 
-  status = solver->setMatrix(A);
+  status = solver.setMatrix(A);
   error_sum += status;
 
   // Solve the first system using KLU
-  status = solver->analyze();
+  status = solver.analyze();
   error_sum += status;
 
-  status = solver->factorize();
+  status = solver.factorize();
   error_sum += status;
 
-  status = solver->solve(vec_rhs, vec_x);
+  status = solver.solve(vec_rhs, vec_x);
   error_sum += status;
 
   vector_type* vec_test;
@@ -113,38 +114,38 @@ int main(int argc, char *argv[])
   vec_r->update(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE);
   vec_diff->update(x_data, ReSolve::memory::HOST, ReSolve::memory::DEVICE);
 
-  // real_type normXmatrix1 = sqrt(vector_handler->dot(vec_test, vec_test, ReSolve::memory::DEVICE));
-  matrix_handler->setValuesChanged(true, "hip");
+  // real_type normXmatrix1 = sqrt(vector_handler.dot(vec_test, vec_test, ReSolve::memory::DEVICE));
+  matrix_handler.setValuesChanged(true, "hip");
   //evaluate the residual ||b-Ax||
-  status = matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE,"csr","hip"); 
+  status = matrix_handler.matvec(A, vec_x, vec_r, &ONE, &MINUSONE,"csr","hip"); 
   error_sum += status;
 
-  real_type normRmatrix1 = sqrt(vector_handler->dot(vec_r, vec_r, "hip"));
+  real_type normRmatrix1 = sqrt(vector_handler.dot(vec_r, vec_r, "hip"));
 
 
   //for testing only - control
 
-  real_type normXtrue = sqrt(vector_handler->dot(vec_x, vec_x, "hip"));
-  real_type normB1 = sqrt(vector_handler->dot(vec_rhs, vec_rhs, "hip"));
+  real_type normXtrue = sqrt(vector_handler.dot(vec_x, vec_x, "hip"));
+  real_type normB1 = sqrt(vector_handler.dot(vec_rhs, vec_rhs, "hip"));
 
   //compute x-x_true
-  vector_handler->axpy(&MINUSONE, vec_x, vec_diff, "hip");
+  vector_handler.axpy(&MINUSONE, vec_x, vec_diff, "hip");
   //evaluate its norm
-  real_type normDiffMatrix1 = sqrt(vector_handler->dot(vec_diff, vec_diff, "hip"));
+  real_type normDiffMatrix1 = sqrt(vector_handler.dot(vec_diff, vec_diff, "hip"));
 
   //compute the residual using exact solution
   vec_r->update(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE);
-  status = matrix_handler->matvec(A, vec_test, vec_r, &ONE, &MINUSONE,"csr", "hip"); 
+  status = matrix_handler.matvec(A, vec_test, vec_r, &ONE, &MINUSONE,"csr", "hip"); 
   error_sum += status;
-  real_type exactSol_normRmatrix1 = sqrt(vector_handler->dot(vec_r, vec_r, "hip"));
+  real_type exactSol_normRmatrix1 = sqrt(vector_handler.dot(vec_r, vec_r, "hip"));
   //evaluate the residual ON THE CPU using COMPUTED solution
 
   vec_r->update(rhs, ReSolve::memory::HOST, ReSolve::memory::HOST);
 
-  status = matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE,"csr", "cpu");
+  status = matrix_handler.matvec(A, vec_x, vec_r, &ONE, &MINUSONE,"csr", "cpu");
   error_sum += status;
 
-  real_type normRmatrix1CPU = sqrt(vector_handler->dot(vec_r, vec_r, "hip"));
+  real_type normRmatrix1CPU = sqrt(vector_handler.dot(vec_r, vec_r, "hip"));
 
   std::cout<<"Results (first matrix): "<<std::endl<<std::endl;
   std::cout<<"\t ||b-A*x||_2                 : " << std::setprecision(16) << normRmatrix1    << " (residual norm)" << std::endl;
@@ -156,7 +157,7 @@ int main(int argc, char *argv[])
 
 
   // Now prepare the Rf solver
-  status = solver->refactorizationSetup();
+  status = solver.refactorizationSetup();
   error_sum += status;
 
   // Load the second matrix
@@ -179,52 +180,52 @@ int main(int argc, char *argv[])
   ReSolve::io::readAndUpdateRhs(rhs2_file, &rhs);
   rhs2_file.close();
 
-  matrix_handler->coo2csr(A_coo, A, "hip");
+  matrix_handler.coo2csr(A_coo, A, "hip");
   vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE);
 
-  status = solver->refactorize();
+  status = solver.refactorize();
   error_sum += status;
   
   vec_x->update(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE);
-  status = solver->solve(vec_rhs, vec_x);
+  status = solver.solve(vec_rhs, vec_x);
   error_sum += status;
   
   vec_rhs->update(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE);
-  status = solver->refine(vec_rhs, vec_x);
+  status = solver.refine(vec_rhs, vec_x);
   error_sum += status;
 
   vec_r->update(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE);
-  matrix_handler->setValuesChanged(true, "hip");
+  matrix_handler.setValuesChanged(true, "hip");
 
   //evaluate final residual
-  status = matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE, "csr", "hip"); 
+  status = matrix_handler.matvec(A, vec_x, vec_r, &ONE, &MINUSONE, "csr", "hip"); 
   error_sum += status;
 
-  real_type normRmatrix2 = sqrt(vector_handler->dot(vec_r, vec_r, "hip"));
+  real_type normRmatrix2 = sqrt(vector_handler.dot(vec_r, vec_r, "hip"));
 
 
   //for testing only - control
-  real_type normB2 = sqrt(vector_handler->dot(vec_rhs, vec_rhs, "hip"));
+  real_type normB2 = sqrt(vector_handler.dot(vec_rhs, vec_rhs, "hip"));
   //compute x-x_true
   vec_diff->update(x_data, ReSolve::memory::HOST, ReSolve::memory::DEVICE);
-  vector_handler->axpy(&MINUSONE, vec_x, vec_diff, "hip");
+  vector_handler.axpy(&MINUSONE, vec_x, vec_diff, "hip");
   //evaluate its norm
-  real_type normDiffMatrix2 = sqrt(vector_handler->dot(vec_diff, vec_diff, "hip"));
+  real_type normDiffMatrix2 = sqrt(vector_handler.dot(vec_diff, vec_diff, "hip"));
 
   //compute the residual using exact solution
   vec_r->update(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE);
-  status = matrix_handler->matvec(A, vec_test, vec_r, &ONE, &MINUSONE, "csr", "hip"); 
+  status = matrix_handler.matvec(A, vec_test, vec_r, &ONE, &MINUSONE, "csr", "hip"); 
   error_sum += status;
-  real_type exactSol_normRmatrix2 = sqrt(vector_handler->dot(vec_r, vec_r, "hip"));
+  real_type exactSol_normRmatrix2 = sqrt(vector_handler.dot(vec_r, vec_r, "hip"));
   std::cout<<"Results (second matrix): "<<std::endl<<std::endl;
   std::cout<<"\t ||b-A*x||_2                 : "<<normRmatrix2<<" (residual norm)"<<std::endl;
   std::cout<<"\t ||b-A*x||_2/||b||_2         : "<<normRmatrix2/normB2<<" (scaled residual norm)"<<std::endl;
   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               : "<<solver->getIterativeSolver().getNumIter()<<" (max 200, restart 100)"<<std::endl;
-  std::cout<<"\t IR starting res. norm       : "<<solver->getIterativeSolver().getInitResidualNorm() <<" "<<std::endl;
-  std::cout<<"\t IR final res. norm          : "<<solver->getIterativeSolver().getFinalResidualNorm() <<" (tol 1e-14)"<<std::endl<<std::endl;
+  std::cout<<"\t IR iterations               : "<<solver.getIterativeSolver().getNumIter()<<" (max 200, restart 100)"<<std::endl;
+  std::cout<<"\t IR starting res. norm       : "<<solver.getIterativeSolver().getInitResidualNorm() <<" "<<std::endl;
+  std::cout<<"\t IR final res. norm          : "<<solver.getIterativeSolver().getFinalResidualNorm() <<" (tol 1e-14)"<<std::endl<<std::endl;
 
   if ((normRmatrix1/normB1 > 1e-12 ) || (normRmatrix2/normB2 > 1e-9)) {
     std::cout << "Result inaccurate!\n";
@@ -237,14 +238,10 @@ int main(int argc, char *argv[])
   }
 
   delete A;
-  delete solver;
   delete [] x;
   delete [] rhs;
   delete vec_r;
   delete vec_x;
-  delete workspace_HIP;
-  delete matrix_handler;
-  delete vector_handler;
 
   return error_sum;
 }