Disallow sloppy calls to functions that allocate memory when they shouldn't.

examples/ExampleHelper.hpp

@@ -142,14 +142,10 @@ namespace ReSolve
                        ReSolve::vector::Vector* r,
                        ReSolve::vector::Vector* x)
       {
+        assert(res_ != nullptr && "resetSystem should be called after setSystem");
         A_ = A;
         r_ = r;
         x_ = x;
-        if (res_ == nullptr)
-        {
-          res_ = new ReSolve::vector::Vector(A->getNumRows());
-        }
-
         computeNorms();
       }
 

examples/gluRefactor.cpp

@@ -264,7 +264,14 @@ int gluRefactor(int argc, char* argv[])
     RESOLVE_RANGE_POP("Triangular solve");
 
     // Print summary of the results
-    helper.resetSystem(A, vec_rhs, vec_x);
+    if (i == 0)
+    {
+      helper.setSystem(A, vec_rhs, vec_x);
+    }
+    else
+    {
+      helper.resetSystem(A, vec_rhs, vec_x);
+    }
     helper.printSummary();
 
   } // for (int i = 0; i < num_systems; ++i)

examples/gpuRefactor.cpp

@@ -260,7 +260,7 @@ int gpuRefactor(int argc, char* argv[])
       std::cout << "KLU solve status: " << status << std::endl;
 
       // Print summary of results
-      helper.resetSystem(A, vec_rhs, vec_x);
+      helper.setSystem(A, vec_rhs, vec_x);
       helper.printShortSummary();
 
       if (i == 1)

examples/kluFactor.cpp

@@ -184,8 +184,14 @@ int main(int argc, char* argv[])
 
     status = KLU.solve(vec_rhs, vec_x);
     std::cout << "KLU solve status: " << status << std::endl;
-
-    helper.resetSystem(A, vec_rhs, vec_x);
+    if (i == 0)
+    {
+      helper.setSystem(A, vec_rhs, vec_x);
+    }
+    else
+    {
+      helper.resetSystem(A, vec_rhs, vec_x);
+    }
     helper.printShortSummary();
     if (is_iterative_refinement)
     {

examples/kluRefactor.cpp

@@ -191,7 +191,14 @@ int main(int argc, char* argv[])
     status = KLU->solve(vec_rhs, vec_x);
     std::cout << "KLU solve status: " << status << std::endl;
 
-    helper.resetSystem(A, vec_rhs, vec_x);
+    if (i == 0)
+    {
+      helper.setSystem(A, vec_rhs, vec_x);
+    }
+    else
+    {
+      helper.resetSystem(A, vec_rhs, vec_x);
+    }
     helper.printShortSummary();
     if (is_iterative_refinement)
     {

examples/randGmres.cpp

@@ -183,7 +183,7 @@ int runGmresExample(int argc, char* argv[])
   FGMRES.solve(vec_rhs, vec_x);
 
   // Print summary of results
-  helper.resetSystem(A, vec_rhs, vec_x);
+  helper.setSystem(A, vec_rhs, vec_x);
   std::cout << "\nRandomized GMRES result on " << hwbackend << "\n";
   std::cout << "---------------------------------\n";
   helper.printIrSummary(&FGMRES);

examples/sysRefactor.cpp

@@ -320,7 +320,14 @@ int sysRefactor(int argc, char* argv[])
     std::cout << "Triangular solve status: " << status << std::endl;
 
     // Print summary of results
-    helper.resetSystem(A, vec_rhs, vec_x);
+    if (i == 0)
+    {
+      helper.setSystem(A, vec_rhs, vec_x);
+    }
+    else
+    {
+      helper.resetSystem(A, vec_rhs, vec_x);
+    }
     helper.printShortSummary();
     if ((i > 1) && is_iterative_refinement)
     {

resolve/LinSolverDirect.cpp

@@ -55,6 +55,7 @@ namespace ReSolve
   {
     if (A == nullptr)
     {
+      out::error() << "LinSolverDirect::setup: input matrix A is nullptr" << std::endl;
       return 1;
     }
     A_ = A;

resolve/LinSolverDirectCuSolverRf.cpp

@@ -50,6 +50,7 @@ namespace ReSolve
    *
    * Sets up the cuSolverRf factorization for the given matrix A and its
    * L and U factors. The permutation vectors P and Q are also set up.
+   * This function should not be called more than once for the same object.
    *
    * @param[in] A - pointer to the matrix A
    * @param[in] L - pointer to the lower triangular factor L in CSR
@@ -75,31 +76,41 @@ namespace ReSolve
     this->A_      = A;
     index_type n  = A_->getNumRows();
 
-    // Remember - P and Q are generally CPU variables!
-    // Factorization data is stored in the handle.
-    // If function is called again, destroy the old handle to get rid of old data.
     if (setup_completed_)
     {
-      cusolverRfDestroy(handle_cusolverrf_);
-      cusolverRfCreate(&handle_cusolverrf_);
+      out::error() << "Trying to setup LinSolverDirectCuSolverRf, but the setup has been already done! " << std::endl;
+      return 1;
     }
 
     if (d_P_ == nullptr)
     {
       mem_.allocateArrayOnDevice(&d_P_, n);
     }
+    else
+    {
+      out::error() << "Trying to allocate permutation vector P in " << __func__ << " in LinSolverDirectCuSolverRf, but the permutation vector P is already allocated! " << std::endl;
+      return 1;
+    }
 
     if (d_Q_ == nullptr)
     {
       mem_.allocateArrayOnDevice(&d_Q_, n);
     }
-
-    if (d_T_ != nullptr)
+    else
     {
-      mem_.deleteOnDevice(d_T_);
+      out::error() << "Trying to allocate permutation vector Q in " << __func__ << " in LinSolverDirectCuSolverRf, but the permutation vector Q is already allocated! " << std::endl;
+      return 1;
     }
 
-    mem_.allocateArrayOnDevice(&d_T_, n);
+    if (d_T_ == nullptr)
+    {
+      mem_.allocateArrayOnDevice(&d_T_, n);
+    }
+    else
+    {
+      out::error() << "Trying to allocate temporary vector T in " << __func__ << " in LinSolverDirectCuSolverRf, but the temporary vector T is already allocated! " << std::endl;
+      return 1;
+    }
 
     mem_.copyArrayHostToDevice(d_P_, P, n);
     mem_.copyArrayHostToDevice(d_Q_, Q, n);

resolve/LinSolverDirectKLU.cpp

@@ -127,7 +127,7 @@ namespace ReSolve
     if (Symbolic_ == nullptr)
     {
       out::error() << "Symbolic_ factorization failed with Common_.status = "
-                   << Common_.status << "\n";
+                   << Common_.status << std::endl;
       return 1;
     }
     return 0;
@@ -167,6 +167,8 @@ namespace ReSolve
 
     if (Numeric_ == nullptr)
     {
+      out::error() << "Numeric_ factorization failed with Common_.status = "
+                   << Common_.status << std::endl;
       return 1;
     }
     else

resolve/LinSolverDirectRocSolverRf.cpp

@@ -82,11 +82,21 @@ namespace ReSolve
     {
       mem_.allocateArrayOnDevice(&d_P_, n);
     }
+    else
+    {
+      out::error() << "Trying to allocate permutation vector P in " << __func__ << " in LinSolverDirectRocSolverRf, but the permutation vector P is already allocated! " << std::endl;
+      return 1;
+    }
 
     if (d_Q_ == nullptr)
     {
       mem_.allocateArrayOnDevice(&d_Q_, n);
     }
+    else
+    {
+      out::error() << "Trying to allocate permutation vector Q in " << __func__ << " in LinSolverDirectRocSolverRf, but the permutation vector Q is already allocated! " << std::endl;
+      return 1;
+    }
     mem_.copyArrayHostToDevice(d_P_, P, n);
     mem_.copyArrayHostToDevice(d_Q_, Q, n);
 

resolve/LinSolverIterative.cpp

@@ -26,6 +26,7 @@ namespace ReSolve
   {
     if (A == nullptr)
     {
+      out::error() << "LinSolverIterative::setup: input matrix A is nullptr" << std::endl;
       return 1;
     }
     this->A_ = A;

resolve/MemoryUtils.hpp

@@ -3,9 +3,12 @@
 #include <cstring> // <- declares `memcpy`
 
 #include <resolve/resolve_defs.hpp>
+#include <resolve/utilities/logger/Logger.hpp>
 
 namespace ReSolve
 {
+  using out = ReSolve::io::Logger;
+
   namespace memory
   {
     enum MemorySpace
@@ -86,12 +89,22 @@ namespace ReSolve
     {
       std::size_t arraysize = static_cast<std::size_t>(n) * sizeof(T);
       *v                    = new T[arraysize];
-      return *v == nullptr ? 1 : 0;
+      if (*v == nullptr)
+      {
+        out::error() << "Memory allocation on host failed for size " << arraysize << " bytes." << std::endl;
+        return 1;
+      }
+      return 0;
     }
 
     template <typename T>
     int deleteOnHost(T* v)
     {
+      if (v == nullptr)
+      {
+        out::error() << "Trying to delete nullptr on host!" << std::endl;
+        return 1;
+      }
       delete[] v;
       v = nullptr;
       return 0;

resolve/SystemSolver.cpp

@@ -372,10 +372,6 @@ namespace ReSolve
                                                       vectorHandler_,
                                                       gs_);
     }
-    else
-    {
-      // do nothing
-    }
 
     return 0;
   }
@@ -393,32 +389,46 @@ namespace ReSolve
       factorizationSolver_->setup(A_);
       return factorizationSolver_->analyze();
     }
+    out::error() << "Analysis method " << factorizationSolver_ << " not recognized!\n"
+                 << "Available options: klu.\n";
     return 1;
   }
 
   int SystemSolver::factorize()
   {
+    if (A_ == nullptr)
+    {
+      out::error() << "System matrix not set!\n";
+      return 1;
+    }
     if (factorizationMethod_ == "klu")
     {
       is_solve_on_device_ = false;
       return factorizationSolver_->factorize();
     }
+    out::error() << "Factorization method << factorizationSolver_ << not recognized!\n"
+                 << "Available options: klu.\n";
     return 1;
   }
 
   int SystemSolver::refactorize()
   {
+    if (A_ == nullptr)
+    {
+      out::error() << "System matrix not set!\n";
+      return 1;
+    }
     if (refactorizationMethod_ == "klu")
     {
       return factorizationSolver_->refactorize();
     }
-
     if (refactorizationMethod_ == "glu" || refactorizationMethod_ == "cusolverrf" || refactorizationMethod_ == "rocsolverrf")
     {
       is_solve_on_device_ = true;
       return refactorizationSolver_->refactorize();
     }
-
+    out::error() << "Refactorization method " << refactorizationSolver_ << " not recognized!\n"
+                 << "Available options: klu, glu, cusolverrf, rocsolverrf.\n";
     return 1;
   }
 

resolve/matrix/io.cpp

@@ -383,10 +383,14 @@ namespace ReSolve
     {
       if (!file)
       {
-        Logger::error() << "Empty input to updateArrayFromFile function ..." << std::endl;
+        Logger::error() << "Empty input to updateArrayFromFile function." << std::endl;
+        return;
+      }
+      if (p_rhs == nullptr)
+      {
+        Logger::error() << "Memory not allocated for updateArrayFromFile function." << std::endl;
         return;
       }
-
       real_type*        rhs = *p_rhs;
       std::stringstream ss;
       std::string       line;
@@ -399,11 +403,6 @@ namespace ReSolve
       }
       ss << line;
       ss >> n >> m;
-
-      if (rhs == nullptr)
-      {
-        rhs = new real_type[n];
-      }
       real_type  a;
       index_type i = 0;
       while (file >> a)