add a LUSOL functionality test

.github/workflows/spack_cpu_build.yaml

@@ -82,6 +82,7 @@ jobs:
  # Minimal Build(s) - GHCR mirror speeds these up a lot!
  spack_spec:
  - resolve@develop~klu~lusol
+ - resolve@develop~klu+lusol
  - resolve@develop+klu~lusol
  - resolve@develop+klu+lusol
 

resolve/matrix/CMakeLists.txt

@@ -15,6 +15,7 @@ set(Matrix_SRC
  Coo.cpp
  MatrixHandler.cpp
  MatrixHandlerCpu.cpp
+ Utilities.cpp
 )
 
 # C++ code that depends on CUDA SDK libraries
@@ -35,6 +36,7 @@ set(Matrix_HEADER_INSTALL
  Csr.hpp
  Csc.hpp
  MatrixHandler.hpp
+ Utilities.hpp
 )
 
 # Build shared library ReSolve::matrix

resolve/matrix/Coo.cpp

@@ -1,12 +1,12 @@
 #include <cstring> // <-- includes memcpy
 #include <iostream>
-#include <iomanip> 
+#include <iomanip>
 
 #include <resolve/utilities/logger/Logger.hpp>
 #include "Coo.hpp"
 
 
-namespace ReSolve 
+namespace ReSolve
 {
  using out = io::Logger;
 
@@ -17,15 +17,123 @@ namespace ReSolve
  matrix::Coo::Coo(index_type n, index_type m, index_type nnz) : Sparse(n, m, nnz)
  {
  }
- 
- matrix::Coo::Coo(index_type n, 
- index_type m, 
+
+ matrix::Coo::Coo(index_type n,
+ index_type m,
  index_type nnz,
  bool symmetric,
  bool expanded) : Sparse(n, m, nnz, symmetric, expanded)
  {
  }
-
+
+ /**
+ * @brief Hijacking constructor
+ */
+ matrix::Coo::Coo(index_type n,
+ index_type m,
+ index_type nnz,
+ bool symmetric,
+ bool expanded,
+ index_type** rows,
+ index_type** cols,
+ real_type** vals,
+ memory::MemorySpace memspaceSrc,
+ memory::MemorySpace memspaceDst)
+ : Sparse(n, m, nnz, symmetric, expanded)
+ {
+ using namespace memory;
+ int control = -1;
+ if ((memspaceSrc == memory::HOST) && (memspaceDst == memory::HOST)) { control = 0;}
+ if ((memspaceSrc == memory::HOST) && (memspaceDst == memory::DEVICE)){ control = 1;}
+ if ((memspaceSrc == memory::DEVICE) && (memspaceDst == memory::HOST)) { control = 2;}
+ if ((memspaceSrc == memory::DEVICE) && (memspaceDst == memory::DEVICE)){ control = 3;}
+
+ switch (control)
+ {
+ case 0: // cpu->cpu
+ // Set host data
+ h_row_data_ = *rows;
+ h_col_data_ = *cols;
+ h_val_data_ = *vals;
+ h_data_updated_ = true;
+ owns_cpu_vals_ = true;
+ owns_cpu_data_ = true;
+ // Set device data to null
+ if (d_row_data_ || d_col_data_ || d_val_data_) {
+ out::error() << "Device data unexpectedly allocated. "
+ << "Possible bug in matrix::Sparse class.\n";
+ }
+ d_row_data_ = nullptr;
+ d_col_data_ = nullptr;
+ d_val_data_ = nullptr;
+ d_data_updated_ = false;
+ owns_gpu_vals_ = true;
+ owns_gpu_data_ = false;
+ // Hijack data from the source
+ *rows = nullptr;
+ *cols = nullptr;
+ *vals = nullptr;
+ break;
+ case 2: // gpu->cpu
+ // Set device data and copy it to host
+ d_row_data_ = *rows;
+ d_col_data_ = *cols;
+ d_val_data_ = *vals;
+ d_data_updated_ = true;
+ owns_gpu_vals_ = true;
+ owns_gpu_data_ = true;
+ copyData(memspaceDst);
+ // Hijack data from the source
+ *rows = nullptr;
+ *cols = nullptr;
+ *vals = nullptr;
+ break;
+ case 1: // cpu->gpu
+ // Set host data and copy it to device
+ h_row_data_ = *rows;
+ h_col_data_ = *cols;
+ h_val_data_ = *vals;
+ h_data_updated_ = true;
+ owns_cpu_vals_ = true;
+ owns_cpu_data_ = true;
+ copyData(memspaceDst);
+
+ // Hijack data from the source
+ *rows = nullptr;
+ *cols = nullptr;
+ *vals = nullptr;
+ break;
+ case 3: // gpu->gpu
+ // Set device data
+ d_row_data_ = *rows;
+ d_col_data_ = *cols;
+ d_val_data_ = *vals;
+ d_data_updated_ = true;
+ owns_gpu_vals_ = true;
+ owns_gpu_data_ = true;
+ // Set host data to null
+ if (h_row_data_ || h_col_data_ || h_val_data_) {
+ out::error() << "Host data unexpectedly allocated. "
+ << "Possible bug in matrix::Sparse class.\n";
+ }
+ h_row_data_ = nullptr;
+ h_col_data_ = nullptr;
+ h_val_data_ = nullptr;
+ h_data_updated_ = false;
+ owns_cpu_vals_ = true;
+ owns_cpu_data_ = false;
+ // Hijack data from the source
+ *rows = nullptr;
+ *cols = nullptr;
+ *vals = nullptr;
+ break;
+ default:
+ out::error() << "Coo constructor failed! "
+ << "Possible bug in memory spaces setting.\n";
+ break;
+ }
+ }
+
  matrix::Coo::~Coo()
  {
  }
@@ -90,7 +198,7 @@ namespace ReSolve
  if (((memspaceIn == memory::DEVICE)) && ((memspaceOut == memory::DEVICE))){ control = 3;}
 
  if (memspaceOut == memory::HOST) {
- //check if cpu data allocated 
+ //check if cpu data allocated
  if ((h_row_data_ == nullptr) != (h_col_data_ == nullptr)) {
  out::error() << "In Coo::updateData one of host row or column data is null!\n";
  }
@@ -150,15 +258,15 @@ namespace ReSolve
  return -1;
  }
  return 0;
- } 
+ }
 
  int matrix::Coo::updateData(index_type* row_data, index_type* col_data, real_type* val_data, index_type new_nnz, memory::MemorySpace memspaceIn, memory::MemorySpace memspaceOut)
  {
  this->destroyMatrixData(memspaceOut);
  this->nnz_ = new_nnz;
  int i = this->updateData(row_data, col_data, val_data, memspaceIn, memspaceOut);
  return i;
- } 
+ }
 
  int matrix::Coo::allocateMatrixData(memory::MemorySpace memspace)
  {
@@ -168,20 +276,20 @@ namespace ReSolve
 
  if (memspace == memory::HOST) {
  this->h_row_data_ = new index_type[nnz_current];
- std::fill(h_row_data_, h_row_data_ + nnz_current, 0); 
+ std::fill(h_row_data_, h_row_data_ + nnz_current, 0);
  this->h_col_data_ = new index_type[nnz_current];
- std::fill(h_col_data_, h_col_data_ + nnz_current, 0); 
+ std::fill(h_col_data_, h_col_data_ + nnz_current, 0);
  this->h_val_data_ = new real_type[nnz_current];
- std::fill(h_val_data_, h_val_data_ + nnz_current, 0.0); 
+ std::fill(h_val_data_, h_val_data_ + nnz_current, 0.0);
  owns_cpu_data_ = true;
  owns_cpu_vals_ = true;
  return 0;
  }
 
  if (memspace == memory::DEVICE) {
- mem_.allocateArrayOnDevice(&d_row_data_, nnz_current); 
- mem_.allocateArrayOnDevice(&d_col_data_, nnz_current); 
- mem_.allocateArrayOnDevice(&d_val_data_, nnz_current); 
+ mem_.allocateArrayOnDevice(&d_row_data_, nnz_current);
+ mem_.allocateArrayOnDevice(&d_col_data_, nnz_current);
+ mem_.allocateArrayOnDevice(&d_val_data_, nnz_current);
  owns_gpu_data_ = true;
  owns_gpu_vals_ = true;
  return 0;
@@ -205,12 +313,12 @@ namespace ReSolve
  out::error() << "In Coo::copyData one of host row or column data is null!\n";
  }
  if ((h_row_data_ == nullptr) && (h_col_data_ == nullptr)) {
- h_row_data_ = new index_type[nnz_current]; 
- h_col_data_ = new index_type[nnz_current]; 
+ h_row_data_ = new index_type[nnz_current];
+ h_col_data_ = new index_type[nnz_current];
  owns_cpu_data_ = true;
  }
  if (h_val_data_ == nullptr) {
- h_val_data_ = new real_type[nnz_current]; 
+ h_val_data_ = new real_type[nnz_current];
  owns_cpu_vals_ = true;
  }
  mem_.copyArrayDeviceToHost(h_row_data_, d_row_data_, nnz_current);
@@ -246,7 +354,7 @@ namespace ReSolve
 
  /**
  * @brief Prints matrix data.
- * 
+ *
  * @param out - Output stream where the matrix data is printed
  */
  void matrix::Coo::print(std::ostream& out)

resolve/matrix/Coo.hpp

@@ -8,11 +8,21 @@ namespace ReSolve { namespace matrix {
  public:
  Coo();
  Coo(index_type n, index_type m, index_type nnz);
- Coo(index_type n, 
- index_type m, 
- index_type nnz,
- bool symmetric,
- bool expanded);
+ Coo(index_type n,
+ index_type m,
+ index_type nnz,
+ bool symmetric,
+ bool expanded);
+ Coo(index_type n,
+ index_type m,
+ index_type nnz,
+ bool symmetric,
+ bool expanded,
+ index_type** rows,
+ index_type** cols,
+ real_type** vals,
+ memory::MemorySpace memspaceSrc,
+ memory::MemorySpace memspaceDst);
  ~Coo();
 
  virtual index_type* getRowData(memory::MemorySpace memspace);