Make CUDA optional dependency

CMakeLists.txt

@@ -27,6 +27,13 @@ option(RESOLVE_USE_GPU  "Use GPU device for computations" ON)
 option(RESOLVE_USE_CUDA "Use CUDA language and SDK" ON)
 set(RESOLVE_CTEST_OUTPUT_DIR ${PROJECT_BINARY_DIR} CACHE PATH "Directory where CTest outputs are saved")
 
+if(RESOLVE_USE_CUDA)
+  set(RESOLVE_USE_GPU On CACHE BOOL "Using GPU!" FORCE)
+else()
+  set(RESOLVE_USE_GPU Off CACHE BOOL "Using GPU!" FORCE)
+endif()
+
+
 set(CMAKE_MACOSX_RPATH 1)
 # set(CMAKE_INSTALL_RPATH "${CMAKE_INSTALL_PREFIX}/lib")
 # https://gitlab.kitware.com/cmake/community/-/wikis/doc/cmake/RPATH-handling#always-full-rpath
@@ -130,10 +137,8 @@ install(FILES "${CMAKE_CURRENT_BINARY_DIR}/ReSolveConfig.cmake"
         "${CMAKE_CURRENT_BINARY_DIR}/ReSolveConfigVersion.cmake"
         DESTINATION share/resolve/cmake)
 
-# Add examples (for now only CUDA examples are available)
-if(RESOLVE_USE_CUDA)
-  add_subdirectory(examples)
-endif(RESOLVE_USE_CUDA)
+# Add usage examples
+add_subdirectory(examples)
 
 # Add tests
 add_subdirectory(tests)

examples/CMakeLists.txt

@@ -14,30 +14,40 @@ target_link_libraries(klu_klu.exe PRIVATE ReSolve)
 add_executable(klu_klu_standalone.exe r_KLU_KLU_standalone.cpp)
 target_link_libraries(klu_klu_standalone.exe PRIVATE ReSolve)
 
-# Build example with KLU factorization and GLU refactorization
-add_executable(klu_glu.exe r_KLU_GLU.cpp)
-target_link_libraries(klu_glu.exe PRIVATE ReSolve)
+# Create CUDA examples
+if(RESOLVE_USE_CUDA)
+
+  # Build example with KLU factorization and GLU refactorization
+  add_executable(klu_glu.exe r_KLU_GLU.cpp)
+  target_link_libraries(klu_glu.exe PRIVATE ReSolve)
 
-# Build example with KLU factorization and Rf refactorization
-add_executable(klu_rf.exe r_KLU_rf.cpp)
-target_link_libraries(klu_rf.exe PRIVATE ReSolve)
+  # Build example with KLU factorization and Rf refactorization
+  add_executable(klu_rf.exe r_KLU_rf.cpp)
+  target_link_libraries(klu_rf.exe PRIVATE ReSolve)
 
-# Build example with KLU factorization, Rf refactorization, and FGMRES iterative refinement
-add_executable(klu_rf_fgmres.exe r_KLU_rf_FGMRES.cpp)
-target_link_libraries(klu_rf_fgmres.exe PRIVATE ReSolve)
+  # Build example with KLU factorization, Rf refactorization, and FGMRES iterative refinement
+  add_executable(klu_rf_fgmres.exe r_KLU_rf_FGMRES.cpp)
+  target_link_libraries(klu_rf_fgmres.exe PRIVATE ReSolve)
 
-# Build example where matrix is factorized once, refactorized once and then the preconditioner is REUSED
-add_executable(klu_rf_fgmres_reuse_refactorization.exe r_KLU_rf_FGMRES_reuse_factorization.cpp)
-target_link_libraries(klu_rf_fgmres_reuse_refactorization.exe PRIVATE ReSolve)
+  # Build example where matrix is factorized once, refactorized once and then the preconditioner is REUSED
+  add_executable(klu_rf_fgmres_reuse_refactorization.exe r_KLU_rf_FGMRES_reuse_factorization.cpp)
+  target_link_libraries(klu_rf_fgmres_reuse_refactorization.exe PRIVATE ReSolve)
 
-# Build example where matrix data is updated 
-add_executable(klu_glu_values_update.exe r_KLU_GLU_matrix_values_update.cpp)
-target_link_libraries(klu_glu_values_update.exe PRIVATE ReSolve)
+  # Build example where matrix data is updated 
+  add_executable(klu_glu_values_update.exe r_KLU_GLU_matrix_values_update.cpp)
+  target_link_libraries(klu_glu_values_update.exe PRIVATE ReSolve)
+
+endif(RESOLVE_USE_CUDA)
 
 
 
 # Install all examples in bin directory
-set(installable_executables klu_klu.exe klu_klu_standalone.exe klu_glu.exe klu_rf.exe klu_rf_fgmres.exe klu_glu_values_update.exe)
+set(installable_executables klu_klu.exe klu_klu_standalone.exe)
+
+if(RESOLVE_USE_CUDA)
+  set(installable_executables ${installable_executables} klu_glu.exe klu_rf.exe klu_rf_fgmres.exe klu_glu_values_update.exe)        
+endif(RESOLVE_USE_CUDA)
+
 install(TARGETS ${installable_executables} 
         RUNTIME DESTINATION bin)
 

examples/r_KLU_GLU.cpp

@@ -9,7 +9,7 @@
 #include <resolve/vector/VectorHandler.hpp>
 #include <resolve/LinSolverDirectKLU.hpp>
 #include <resolve/LinSolverDirectCuSolverGLU.hpp>
-#include <resolve/LinAlgWorkspace.hpp>
+#include <resolve/workspace/LinAlgWorkspace.hpp>
 
 using namespace ReSolve::constants;
 
@@ -145,7 +145,7 @@ int main(int argc, char *argv[])
     vec_r->update(rhs, "cpu", "cuda");
 
 
-    matrix_handler->setValuesChanged(true);
+    matrix_handler->setValuesChanged(true, "cuda");
     matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE,"csr", "cuda"); 
 
     printf("\t 2-Norm of the residual: %16.16e\n", sqrt(vector_handler->dot(vec_r, vec_r, "cuda")));

examples/r_KLU_GLU_matrix_values_update.cpp

@@ -9,7 +9,7 @@
 #include <resolve/vector/VectorHandler.hpp>
 #include <resolve/LinSolverDirectKLU.hpp>
 #include <resolve/LinSolverDirectCuSolverGLU.hpp>
-#include <resolve/LinAlgWorkspace.hpp>
+#include <resolve/workspace/LinAlgWorkspace.hpp>
 
 // this updates the matrix values to simulate what CFD/optimization software does. 
 
@@ -155,7 +155,7 @@ int main(int argc, char *argv[])
       vec_r->update(rhs, "cpu", "cuda");
 
 
-      matrix_handler->setValuesChanged(true);
+      matrix_handler->setValuesChanged(true, "cuda");
       matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE,"csr", "cuda"); 
 
       printf("\t 2-Norm of the residual: %16.16e\n", sqrt(vector_handler->dot(vec_r, vec_r, "cuda")));

examples/r_KLU_KLU.cpp

@@ -1,5 +1,6 @@
 #include <string>
 #include <iostream>
+#include <cmath>
 
 #include <resolve/matrix/Coo.hpp>
 #include <resolve/matrix/Csr.hpp>
@@ -9,7 +10,7 @@
 #include <resolve/matrix/MatrixHandler.hpp>
 #include <resolve/vector/VectorHandler.hpp>
 #include <resolve/LinSolverDirectKLU.hpp>
-#include <resolve/LinAlgWorkspace.hpp>
+#include <resolve/workspace/LinAlgWorkspace.hpp>
 
 using namespace ReSolve::constants;
 
@@ -35,10 +36,9 @@ int main(int argc, char *argv[])
 
   ReSolve::matrix::Coo* A_coo;
   ReSolve::matrix::Csr* A;
-  ReSolve::LinAlgWorkspaceCUDA* workspace_CUDA = new ReSolve::LinAlgWorkspaceCUDA;
-  workspace_CUDA->initializeHandles();
-  ReSolve::MatrixHandler* matrix_handler =  new ReSolve::MatrixHandler(workspace_CUDA);
-  ReSolve::VectorHandler* vector_handler =  new ReSolve::VectorHandler(workspace_CUDA);
+  ReSolve::LinAlgWorkspaceCpu* workspace = new ReSolve::LinAlgWorkspaceCpu();
+  ReSolve::MatrixHandler* matrix_handler = new ReSolve::MatrixHandler(workspace);
+  ReSolve::VectorHandler* vector_handler = new ReSolve::VectorHandler(workspace);
   real_type* rhs;
   real_type* x;
 
@@ -96,7 +96,11 @@ int main(int argc, char *argv[])
       ReSolve::io::readAndUpdateMatrix(mat_file, A_coo);
       ReSolve::io::readAndUpdateRhs(rhs_file, &rhs);
     }
-    std::cout<<"Finished reading the matrix and rhs, size: "<<A->getNumRows()<<" x "<<A->getNumColumns()<< ", nnz: "<< A->getNnz()<< ", symmetric? "<<A->symmetric()<< ", Expanded? "<<A->expanded()<<std::endl;
+    std::cout << "Finished reading the matrix and rhs, size: " << A->getNumRows()
+              << " x "           << A->getNumColumns()
+              << ", nnz: "       << A->getNnz()
+              << ", symmetric? " << A->symmetric()
+              << ", Expanded? "  << A->expanded() << std::endl;
     mat_file.close();
     rhs_file.close();
 
@@ -106,8 +110,8 @@ int main(int argc, char *argv[])
       vec_rhs->update(rhs, "cpu", "cpu");
       vec_rhs->setDataUpdated("cpu");
     } else { 
-      matrix_handler->coo2csr(A_coo, A, "cuda");
-      vec_rhs->update(rhs, "cpu", "cuda");
+      matrix_handler->coo2csr(A_coo, A, "cpu");
+      vec_rhs->update(rhs, "cpu", "cpu");
     }
     std::cout<<"COO to CSR completed. Expanded NNZ: "<< A->getNnzExpanded()<<std::endl;
     //Now call direct solver
@@ -129,15 +133,15 @@ int main(int argc, char *argv[])
       status = KLU->solve(vec_rhs, vec_x);
       std::cout<<"KLU solve status: "<<status<<std::endl;      
     }
-    vec_r->update(rhs, "cpu", "cuda");
+    vec_r->update(rhs, "cpu", "cpu");
 
-    matrix_handler->setValuesChanged(true);
+    matrix_handler->setValuesChanged(true, "cpu");
 
-    matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE, "csr", "cuda"); 
+    matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE, "csr", "cpu"); 
     real_type* test = vec_r->getData("cpu");
     (void) test; // TODO: Do we need `test` variable in this example?
 
-    printf("\t 2-Norm of the residual: %16.16e\n", sqrt(vector_handler->dot(vec_r, vec_r, "cuda")));
+    printf("\t 2-Norm of the residual: %16.16e\n", sqrt(vector_handler->dot(vec_r, vec_r, "cpu")));
 
 
   }

examples/r_KLU_KLU_standalone.cpp

@@ -1,5 +1,6 @@
 #include <string>
 #include <iostream>
+#include <cmath>
 
 #include <resolve/matrix/Coo.hpp>
 #include <resolve/matrix/Csr.hpp>
@@ -9,7 +10,7 @@
 #include <resolve/matrix/MatrixHandler.hpp>
 #include <resolve/vector/VectorHandler.hpp>
 #include <resolve/LinSolverDirectKLU.hpp>
-#include <resolve/LinAlgWorkspace.hpp>
+#include <resolve/workspace/LinAlgWorkspace.hpp>
 
 using namespace ReSolve::constants;
 
@@ -31,10 +32,9 @@ int main(int argc, char *argv[])
 
   ReSolve::matrix::Coo* A_coo;
   ReSolve::matrix::Csr* A;
-  ReSolve::LinAlgWorkspaceCUDA* workspace_CUDA = new ReSolve::LinAlgWorkspaceCUDA;
-  workspace_CUDA->initializeHandles();
-  ReSolve::MatrixHandler* matrix_handler =  new ReSolve::MatrixHandler(workspace_CUDA);
-  ReSolve::VectorHandler* vector_handler =  new ReSolve::VectorHandler(workspace_CUDA);
+  ReSolve::LinAlgWorkspaceCpu* workspace = new ReSolve::LinAlgWorkspaceCpu();
+  ReSolve::MatrixHandler* matrix_handler = new ReSolve::MatrixHandler(workspace);
+  ReSolve::VectorHandler* vector_handler = new ReSolve::VectorHandler(workspace);
   real_type* rhs;
   real_type* x;
 
@@ -95,15 +95,15 @@ int main(int argc, char *argv[])
   std::cout << "KLU factorization status: " << status << std::endl;
   status = KLU->solve(vec_rhs, vec_x);
   std::cout << "KLU solve status: " << status << std::endl;      
-  vec_r->update(rhs, "cpu", "cuda");
+  vec_r->update(rhs, "cpu", "cpu");
 
-  matrix_handler->setValuesChanged(true);
+  matrix_handler->setValuesChanged(true, "cpu");
 
-  matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE, "csr", "cuda"); 
+  matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE, "csr", "cpu"); 
   real_type* test = vec_r->getData("cpu");
   (void) test; // TODO: Do we need `test` variable in this example?
 
-  printf("\t 2-Norm of the residual: %16.16e\n", sqrt(vector_handler->dot(vec_r, vec_r, "cuda")));
+  printf("\t 2-Norm of the residual: %16.16e\n", sqrt(vector_handler->dot(vec_r, vec_r, "cpu")));
 
 
 

examples/r_KLU_rf.cpp

@@ -10,7 +10,7 @@
 #include <resolve/vector/VectorHandler.hpp>
 #include <resolve/LinSolverDirectKLU.hpp>
 #include <resolve/LinSolverDirectCuSolverRf.hpp>
-#include <resolve/LinAlgWorkspace.hpp>
+#include <resolve/workspace/LinAlgWorkspace.hpp>
 
 using namespace ReSolve::constants;
 
@@ -158,7 +158,7 @@ int main(int argc, char *argv[] )
     }
     vec_r->update(rhs, "cpu", "cuda");
 
-    matrix_handler->setValuesChanged(true);
+    matrix_handler->setValuesChanged(true, "cuda");
 
     matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE,"csr", "cuda"); 
 

examples/r_KLU_rf_FGMRES.cpp

@@ -11,7 +11,7 @@
 #include <resolve/LinSolverDirectKLU.hpp>
 #include <resolve/LinSolverDirectCuSolverRf.hpp>
 #include <resolve/LinSolverIterativeFGMRES.hpp>
-#include <resolve/LinAlgWorkspace.hpp>
+#include <resolve/workspace/LinAlgWorkspace.hpp>
 
 using namespace ReSolve::constants;
 
@@ -125,7 +125,7 @@ int main(int argc, char *argv[])
     real_type norm_b;
     if (i < 2){
       KLU->setup(A);
-      matrix_handler->setValuesChanged(true);
+      matrix_handler->setValuesChanged(true, "cuda");
       status = KLU->analyze();
       std::cout<<"KLU analysis status: "<<status<<std::endl;
       status = KLU->factorize();
@@ -135,7 +135,7 @@ int main(int argc, char *argv[])
       vec_r->update(rhs, "cpu", "cuda");
       norm_b = vector_handler->dot(vec_r, vec_r, "cuda");
       norm_b = sqrt(norm_b);
-      matrix_handler->setValuesChanged(true);
+      matrix_handler->setValuesChanged(true, "cuda");
       matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE,"csr", "cuda"); 
       printf("\t 2-Norm of the residual : %16.16e\n", sqrt(vector_handler->dot(vec_r, vec_r, "cuda"))/norm_b);
       if (i == 1) {
@@ -165,7 +165,7 @@ int main(int argc, char *argv[])
        norm_b = vector_handler->dot(vec_r, vec_r, "cuda");
       norm_b = sqrt(norm_b);
 
-      //matrix_handler->setValuesChanged(true);
+      //matrix_handler->setValuesChanged(true, "cuda");
       FGMRES->resetMatrix(A);
       FGMRES->setupPreconditioner("CuSolverRf", Rf);
 

examples/r_KLU_rf_FGMRES_reuse_factorization.cpp

@@ -11,7 +11,7 @@
 #include <resolve/LinSolverDirectKLU.hpp>
 #include <resolve/LinSolverDirectCuSolverRf.hpp>
 #include <resolve/LinSolverIterativeFGMRES.hpp>
-#include <resolve/LinAlgWorkspace.hpp>
+#include <resolve/workspace/LinAlgWorkspace.hpp>
 
 using namespace ReSolve::constants;
 
@@ -126,7 +126,7 @@ int main(int argc, char *argv[])
     real_type norm_b;
     if (i < 2){
       KLU->setup(A);
-      matrix_handler->setValuesChanged(true);
+      matrix_handler->setValuesChanged(true, "cuda");
       status = KLU->analyze();
       std::cout<<"KLU analysis status: "<<status<<std::endl;
       status = KLU->factorize();
@@ -136,7 +136,7 @@ int main(int argc, char *argv[])
       vec_r->update(rhs, "cpu", "cuda");
       norm_b = vector_handler->dot(vec_r, vec_r, "cuda");
       norm_b = sqrt(norm_b);
-      matrix_handler->setValuesChanged(true);
+      matrix_handler->setValuesChanged(true, "cuda");
       matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE,"csr", "cuda"); 
       printf("\t 2-Norm of the residual : %16.16e\n", sqrt(vector_handler->dot(vec_r, vec_r, "cuda"))/norm_b);
       if (i == 1) {
@@ -178,7 +178,7 @@ int main(int argc, char *argv[])
        norm_b = vector_handler->dot(vec_r, vec_r, "cuda");
       norm_b = sqrt(norm_b);
 
-      matrix_handler->setValuesChanged(true);
+      matrix_handler->setValuesChanged(true, "cuda");
       FGMRES->resetMatrix(A);
 
       matrix_handler->matvec(A, vec_x, vec_r, &ONE, &MINUSONE,"csr", "cuda");