Require explicit memory management from the user and remove sloppy memory management.

.gitignore

@@ -2,6 +2,8 @@
 .DS_Store
 build/*
 _build
+LastTest.log
+CTestCostData.txt
 install/*
 *~
 *.swp

CHANGELOG.md

@@ -38,6 +38,8 @@
 
 - Added Developer Guide documentation for vector and matrix classes and handlers.
 
+- Remove automatic memory allocation and memory syncing in various memory management functions and require the user to explicitly manage their own memory.
+
 ## Changes to Re::Solve in release 0.99.2
 
 ### Major Features

examples/ExampleHelper.hpp

@@ -130,13 +130,15 @@ namespace ReSolve
         if (res_ == nullptr)
         {
           res_ = new ReSolve::vector::Vector(A->getNumRows());
+          res_->allocate(memspace_);
         }
         else
         {
           if (res_->getSize() != A->getNumRows())
           {
             delete res_;
             res_ = new ReSolve::vector::Vector(A->getNumRows());
+            res_->allocate(memspace_);
           }
         }
 

examples/experimental/r_KLU_GLU_matrix_values_update.cpp

@@ -121,6 +121,7 @@ int main(int argc, char* argv[])
       ReSolve::io::updateArrayFromFile(rhs_file, &rhs);
     }
     // Copy matrix data to device
+    A->allocateMatrixData(ReSolve::memory::DEVICE);
     A->syncData(ReSolve::memory::DEVICE);
 
     std::cout << "Finished reading the matrix and rhs, size: " << A->getNumRows() << " x " << A->getNumColumns()
@@ -133,10 +134,12 @@ int main(int argc, char* argv[])
     // Update host and device data.
     if (i < 1)
     {
+      vec_rhs->allocate(ReSolve::memory::HOST);
       vec_rhs->copyFromExternal(rhs, ReSolve::memory::HOST, ReSolve::memory::HOST);
     }
     else
     {
+      vec_rhs->allocate(ReSolve::memory::DEVICE);
       vec_rhs->copyFromExternal(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE);
     }
     std::cout << "CSR matrix loaded. Expanded NNZ: " << A->getNnz() << std::endl;
@@ -173,6 +176,7 @@ int main(int argc, char* argv[])
       status = GLU->solve(vec_rhs, vec_x);
       std::cout << "CUSOLVER GLU solve status: " << status << std::endl;
     }
+    vec_r->allocate(ReSolve::memory::DEVICE);
     vec_r->copyFromExternal(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE);
 
     matrix_handler->setValuesChanged(true, ReSolve::memory::DEVICE);

examples/experimental/r_KLU_cusolverrf_redo_factorization.cpp

@@ -117,6 +117,7 @@ int main(int argc, char* argv[])
       ReSolve::io::updateArrayFromFile(rhs_file, &rhs);
     }
     // Copy matrix data to device
+    A->allocateMatrixData(ReSolve::memory::DEVICE);
     A->syncData(ReSolve::memory::DEVICE);
 
     std::cout << "Finished reading the matrix and rhs, size: " << A->getNumRows() << " x " << A->getNumColumns()
@@ -129,10 +130,12 @@ int main(int argc, char* argv[])
     // Update host and device data.
     if (i < 2)
     {
+      vec_rhs->allocate(ReSolve::memory::HOST);
       vec_rhs->copyFromExternal(rhs, ReSolve::memory::HOST, ReSolve::memory::HOST);
     }
     else
     {
+      vec_rhs->allocate(ReSolve::memory::DEVICE);
       vec_rhs->copyFromExternal(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE);
     }
     std::cout << "CSR matrix loaded. Expanded NNZ: " << A->getNnz() << std::endl;
@@ -182,6 +185,7 @@ int main(int argc, char* argv[])
     }
 
     // Make sure vec_r is properly initialized before using it.
+    vec_r->allocate(ReSolve::memory::DEVICE);
     vec_r->copyFromExternal(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE);
 
     matrix_handler->setValuesChanged(true, ReSolve::memory::DEVICE);

examples/experimental/r_KLU_rf_FGMRES_reuse_factorization.cpp

@@ -107,13 +107,15 @@ int main(int argc, char* argv[])
       vec_x->allocate(ReSolve::memory::HOST); // for KLU
       vec_x->allocate(ReSolve::memory::DEVICE);
       vec_r = new vector_type(A->getNumRows());
+      vec_r->allocate(ReSolve::memory::DEVICE);
     }
     else
     {
       ReSolve::io::updateMatrixFromFile(mat_file, A);
       ReSolve::io::updateArrayFromFile(rhs_file, &rhs);
     }
     // Copy matrix data to device
+    A->allocateMatrixData(ReSolve::memory::DEVICE);
     A->syncData(ReSolve::memory::DEVICE);
 
     std::cout << "Finished reading the matrix and rhs, size: " << A->getNumRows() << " x " << A->getNumColumns()
@@ -126,11 +128,14 @@ int main(int argc, char* argv[])
     // Update host and device data.
     if (i < 2)
     {
+      vec_rhs->allocate(ReSolve::memory::HOST);
       vec_rhs->copyFromExternal(rhs, ReSolve::memory::HOST, ReSolve::memory::HOST);
       vec_rhs->setDataUpdated(ReSolve::memory::HOST);
     }
     else
     {
+      vec_rhs->allocate(ReSolve::memory::HOST);
+      vec_rhs->allocate(ReSolve::memory::DEVICE);
       A->setUpdated(ReSolve::memory::HOST);
       A->syncData(ReSolve::memory::DEVICE);
       vec_rhs->copyFromExternal(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE);

examples/experimental/r_KLU_rocsolverrf_redo_factorization.cpp

@@ -106,6 +106,7 @@ int main(int argc, char* argv[])
       ReSolve::io::updateArrayFromFile(rhs_file, &rhs);
     }
     // Copy matrix data to device
+    A->allocateMatrixData(ReSolve::memory::DEVICE);
     A->syncData(ReSolve::memory::DEVICE);
 
     std::cout << "Finished reading the matrix and rhs, size: " << A->getNumRows() << " x " << A->getNumColumns()
@@ -118,10 +119,12 @@ int main(int argc, char* argv[])
     // Update host and device data.
     if (i < 2)
     {
+      vec_rhs->allocate(ReSolve::memory::HOST);
       vec_rhs->copyFromExternal(rhs, ReSolve::memory::HOST, ReSolve::memory::HOST);
     }
     else
     {
+      vec_rhs->allocate(ReSolve::memory::DEVICE);
       vec_rhs->copyFromExternal(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE);
     }
     std::cout << "CSR matrix loaded. Expanded NNZ: " << A->getNnz() << std::endl;
@@ -158,6 +161,7 @@ int main(int argc, char* argv[])
     }
 
     // Check accuracy of the solution
+    vec_rhs->allocate(ReSolve::memory::DEVICE);
     vec_r->copyFromExternal(rhs, ReSolve::memory::HOST, ReSolve::memory::DEVICE);
     matrix_handler->setValuesChanged(true, ReSolve::memory::DEVICE);
     matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUS_ONE, ReSolve::memory::DEVICE);

examples/gluRefactor.cpp

@@ -221,7 +221,9 @@ int gluRefactor(int argc, char* argv[])
     rhs_file.close();
 
     // Copy data to device
+    A->allocateMatrixData(memory::DEVICE);
     A->syncData(memory::DEVICE);
+    vec_rhs->allocate(memory::DEVICE);
     vec_rhs->syncData(memory::DEVICE);
     RESOLVE_RANGE_POP("File input");
 

examples/gpuRefactor.cpp

@@ -249,7 +249,9 @@ int gpuRefactor(int argc, char* argv[])
     rhs_file.close();
 
     // Copy data to device
+    A->allocateMatrixData(memory::DEVICE);
     A->syncData(memory::DEVICE);
+    vec_rhs->allocate(memory::DEVICE);
     vec_rhs->syncData(memory::DEVICE);
     RESOLVE_RANGE_POP("File input");
 

examples/randGmres.cpp

@@ -162,8 +162,10 @@ int runGmresExample(int argc, char* argv[])
   vec_x->allocate(memspace);
   if (memspace == memory::DEVICE)
   {
-    A->syncData(memspace);
-    vec_rhs->syncData(memspace);
+    A->allocateMatrixData(memory::DEVICE);
+    A->syncData(memory::DEVICE);
+    vec_rhs->allocate(memory::DEVICE);
+    vec_rhs->syncData(memory::DEVICE);
   }
 
   printSystemInfo(matrix_pathname, A);

examples/sysGmres.cpp

@@ -223,8 +223,10 @@ int sysGmres(int argc, char* argv[])
   if (memspace == memory::DEVICE)
   {
     // Copy data to the device
-    A->syncData(memspace);
-    vec_rhs->syncData(memspace);
+    A->allocateMatrixData(memory::DEVICE);
+    A->syncData(memory::DEVICE);
+    vec_rhs->allocate(memory::DEVICE);
+    vec_rhs->syncData(memory::DEVICE);
   }
 
   mat_file.close();

examples/sysRefactor.cpp

@@ -306,7 +306,9 @@ int sysRefactor(int argc, char* argv[])
     // Ensure matrix data is synced to the device before any GPU operations
     if (hw_backend == "CUDA" || hw_backend == "HIP")
     {
+      A->allocateMatrixData(memory::DEVICE);
       A->syncData(memory::DEVICE);
+      vec_rhs->allocate(memory::DEVICE);
       vec_rhs->syncData(memory::DEVICE);
     }
 
@@ -315,67 +317,71 @@ int sysRefactor(int argc, char* argv[])
 
     double solve_time_ms = 0.0;
 
-    syncDevice(hw_backend);
-    auto solve_start = std::chrono::high_resolution_clock::now();
-
-    // Now call direct solver
-    if (i == 1)
+    if (hw_backend == "CUDA" || hw_backend == "HIP")
+      syncDevice(hw_backend);
     {
-      // Set matrix in solver after the initial matrix is loaded
-      status = solver.setMatrix(A);
-      if (status != 0)
+      auto solve_start = std::chrono::high_resolution_clock::now();
+
+      // Now call direct solver
+      if (i == 1)
       {
-        std::cout << "Failed to set matrix in solver. Status: " << status << std::endl;
-        return 1;
+        // Set matrix in solver after the initial matrix is loaded
+        status = solver.setMatrix(A);
+        if (status != 0)
+        {
+          std::cout << "Failed to set matrix in solver. Status: " << status << std::endl;
+          return 1;
+        }
+
+        // Analysis (symbolic factorization)
+        status = solver.analyze();
+        std::cout << "Analysis on the host status: " << status << std::endl;
+
+        // Numeric factorization on the host
+        status = solver.factorize();
+        std::cout << "Numeric factorization on the host status: " << status << std::endl;
+        // Set up refactorization solver
+        status = solver.refactorizationSetup();
+        std::cout << "Refactorization setup status: " << status << std::endl;
+      }
+      else
+      {
+        // Refactorize on the device
+        status = solver.refactorize();
+        std::cout << "Refactorization on the device status: " << status << std::endl;
       }
 
-      // Analysis (symbolic factorization)
-      status = solver.analyze();
-      std::cout << "Analysis on the host status: " << status << std::endl;
-
-      // Numeric factorization on the host
-      status = solver.factorize();
-      std::cout << "Numeric factorization on the host status: " << status << std::endl;
-      // Set up refactorization solver
-      status = solver.refactorizationSetup();
-      std::cout << "Refactorization setup status: " << status << std::endl;
-    }
-    else
-    {
-      // Refactorize on the device
-      status = solver.refactorize();
-      std::cout << "Refactorization on the device status: " << status << std::endl;
-    }
-
-    status = solver.solve(vec_rhs, vec_x);
-    syncDevice(hw_backend);
-    auto solve_end = std::chrono::high_resolution_clock::now();
-    solve_time_ms  = std::chrono::duration<double, std::milli>(solve_end - solve_start).count();
-    std::cout << "Triangular solve status: " << status << std::endl;
+      status = solver.solve(vec_rhs, vec_x);
+      if (hw_backend == "CUDA" || hw_backend == "HIP")
+        syncDevice(hw_backend);
+      {
+        auto solve_end = std::chrono::high_resolution_clock::now();
+        solve_time_ms  = std::chrono::duration<double, std::milli>(solve_end - solve_start).count();
+        std::cout << "Triangular solve status: " << status << std::endl;
+
+        // Print summary of results
+        helper.printShortSummary(A, vec_rhs, vec_x);
+        if ((i > 1) && is_iterative_refinement)
+        {
+          helper.printIrSummary(&(solver.getIterativeSolver()));
+        }
+
+        if (is_timing)
+        {
+          std::cout << "TIMING,"
+                    << "sysRefactor,"
+                    << hw_backend << ","
+                    << is_iterative_refinement << ","
+                    << i << ","
+                    << solve_time_ms
+                    << std::endl;
+        }
+      }
 
-    // Print summary of results
-    helper.printShortSummary(A, vec_rhs, vec_x);
-    if ((i > 1) && is_iterative_refinement)
-    {
-      helper.printIrSummary(&(solver.getIterativeSolver()));
-    }
+      // Delete objects created on heap
+      delete A;
+      delete vec_x;   // Delete the solution vector
+      delete vec_rhs; // Delete the RHS vector
 
-    if (is_timing)
-    {
-      std::cout << "TIMING,"
-                << "sysRefactor,"
-                << hw_backend << ","
-                << is_iterative_refinement << ","
-                << i << ","
-                << solve_time_ms
-                << std::endl;
+      return 0;
     }
-  }
-
-  // Delete objects created on heap
-  delete A;
-  delete vec_x;   // Delete the solution vector
-  delete vec_rhs; // Delete the RHS vector
-
-  return 0;
-}

resolve/GramSchmidt.cpp

@@ -191,6 +191,10 @@ namespace ReSolve
       // copy H_col to aux, we will need it later
       vec_Hcolumn_->setDataUpdated(memspace_);
       vec_Hcolumn_->resize(i + 1);
+      if (memspace_ == memory::DEVICE)
+      {
+        vec_Hcolumn_->syncData(memory::HOST);
+      }
       vec_Hcolumn_->copyToExternal(h_aux_, 0, memory::HOST, memory::HOST);
 
       // Hcol = V(:,1:i)^T*V(:,i+1);
@@ -201,6 +205,10 @@ namespace ReSolve
 
       // copy H_col to H
       vec_Hcolumn_->setDataUpdated(memspace_);
+      if (memspace_ == memory::DEVICE)
+      {
+        vec_Hcolumn_->syncData(memory::HOST);
+      }
       vec_Hcolumn_->copyToExternal(&H[idxmap(i, 0, num_vecs_ + 1)], 0, memory::HOST, memory::HOST);
 
       // add both pieces together (unstable otherwise, careful here!!)
@@ -379,6 +387,10 @@ namespace ReSolve
 
       // copy H_col to H
       vec_Hcolumn_->setDataUpdated(memspace_);
+      if (memspace_ == memory::DEVICE)
+      {
+        vec_Hcolumn_->syncData(memory::HOST);
+      }
       vec_Hcolumn_->resize(i + 1);
       vec_Hcolumn_->copyToExternal(&H[idxmap(i, 0, num_vecs_ + 1)], 0, memory::HOST, memory::HOST);
 

resolve/LinSolverDirectCuSolverGLU.cpp

@@ -184,6 +184,9 @@ namespace ReSolve
     return error_sum;
   }
 
+  /**
+   * @brief Solve the system. Input vector data must be synced on the device.
+   */
   int LinSolverDirectCuSolverGLU::solve(vector_type* rhs, vector_type* x)
   {
     RESOLVE_RANGE_PUSH(__FUNCTION__);

resolve/LinSolverDirectRocSolverRf.cpp

@@ -77,6 +77,7 @@ namespace ReSolve
     combineFactors(L, U);
 
     M_->setUpdated(ReSolve::memory::HOST);
+    M_->allocateMatrixData(ReSolve::memory::DEVICE);
     M_->syncData(ReSolve::memory::DEVICE);
 
     // remember - P and Q are generally CPU variables