ORNL · shakedregev · Jun 10, 2026 · May 15, 2026 · May 15, 2026 · May 15, 2026
diff --git a/.gitignore b/.gitignore
@@ -2,6 +2,8 @@
 .DS_Store
 build/*
 _build
+LastTest.log
+CTestCostData.txt
 install/*
 *~
 *.swp

diff --git a/CHANGELOG.md b/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

diff --git a/examples/ExampleHelper.hpp b/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_);
           }
         }
 

diff --git a/examples/experimental/r_KLU_GLU_matrix_values_update.cpp b/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);

diff --git a/examples/experimental/r_KLU_cusolverrf_redo_factorization.cpp b/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);

diff --git a/examples/experimental/r_KLU_rf_FGMRES_reuse_factorization.cpp b/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);

diff --git a/examples/experimental/r_KLU_rocsolverrf_redo_factorization.cpp b/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);

diff --git a/examples/gluRefactor.cpp b/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");
 

diff --git a/examples/gpuRefactor.cpp b/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");
 

diff --git a/examples/randGmres.cpp b/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);

diff --git a/examples/sysGmres.cpp b/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();

diff --git a/examples/sysRefactor.cpp b/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);
     }
 

diff --git a/resolve/GramSchmidt.cpp b/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);
 

diff --git a/resolve/LinSolverDirectCuSolverGLU.cpp b/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__);

diff --git a/resolve/LinSolverDirectCuSolverRf.cpp b/resolve/LinSolverDirectCuSolverRf.cpp
@@ -234,7 +234,7 @@ namespace ReSolve
    * @return 0 if successful, 1 otherwise
    */
   int LinSolverDirectCuSolverRf::solve(vector_type* rhs, vector_type* x)
-  {
+  { // CAREFUL: ANYTHING THAT CALLS THIS. DELETE THIS COMMENT
     x->copyFromExternal(rhs->getData(memory::DEVICE), memory::DEVICE, memory::DEVICE);
     x->setDataUpdated(memory::DEVICE);
     status_cusolverrf_ = cusolverRfSolve(handle_cusolverrf_,

diff --git a/resolve/LinSolverDirectRocSolverRf.cpp b/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

diff --git a/resolve/SystemSolver.cpp b/resolve/SystemSolver.cpp
@@ -453,6 +453,8 @@ namespace ReSolve
 
     L_ = factorizationSolver_->getLFactor();
     U_ = factorizationSolver_->getUFactor();
+    L_->allocateMatrixData(memory::DEVICE);
+    U_->allocateMatrixData(memory::DEVICE);
     P_ = factorizationSolver_->getPOrdering();
     Q_ = factorizationSolver_->getQOrdering();
 
@@ -862,6 +864,7 @@ namespace ReSolve
       else
       {
         resVector_->copyFromExternal(rhs, memory::HOST, memory::DEVICE);
+        resVector_->syncData(memory::HOST);
         norm_b = std::sqrt(vectorHandler_->dot(resVector_, resVector_, memory::HOST));
         // ms = memory::HOST;
       }