Fix alignment errors from cuSPARSE

cmake/ExternalHYPRE.cmake

@@ -141,7 +141,7 @@ if(PALACE_WITH_CUDA)
     "--with-cuda-home=${CUDA_DIR}"
     "--enable-curand"
     "--enable-cusparse"
-    "--enable-device-memory-pool"
+    # "--enable-device-memory-pool"
   )
   if(NOT "${CMAKE_CUDA_ARCHITECTURES}" STREQUAL "")
     list(GET CMAKE_CUDA_ARCHITECTURES 0 HYPRE_CUDA_ARCH)

cmake/ExternalMFEM.cmake

@@ -364,6 +364,7 @@ set(MFEM_PATCH_FILES
   "${CMAKE_SOURCE_DIR}/extern/patch/mfem/patch_par_tet_mesh_fix.diff"
   "${CMAKE_SOURCE_DIR}/extern/patch/mfem/patch_mesh_prism_vtu_fix.diff"
   "${CMAKE_SOURCE_DIR}/extern/patch/mfem/patch_workspace_vectors.diff"
+  "${CMAKE_SOURCE_DIR}/extern/patch/mfem/patch_workspace_vectors_2.diff"
   "${CMAKE_SOURCE_DIR}/extern/patch/mfem/patch_hypre_runtime_compute_policy.diff"
 )
 

extern/patch/mfem/patch_workspace_vectors_2.diff

@@ -0,0 +1,74 @@
+diff --git a/general/workspace.cpp b/general/workspace.cpp
+index ad8210a61..3ec077dfa 100644
+--- a/general/workspace.cpp
++++ b/general/workspace.cpp
+@@ -15,11 +15,11 @@ namespace mfem
+ {
+
+ WorkspaceVector::WorkspaceVector(
+-   internal::WorkspaceChunk &chunk_, int offset_, int n)
++   internal::WorkspaceChunk &chunk_, int offset_, int n, int padding)
+    : Vector(chunk_.GetData(), chunk_.GetOffset(), n),
+      chunk(chunk_),
+      offset(offset_),
+-     original_size(n)
++     original_size(n + padding)
+ {
+    UseDevice(true);
+ }
+@@ -48,8 +48,11 @@ WorkspaceChunk::WorkspaceChunk(int capacity)
+ WorkspaceVector WorkspaceChunk::NewVector(int n)
+ {
+    MFEM_ASSERT(HasCapacityFor(n), "Requested vector is too large.");
+-   WorkspaceVector vector(*this, offset, n);
+-   offset += n;
++   constexpr int alignment = 16;
++   const int s = (n * sizeof(double)) % alignment;
++   const int padding = (s > 0) ? (alignment - s) / sizeof(double) : 0;
++   WorkspaceVector vector(*this, offset, n, padding);
++   offset += n + padding;
+    vector_count += 1;
+    return vector;
+ }
+diff --git a/general/workspace.hpp b/general/workspace.hpp
+index b5c0ebae8..e41e71b11 100644
+--- a/general/workspace.hpp
++++ b/general/workspace.hpp
+@@ -50,7 +50,8 @@ class WorkspaceVector : public Vector
+    bool moved_from = false;
+
+    /// Private constructor, create with Workspace::NewVector() instead.
+-   WorkspaceVector(internal::WorkspaceChunk &chunk_, int offset_, int n);
++   WorkspaceVector(internal::WorkspaceChunk &chunk_, int offset_, int n,
++                   int padding = 0);
+
+ public:
+    /// @brief Move constructor. The moved-from WorkspaceVector has @a
+@@ -104,6 +105,9 @@ public:
+    /// Create a WorkspaceChunk with the given @a capacity.
+    WorkspaceChunk(int capacity);
+
++   /// Return the data offset.
++   int GetOffset() const { return offset; }
++
+    /// @brief Return the available capacity (i.e. the largest vector that will
+    /// fit in this chunk).
+    int GetAvailableCapacity() const { return data.Size() - offset; }
+@@ -115,15 +119,12 @@ public:
+    /// "original capacity" remains unchained.
+    int GetOriginalCapacity() const { return original_capacity; }
+
+-   /// Return the data offset.
+-   int GetOffset() const { return offset; }
++   /// Returns true if this chunk can fit a new vector of size @a n.
++   bool HasCapacityFor(int n) const { return n <= GetAvailableCapacity(); }
+
+    /// Sets whether the chunk is in the front of the list
+    void SetFront(bool front_) { front = front_; }
+
+-   /// Returns true if this chunk can fit a new vector of size @a n.
+-   bool HasCapacityFor(int n) const { return n <= GetAvailableCapacity(); }
+-
+    /// Returns true if this chunk is empty.
+    bool IsEmpty() const { return vector_count == 0; }
+

palace/linalg/vector.cpp

@@ -35,11 +35,6 @@ ComplexVector::ComplexVector(const std::complex<double> *py, int size, bool on_d
   Set(py, size, on_dev);
 }
 
-ComplexVector::ComplexVector(Vector &y, int offset, int size)
-{
-  MakeRef(y, offset, size);
-}
-
 void ComplexVector::UseDevice(bool use_dev)
 {
   xr.UseDevice(use_dev);
@@ -48,17 +43,13 @@ void ComplexVector::UseDevice(bool use_dev)
 
 void ComplexVector::SetSize(int size)
 {
-  xr.SetSize(size);
-  xi.SetSize(size);
-}
-
-void ComplexVector::MakeRef(Vector &y, int offset, int size)
-{
-  MFEM_ASSERT(y.Size() >= offset + 2 * size,
-              "Insufficient storage for ComplexVector alias reference of the given size!");
-  y.ReadWrite();  // Ensure memory is allocated on device before aliasing
-  xr.MakeRef(y, offset, size);
-  xi.MakeRef(y, offset + size, size);
+  if (size != Size())
+  {
+    MFEM_VERIFY(Size() == 0 || (xr.OwnsData() && xi.OwnsData()),
+                "Alias vectors should not be resized!");
+    xr.SetSize(size);
+    xi.SetSize(size);
+  }
 }
 
 void ComplexVector::Set(const ComplexVector &y)

palace/linalg/vector.hpp

@@ -21,7 +21,7 @@ using Vector = mfem::Vector;
 // A complex-valued vector represented as two real vectors, one for each component.
 class ComplexVector
 {
-private:
+protected:
   Vector xr, xi;
 
 public:
@@ -37,10 +37,6 @@ class ComplexVector
   // Copy constructor from an array of complex values.
   ComplexVector(const std::complex<double> *py, int size, bool on_dev);
 
-  // Create a vector referencing the memory of another vector, at the given base offset and
-  // size.
-  ComplexVector(Vector &y, int offset, int size);
-
   // Flag for runtime execution on the mfem::Device. See the documentation for mfem::Vector.
   void UseDevice(bool use_dev);
   bool UseDevice() const { return xr.UseDevice(); }
@@ -52,10 +48,6 @@ class ComplexVector
   // where the new size is less than or greater than Size() or Capacity().
   void SetSize(int size);
 
-  // Set this vector to reference the memory of another vector, at the given base offset and
-  // size.
-  void MakeRef(Vector &y, int offset, int size);
-
   // Get access to the real and imaginary vector parts.
   const Vector &Real() const { return xr; }
   Vector &Real() { return xr; }

palace/utils/workspace.hpp

@@ -18,12 +18,14 @@ namespace palace::workspace
 class ComplexWorkspaceVector : public ComplexVector
 {
 private:
-  mfem::WorkspaceVector data;
+  mfem::WorkspaceVector data_xr, data_xi;
 
 public:
-  ComplexWorkspaceVector(int size) : data(mfem::Workspace::NewVector(2 * size))
+  ComplexWorkspaceVector(int size)
+    : data_xr(mfem::Workspace::NewVector(size)), data_xi(mfem::Workspace::NewVector(size))
   {
-    MakeRef(data, 0, size);
+    xr.MakeRef(data_xr, 0, size);
+    xi.MakeRef(data_xi, 0, size);
   }
 
   // No copy constructor.