Add GPU scan support for Array<T>::PartialSum

Fix comments in scan.hpp
Add more device support to Array<T>
2025-10-21 09:57:37 -07:00 · 2025-10-21 09:57:37 -07:00 · 2025-10-21 09:57:37 -07:00 · 2025-10-21 09:57:37 -07:00 · 2025-10-20 15:38:54 -07:00 · 2025-10-18 10:48:44 -07:00
25 changed files with 1972 additions and 543 deletions
@@ -29,3 +29,47 @@ jobs:
        operations-per-run: 500
        exempt-issue-labels: "bug,WIP,ready-for-review,in-review,in-next"
        exempt-pr-labels: "bug,WIP,ready-for-review,in-review,in-next"
+
+  # Stale action for PRs with "in-review" label.
+  stale-in-review-pr:
+
+    runs-on: ubuntu-latest
+    permissions:
+      issues: write
+      pull-requests: write
+      actions: write
+
+    steps:
+    - uses: actions/stale@v9
+      with:
+        repo-token: ${{ secrets.GITHUB_TOKEN }}
+        stale-pr-message: ':warning: This PR has been automatically marked as stale because it has not had any activity in the last 150 days. *If no activity occurs in the next 30 days, it will be automatically closed.* Thank you for your contributions.'
+        only-pr-labels: "in-review"
+        days-before-pr-stale: 150
+        days-before-pr-close: 30
+        days-before-issue-stale: -1
+        days-before-issue-close: -1
+        stale-pr-label: 'stale'
+        operations-per-run: 500
+
+  # Stale action for PRs with "WIP" label.
+  stale-wip-pr:
+
+    runs-on: ubuntu-latest
+    permissions:
+      issues: write
+      pull-requests: write
+      actions: write
+
+    steps:
+    - uses: actions/stale@v9
+      with:
+        repo-token: ${{ secrets.GITHUB_TOKEN }}
+        stale-pr-message: ':warning: This PR has been automatically marked as stale because it has not had any activity in the last 300 days. *If no activity occurs in the next 30 days, it will be automatically closed.* Thank you for your contributions.'
+        only-pr-labels: "WIP"
+        days-before-pr-stale: 300
+        days-before-pr-close: 30
+        days-before-issue-stale: -1
+        days-before-issue-close: -1
+        stale-pr-label: 'stale'
+        operations-per-run: 500
@@ -197,16 +197,6 @@ if (HYPRE_FOUND AND HYPRE_USING_HIP)
  message(STATUS "Updated HYPRE_LIBRARIES: ${HYPRE_LIBRARIES}")
 endif()

-# Hypre+Umpire check
-if (HYPRE_FOUND AND (HYPRE_USING_CUDA OR HYPRE_USING_HIP) AND NOT MFEM_USE_UMPIRE)
-  message(WARNING
-"===============================================================
- Detected GPU-enabled HYPRE build without Umpire support.
- This is not recommended for performance reasons!
- Consider rebuilding HYPRE with Umpire support.
-===============================================================")
-endif()
-
 find_package_handle_standard_args(HYPRE
  REQUIRED_VARS HYPRE_LIBRARIES HYPRE_INCLUDE_DIRS HYPRE_VERSION
 )
@@ -57,7 +57,7 @@ CUDA_DIR = $(or $(CUDA_HOME),$(patsubst %/,%,$(dir \
 CLANG_CUDA_FLAGS = -xcuda --cuda-path=$(CUDA_DIR) --cuda-gpu-arch=$(CUDA_ARCH)
 # flags for nvcc
 NVCC_FLAGS = -x=cu --expt-extended-lambda --expt-relaxed-constexpr \
- -arch=$(CUDA_ARCH)
+ -arch=$(CUDA_ARCH) -isystem "$(CUDA_DIR)/include"
 # Prefixes for passing flags to the host compiler and linker when using
 # CUDA_CXX=nvcc
 CUDA_XCOMPILER = -Xcompiler=
@@ -0,0 +1,593 @@
+MFEM mesh v1.0
+
+# Created by: Pointwise
+
+# MFEM Geometry Types:
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+# PRISM       = 6
+
+dimension
+2
+
+elements
+160
+1 3 1 164 163 0
+1 3 164 165 162 163
+1 3 2 166 164 1
+1 3 166 132 165 164
+1 3 3 167 166 2
+1 3 167 131 132 166
+1 3 4 168 167 3
+1 3 168 130 131 167
+1 3 5 169 168 4
+1 3 169 129 130 168
+1 3 6 170 169 5
+1 3 170 128 129 169
+1 3 171 172 170 6
+1 3 172 127 128 170
+1 3 124 125 172 171
+1 3 125 126 127 172
+1 3 162 165 173 161
+1 3 165 132 133 173
+1 3 161 173 174 160
+1 3 173 133 134 174
+1 3 160 174 175 159
+1 3 174 134 135 175
+1 3 6 7 176 171
+1 3 7 8 177 176
+1 3 171 176 123 124
+1 3 176 177 122 123
+1 3 159 175 178 158
+1 3 175 135 136 178
+1 3 158 178 179 157
+1 3 178 136 137 179
+1 3 157 179 180 156
+1 3 179 137 138 180
+1 3 122 177 181 121
+1 3 177 8 182 181
+1 3 8 9 183 182
+1 3 9 10 184 183
+1 3 10 11 185 184
+1 3 11 12 186 185
+1 3 12 13 187 186
+1 3 13 14 15 187
+1 3 121 181 119 120
+1 3 181 182 118 119
+1 3 182 183 117 118
+1 3 183 184 188 117
+1 3 184 185 109 188
+1 3 185 186 108 109
+1 3 186 187 189 108
+1 3 187 15 16 189
+1 3 109 110 190 188
+1 3 110 111 191 190
+1 3 111 112 113 191
+1 3 188 190 116 117
+1 3 190 191 115 116
+1 3 191 113 114 115
+1 3 189 192 107 108
+1 3 192 193 106 107
+1 3 193 194 105 106
+1 3 194 195 104 105
+1 3 195 196 103 104
+1 3 16 17 192 189
+1 3 17 18 193 192
+1 3 18 19 194 193
+1 3 19 20 195 194
+1 3 20 21 196 195
+1 3 97 98 197 96
+1 3 98 99 198 197
+1 3 99 100 199 198
+1 3 100 101 200 199
+1 3 101 102 201 200
+1 3 102 103 202 201
+1 3 103 196 203 202
+1 3 196 21 22 203
+1 3 96 197 204 95
+1 3 197 198 39 204
+1 3 198 199 38 39
+1 3 199 200 205 38
+1 3 200 201 32 205
+1 3 201 202 31 32
+1 3 202 203 206 31
+1 3 203 22 23 206
+1 3 32 33 207 205
+1 3 33 34 35 207
+1 3 205 207 37 38
+1 3 207 35 36 37
+1 3 39 40 208 204
+1 3 40 41 209 208
+1 3 41 42 210 209
+1 3 42 43 211 210
+1 3 43 44 212 211
+1 3 204 208 94 95
+1 3 208 209 93 94
+1 3 209 210 92 93
+1 3 210 211 91 92
+1 3 211 212 90 91
+1 3 90 212 213 89
+1 3 212 44 214 213
+1 3 44 45 215 214
+1 3 45 46 216 215
+1 3 46 47 217 216
+1 3 47 48 218 217
+1 3 48 49 219 218
+1 3 49 50 51 219
+1 3 89 213 87 88
+1 3 213 214 86 87
+1 3 214 215 85 86
+1 3 215 216 84 85
+1 3 216 217 83 84
+1 3 217 218 82 83
+1 3 218 219 220 82
+1 3 219 51 52 220
+1 3 53 221 220 52
+1 3 221 81 82 220
+1 3 54 222 221 53
+1 3 222 80 81 221
+1 3 55 223 222 54
+1 3 223 79 80 222
+1 3 26 27 224 25
+1 3 27 28 29 224
+1 3 25 224 225 24
+1 3 224 29 30 225
+1 3 24 225 206 23
+1 3 225 30 31 206
+1 3 154 155 226 153
+1 3 155 156 180 226
+1 3 153 226 227 152
+1 3 226 180 138 227
+1 3 152 227 228 151
+1 3 227 138 139 228
+1 3 151 228 229 150
+1 3 228 139 140 229
+1 3 150 229 230 149
+1 3 229 140 141 230
+1 3 149 230 231 148
+1 3 230 141 142 231
+1 3 148 231 232 147
+1 3 231 142 143 232
+1 3 147 232 145 146
+1 3 232 143 144 145
+1 3 56 233 223 55
+1 3 233 78 79 223
+1 3 57 234 233 56
+1 3 234 77 78 233
+1 3 58 235 234 57
+1 3 235 76 77 234
+1 3 61 236 59 60
+1 3 236 235 58 59
+1 3 62 237 236 61
+1 3 237 76 235 236
+1 3 63 238 237 62
+1 3 238 75 76 237
+1 3 64 239 238 63
+1 3 239 74 75 238
+1 3 65 240 239 64
+1 3 240 73 74 239
+1 3 66 241 240 65
+1 3 241 72 73 240
+1 3 67 242 241 66
+1 3 242 71 72 241
+1 3 68 69 242 67
+1 3 69 70 71 242
+
+boundary
+164
+3 1 0 1
+3 1 1 2
+3 1 2 3
+3 1 3 4
+3 1 4 5
+3 1 5 6
+3 1 6 7
+3 1 7 8
+3 1 8 9
+3 1 9 10
+3 1 10 11
+3 1 11 12
+3 1 12 13
+3 1 13 14
+3 1 16 17
+3 1 17 18
+3 1 18 19
+3 1 19 20
+3 1 20 21
+3 1 21 22
+3 1 22 23
+3 1 23 24
+3 1 24 25
+3 1 25 26
+3 1 26 27
+3 1 27 28
+3 1 28 29
+3 1 29 30
+3 1 30 31
+3 1 31 32
+3 1 32 33
+3 1 33 34
+3 1 34 35
+3 1 35 36
+3 1 36 37
+3 1 37 38
+3 1 38 39
+3 1 39 40
+3 1 40 41
+3 1 41 42
+3 1 42 43
+3 1 43 44
+3 1 49 50
+3 1 48 49
+3 1 47 48
+3 1 46 47
+3 1 45 46
+3 1 44 45
+3 1 52 53
+3 1 53 54
+3 1 54 55
+3 1 57 58
+3 1 56 57
+3 1 55 56
+3 1 60 61
+3 1 61 62
+3 1 62 63
+3 1 63 64
+3 1 64 65
+3 1 65 66
+3 1 66 67
+3 1 67 68
+3 1 75 76
+3 1 74 75
+3 1 73 74
+3 1 72 73
+3 1 71 72
+3 1 70 71
+3 1 76 77
+3 1 77 78
+3 1 78 79
+3 1 81 82
+3 1 80 81
+3 1 79 80
+3 1 82 83
+3 1 83 84
+3 1 84 85
+3 1 85 86
+3 1 86 87
+3 1 87 88
+3 1 94 95
+3 1 93 94
+3 1 92 93
+3 1 91 92
+3 1 90 91
+3 1 96 97
+3 1 95 96
+3 1 97 98
+3 1 98 99
+3 1 99 100
+3 1 100 101
+3 1 101 102
+3 1 102 103
+3 1 107 108
+3 1 106 107
+3 1 105 106
+3 1 104 105
+3 1 103 104
+3 1 108 109
+3 1 109 110
+3 1 110 111
+3 1 111 112
+3 1 112 113
+3 1 113 114
+3 1 114 115
+3 1 115 116
+3 1 116 117
+3 1 119 120
+3 1 118 119
+3 1 117 118
+3 1 131 132
+3 1 130 131
+3 1 129 130
+3 1 128 129
+3 1 127 128
+3 1 126 127
+3 1 132 133
+3 1 133 134
+3 1 134 135
+3 1 137 138
+3 1 136 137
+3 1 135 136
+3 1 138 139
+3 1 139 140
+3 1 140 141
+3 1 141 142
+3 1 142 143
+3 1 143 144
+3 1 147 148
+3 1 146 147
+3 1 153 154
+3 1 152 153
+3 1 151 152
+3 1 150 151
+3 1 149 150
+3 1 148 149
+3 1 156 157
+3 1 157 158
+3 1 158 159
+3 1 161 162
+3 1 160 161
+3 1 159 160
+2 1 69 70
+2 1 68 69
+3 1 88 89
+3 1 89 90
+3 1 121 122
+3 1 120 121
+3 1 123 124
+3 1 122 123
+3 1 125 126
+3 1 124 125
+1 1 144 145
+1 1 145 146
+3 1 15 16
+3 1 14 15
+3 1 50 51
+3 1 51 52
+3 1 59 60
+3 1 58 59
+3 1 154 155
+3 1 155 156
+3 1 163 0
+3 1 162 163
+
+vertices
+243
+2
+4 4
+4 3.5
+4 3
+4 2.5
+4 2
+4 1.5
+4 1
+4.5 1
+5 1
+5 1.5
+5 2
+5 2.5
+5 3
+5 3.5
+5 4
+5.500 4
+6 4
+6.500 4
+7 4
+7.5 4
+8 4
+8.5 4
+9 4
+9.5 4
+10 4
+10.5 4
+11 4
+11 3.5
+11 3
+10.5 3
+10 3
+9.5 3
+9.5 2.5
+10 2.5
+10.5 2.5
+10.5 2
+10.5 1.5
+10 1.5
+9.5 1.5
+9.5 1
+10 1
+10.5 1
+11 1
+11.5 1
+12 1
+12 1.5
+12 2
+12 2.5
+12 3
+12 3.5
+12 4
+12.5 4
+13 4
+13.333 3.75
+13.666 3.5
+14.000 3.25
+14.333 3.5
+14.666 3.75
+15.000 4
+15.500 4
+16.000 4
+16.000 3.5
+16.000 3
+16.000 2.5
+16.000 2
+16.000 1.5
+16.000 1
+16.000 0.5
+16.000 0
+15.500 0
+15.000 0
+15.000 0.5000000000000002
+15.000 1
+15.000 1.5
+15.000 2
+15.000 2.5
+15.000 3
+14.666 2.75
+14.333 2.5
+14.000 2.25
+13.666 2.5
+13.333 2.75
+13 3
+13 2.5
+13 2
+13 1.5
+13 1
+13 0.500
+13 0
+12.5 0
+12 0
+11.5 0
+11 0
+10.5 0
+10 0
+9.5 0
+9 0
+8.5 0
+8.5 0.5
+8.5 1
+8.5 1.5
+8.5 2
+8.5 2.5
+8.5 3
+8 3
+7.5 3
+7 3
+6.500 3
+6 3
+6 2.5
+6.5 2.5
+7 2.5
+7.5 2.5
+7.5 2
+7.5 1.5
+7.000 1.5
+6.5 1.5
+6 1.5
+6 1
+6 0.5
+6 0
+5.5 0
+5 0
+4.5 0
+4 0
+3.5 0
+3 0
+3 0.500
+3 1
+3 1.5
+3 2
+3 2.5
+3 3
+2.666 2.75
+2.333 2.5
+2.000 2.25
+1.666 2.5
+1.333 2.75
+1.000 3
+1.000 2.5
+1.000 2
+1.000 1.5
+1.000 1
+1.000 0.5000
+1.000 0
+0.5000 0
+0.0000 0
+0.0000 0.5
+0.0000 1
+0.0000 1.5
+0.0000 2
+0.0000 2.5
+0.0000 3
+0.0000 3.5
+0.0000 4
+0.5000 4
+1.000 4
+1.333 3.75
+1.666 3.5
+2.000 3.25
+2.333 3.5
+2.666 3.75
+3 4
+3.5 4
+3.5 3.5
+3 3.5
+3.5 3
+3.5 2.5
+3.5 2
+3.5 1.5
+3.5 1
+4 0.5
+3.5 0.5
+2.666 3.25
+2.333 3
+2.000 2.75
+4.5 0.5
+5 0.5
+1.666 3
+1.333 3.25
+1.000 3.5
+5.5 0.5
+5.500 1
+5.500 1.5
+5.500 2
+5.500 2.5
+5.500 3
+5.500 3.5
+6 2
+6 3.5
+6.5 2
+7 2
+6.5 3.5
+7 3.5
+7.5 3.5
+8 3.5
+8.5 3.5
+9 0.5
+9 1
+9 1.5
+9 2
+9 2.5
+9 3
+9 3.5
+9.5 0.5
+9.5 2
+9.5 3.5
+10 2
+10 0.5
+10.5 0.5
+11 0.5
+11.5 0.5
+12 0.5
+12.5 0.500
+12.5 1
+12.5 1.5
+12.5 2
+12.5 2.5
+12.5 3
+12.5 3.5
+13 3.5
+13.333 3.250
+13.666 3
+14.000 2.75
+10.5 3.5
+10 3.5
+0.500 3.5
+0.500 3
+0.500 2.5
+0.500 2
+0.500 1.5
+0.500 1
+0.500 0.5
+14.333 3
+14.666 3.25
+15.000 3.5
+15.500 3.5
+15.500 3
+15.500 2.5
+15.500 2
+15.500 1.5
+15.500 1
+15.500 0.5
@@ -436,7 +436,7 @@ void NonlinearForm::Mult(const Vector &x, Vector &y) const
   // In parallel, the result is in 'py' which is an alias for 'aux2'.
 }

-Operator &NonlinearForm::GetGradient(const Vector &x) const
+Operator &NonlinearForm::GetGradient(const Vector &x, bool finalize) const
 {
   if (ext)
   {
@@ -644,6 +644,8 @@ Operator &NonlinearForm::GetGradient(const Vector &x) const
      }
   }

+   if (!finalize) { return *Grad; }
+
   if (!Grad->Finalized())
   {
      Grad->Finalize(skip_zeros);
@@ -1203,7 +1205,14 @@ const BlockVector &BlockNonlinearForm::Prolongate(const BlockVector &bx) const
      aux1.Update(block_offsets);
      for (int s = 0; s < fes.Size(); s++)
      {
-         P[s]->Mult(bx.GetBlock(s), aux1.GetBlock(s));
+         if (P[s])
+         {
+            P[s]->Mult(bx.GetBlock(s), aux1.GetBlock(s));
+         }
+         else
+         {
+            aux1.GetBlock(s) = bx.GetBlock(s);
+         }
      }
      return aux1;
   }
@@ -1232,11 +1241,16 @@ void BlockNonlinearForm::Mult(const Vector &x, Vector &y) const
      {
         cP[s]->MultTranspose(pby.GetBlock(s), by.GetBlock(s));
      }
+      else if (needs_prolongation)
+      {
+         by.GetBlock(s) = pby.GetBlock(s);
+      }
      by.GetBlock(s).SetSubVector(*ess_tdofs[s], 0.0);
   }
 }

-void BlockNonlinearForm::ComputeGradientBlocked(const BlockVector &bx) const
+void BlockNonlinearForm::ComputeGradientBlocked(const BlockVector &bx,
+                                                bool finalize) const
 {
   const int skip_zeros = 0;
   Array<Array<int> *> vdofs(fes.Size());
@@ -1490,7 +1504,7 @@ void BlockNonlinearForm::ComputeGradientBlocked(const BlockVector &bx) const
      }
   }

-   if (!Grads(0,0)->Finalized())
+   if (finalize && !Grads(0,0)->Finalized())
   {
      for (int i=0; i<fes.Size(); ++i)
      {
@@ -1529,7 +1543,23 @@ Operator &BlockNonlinearForm::GetGradient(const Vector &x) const
         for (int s2 = 0; s2 < fes.Size(); ++s2)
         {
            delete cGrads(s1, s2);
-            cGrads(s1, s2) = RAP(*cP[s1], *Grads(s1, s2), *cP[s2]);
+            if (cP[s1] && cP[s2])
+            {
+               cGrads(s1, s2) = RAP(*cP[s1], *Grads(s1, s2), *cP[s2]);
+            }
+            else if (cP[s1])
+            {
+               cGrads(s1, s2) = TransposeMult(*cP[s1], *Grads(s1, s2));
+            }
+            else if (cP[s2])
+            {
+               cGrads(s1, s2) = mfem::Mult(*Grads(s1, s2), *cP[s2]);
+            }
+            else
+            {
+               cGrads(s1, s2) = NULL;
+               continue;
+            }
            mGrads(s1, s2) = cGrads(s1, s2);
         }
      }
@@ -217,7 +217,12 @@ public:
       In general, @a x may have non-homogeneous essential boundary values.

       The state @a x must be a true-dof vector. */
-   Operator &GetGradient(const Vector &x) const override;
+   Operator &GetGradient(const Vector &x) const override { return GetGradient(x, true); }
+
+   /** @brief Compute the gradient Operator of the NonlinearForm corresponding
+       to the state @a x with optional finalization and elimintaion. */
+   /** @see GetGradient(const Vector &) */
+   Operator &GetGradient(const Vector &x, bool finalize) const;

   /// Update the NonlinearForm to propagate updates of the associated FE space.
   /** After calling this method, the essential boundary conditions need to be
@@ -308,7 +313,7 @@ protected:
   void MultBlocked(const BlockVector &bx, BlockVector &by) const;

   /// Specialized version of GetGradient() for BlockVector
-   void ComputeGradientBlocked(const BlockVector &bx) const;
+   void ComputeGradientBlocked(const BlockVector &bx, bool finalize = true) const;

 public:
   /// Construct an empty BlockNonlinearForm. Initialize with SetSpaces().
@@ -151,6 +151,15 @@ void ParBilinearForm::ParallelRAP(SparseMatrix &loc_A, OperatorHandle &A,
   }
 }

+HypreParMatrix *ParBilinearForm::ParallelAssembleInternalMatrix()
+{
+   if (p_mat.Ptr() == NULL)
+   {
+      ParallelAssemble(p_mat, mat);
+   }
+   return p_mat.As<HypreParMatrix>();
+}
+
 void ParBilinearForm::ParallelAssemble(OperatorHandle &A, SparseMatrix *A_local)
 {
   A.Clear();
@@ -333,6 +342,15 @@ void ParBilinearForm
   A.EliminateRowsCols(dof_list, X, B);
 }

+void ParBilinearForm::ParallelEliminateEssentialBC(
+   const Array<int> &bdr_attr_is_ess, const HypreParVector &X, HypreParVector &B)
+{
+   Array<int> dof_list;
+   pfes->GetEssentialTrueDofs(bdr_attr_is_ess, dof_list);
+
+   p_mat.As<HypreParMatrix>()->EliminateRowsCols(dof_list, X, B);
+}
+
 HypreParMatrix *ParBilinearForm::
 ParallelEliminateEssentialBC(const Array<int> &bdr_attr_is_ess,
                             HypreParMatrix &A) const
@@ -344,6 +362,26 @@ ParallelEliminateEssentialBC(const Array<int> &bdr_attr_is_ess,
   return A.EliminateRowsCols(dof_list);
 }

+void ParBilinearForm::ParallelEliminateEssentialBC(const Array<int>
+                                                   &bdr_attr_is_ess)
+{
+   Array<int> tdofs_list;
+   pfes->GetEssentialTrueDofs(bdr_attr_is_ess, tdofs_list);
+
+   ParallelEliminateTDofs(tdofs_list);
+}
+
+void ParBilinearForm::ParallelEliminateTDofs(const Array<int> &tdofs_list)
+{
+   p_mat_e.EliminateRowsCols(p_mat, tdofs_list);
+}
+
+void ParBilinearForm::ParallelEliminateTDofsInRHS(
+   const Array<int> &tdofs_list, const Vector &x, Vector &b)
+{
+   p_mat.EliminateBC(p_mat_e, tdofs_list, x, b);
+}
+
 void ParBilinearForm::TrueAddMult(const Vector &x, Vector &y, const real_t a)
 const
 {
@@ -485,7 +523,7 @@ void ParBilinearForm::FormLinearSystem(
      HypreParVector true_X(pfes), true_B(pfes);
      P.MultTranspose(b, true_B);
      R.Mult(x, true_X);
-      p_mat.EliminateBC(p_mat_e, ess_tdof_list, true_X, true_B);
+      ParallelEliminateTDofsInRHS(ess_tdof_list, true_X, true_B);
      R.MultTranspose(true_B, b);
      hybridization->ReduceRHS(true_B, B);
      X.SetSize(B.Size());
@@ -498,17 +536,11 @@ void ParBilinearForm::FormLinearSystem(
      B.SetSize(X.Size());
      P.MultTranspose(b, B);
      R.Mult(x, X);
-      p_mat.EliminateBC(p_mat_e, ess_tdof_list, X, B);
+      ParallelEliminateTDofsInRHS(ess_tdof_list, X, B);
      if (!copy_interior) { X.SetSubVectorComplement(ess_tdof_list, 0.0); }
   }
 }

-void ParBilinearForm::EliminateVDofsInRHS(
-   const Array<int> &vdofs, const Vector &x, Vector &b)
-{
-   p_mat.EliminateBC(p_mat_e, vdofs, x, b);
-}
-
 void ParBilinearForm::FormSystemMatrix(const Array<int> &ess_tdof_list,
                                       OperatorHandle &A)
 {
@@ -553,7 +585,7 @@ void ParBilinearForm::FormSystemMatrix(const Array<int> &ess_tdof_list,
         mat = NULL;
         delete mat_e;
         mat_e = NULL;
-         p_mat_e.EliminateRowsCols(p_mat, ess_tdof_list);
+         ParallelEliminateTDofs(ess_tdof_list);
      }
      if (hybridization)
      {
@@ -615,36 +647,180 @@ void ParBilinearForm::Update(FiniteElementSpace *nfes)
   p_mat_e.Clear();
 }

-
-HypreParMatrix *ParMixedBilinearForm::ParallelAssemble()
+void ParMixedBilinearForm::pAllocMat()
 {
-   // construct the block-diagonal matrix A
-   HypreParMatrix *A =
-      new HypreParMatrix(trial_pfes->GetComm(),
-                         test_pfes->GlobalVSize(),
-                         trial_pfes->GlobalVSize(),
-                         test_pfes->GetDofOffsets(),
-                         trial_pfes->GetDofOffsets(),
-                         mat);
+   const int trial_nbr_size = trial_pfes->GetFaceNbrVSize();
+   const int test_nbr_size = test_pfes->GetFaceNbrVSize();

-   HypreParMatrix *rap = RAP(test_pfes->Dof_TrueDof_Matrix(), A,
-                             trial_pfes->Dof_TrueDof_Matrix());
-
-   delete A;
-
-   return rap;
+   if (keep_nbr_block)
+   {
+      mat = new SparseMatrix(height + test_nbr_size, width + trial_nbr_size);
+   }
+   else
+   {
+      mat = new SparseMatrix(height, width + trial_nbr_size);
+   }
 }

-void ParMixedBilinearForm::ParallelAssemble(OperatorHandle &A)
+void ParMixedBilinearForm::AssembleSharedFaces(int skip_zeros)
 {
-   // construct the rectangular block-diagonal matrix dA
-   OperatorHandle dA(A.Type());
-   dA.MakeRectangularBlockDiag(trial_pfes->GetComm(),
-                               test_pfes->GlobalVSize(),
-                               trial_pfes->GlobalVSize(),
-                               test_pfes->GetDofOffsets(),
-                               trial_pfes->GetDofOffsets(),
-                               mat);
+   ParMesh *pmesh = trial_pfes->GetParMesh();
+   FaceElementTransformations *T;
+   Array<int> tr_vdofs1, tr_vdofs2, tr_vdofs_all;
+   Array<int> te_vdofs1, te_vdofs2, te_vdofs_all;
+   DenseMatrix elemmat;
+
+   int nfaces = pmesh->GetNSharedFaces();
+   for (int i = 0; i < nfaces; i++)
+   {
+      T = pmesh->GetSharedFaceTransformations(i);
+      int Elem2NbrNo = T->Elem2No - pmesh->GetNE();
+      trial_pfes->GetElementVDofs(T->Elem1No, tr_vdofs1);
+      test_pfes->GetElementVDofs(T->Elem1No, te_vdofs1);
+      trial_pfes->GetFaceNbrElementVDofs(Elem2NbrNo, tr_vdofs2);
+      test_pfes->GetFaceNbrElementVDofs(Elem2NbrNo, te_vdofs2);
+
+      tr_vdofs1.Copy(tr_vdofs_all);
+      for (int j = 0; j < tr_vdofs2.Size(); j++)
+      {
+         if (tr_vdofs2[j] >= 0)
+         {
+            tr_vdofs2[j] += width;
+         }
+         else
+         {
+            tr_vdofs2[j] -= width;
+         }
+      }
+      tr_vdofs_all.Append(tr_vdofs2);
+
+      if (keep_nbr_block)
+      {
+         te_vdofs1.Copy(te_vdofs_all);
+         for (int j = 0; j < te_vdofs2.Size(); j++)
+         {
+            if (te_vdofs2[j] >= 0)
+            {
+               te_vdofs2[j] += height;
+            }
+            else
+            {
+               te_vdofs2[j] -= height;
+            }
+         }
+         te_vdofs_all.Append(te_vdofs2);
+      }
+
+      for (int k = 0; k < interior_face_integs.Size(); k++)
+      {
+         interior_face_integs[k]->
+         AssembleFaceMatrix(*trial_pfes->GetFE(T->Elem1No),
+                            *test_pfes->GetFE(T->Elem1No),
+                            *trial_pfes->GetFaceNbrFE(Elem2NbrNo),
+                            *test_pfes->GetFaceNbrFE(Elem2NbrNo),
+                            *T, elemmat);
+         if (keep_nbr_block)
+         {
+            mat->AddSubMatrix(te_vdofs_all, tr_vdofs_all, elemmat, skip_zeros);
+         }
+         else
+         {
+            mat->AddSubMatrix(te_vdofs1, tr_vdofs_all, elemmat, skip_zeros);
+         }
+      }
+   }
+}
+
+void ParMixedBilinearForm::Assemble(int skip_zeros)
+{
+   if (interior_face_integs.Size())
+   {
+      trial_pfes->ExchangeFaceNbrData();
+      test_pfes->ExchangeFaceNbrData();
+      if (!ext && mat == NULL)
+      {
+         pAllocMat();
+      }
+   }
+
+   MixedBilinearForm::Assemble(skip_zeros);
+
+   if (!ext && interior_face_integs.Size() > 0)
+   {
+      AssembleSharedFaces(skip_zeros);
+   }
+}
+
+HypreParMatrix *ParMixedBilinearForm::ParallelAssembleInternalMatrix()
+{
+   if (p_mat.Ptr() == NULL)
+   {
+      ParallelAssemble(p_mat, mat);
+   }
+   return p_mat.As<HypreParMatrix>();
+}
+
+HypreParMatrix *ParMixedBilinearForm::ParallelAssemble(SparseMatrix *m)
+{
+   OperatorHandle Mh(Operator::Hypre_ParCSR);
+   ParallelAssemble(Mh, m);
+   Mh.SetOperatorOwner(false);
+   return Mh.As<HypreParMatrix>();
+}
+
+void ParMixedBilinearForm::ParallelAssemble(OperatorHandle &A,
+                                            SparseMatrix *A_local)
+{
+   A.Clear();
+
+   if (A_local == NULL) { return; }
+   MFEM_VERIFY(A_local->Finalized(), "the local matrix must be finalized");
+
+   OperatorHandle dA(A.Type()), hdA;
+
+   if (interior_face_integs.Size() == 0)
+   {
+      // construct the rectangular block-diagonal matrix dA
+      dA.MakeRectangularBlockDiag(trial_pfes->GetComm(),
+                                  test_pfes->GlobalVSize(),
+                                  trial_pfes->GlobalVSize(),
+                                  test_pfes->GetDofOffsets(),
+                                  trial_pfes->GetDofOffsets(),
+                                  A_local);
+   }
+   else
+   {
+      // handle the case when 'a' contains off-diagonal
+      const int lvrows = test_pfes->GetVSize();
+      const int lvcols = trial_pfes->GetVSize();
+      const HYPRE_BigInt *face_nbr_glob_lcol = trial_pfes->GetFaceNbrGlobalDofMap();
+      const HYPRE_BigInt lcol_offset = trial_pfes->GetMyDofOffset();
+
+      Array<HYPRE_BigInt> glob_J(A_local->NumNonZeroElems());
+      const int *J = A_local->GetJ();
+      for (int i = 0; i < glob_J.Size(); i++)
+      {
+         if (J[i] < lvcols)
+         {
+            glob_J[i] = J[i] + lcol_offset;
+         }
+         else
+         {
+            glob_J[i] = face_nbr_glob_lcol[J[i] - lvcols];
+         }
+      }
+
+      // TODO - construct dA directly in the A format
+      hdA.Reset(
+         new HypreParMatrix(trial_pfes->GetComm(), lvrows, test_pfes->GlobalVSize(),
+                            trial_pfes->GlobalVSize(), A_local->GetI(), glob_J,
+                            A_local->GetData(), test_pfes->GetDofOffsets(),
+                            trial_pfes->GetDofOffsets()));
+      // - hdA owns the new HypreParMatrix
+      // - the above constructor copies all input arrays
+      glob_J.DeleteAll();
+      dA.ConvertFrom(hdA);
+   }

   OperatorHandle P_test(A.Type()), P_trial(A.Type());

@@ -670,6 +846,44 @@ void ParMixedBilinearForm::TrueAddMult(const Vector &x, Vector &y,
   test_pfes->Dof_TrueDof_Matrix()->MultTranspose(a, Yaux, 1.0, y);
 }

+void ParMixedBilinearForm::ParallelEliminateTrialEssentialBC(
+   const Array<int> &bdr_attr_is_ess)
+{
+   Array<int> trial_tdof_list;
+   trial_pfes->GetEssentialTrueDofs(bdr_attr_is_ess, trial_tdof_list);
+
+   ParallelEliminateTrialTDofs(trial_tdof_list);
+}
+
+void ParMixedBilinearForm::ParallelEliminateTrialTDofs(
+   const Array<int> &trial_tdof_list)
+{
+   HypreParMatrix *temp = p_mat.As<HypreParMatrix>()->EliminateCols(
+                             trial_tdof_list);
+   p_mat_e.Reset(temp, true);
+}
+
+void ParMixedBilinearForm::ParallelEliminateTrialTDofsInRHS(
+   const Array<int> &trial_tdof_list, const Vector &x, Vector &b)
+{
+   p_mat_e.As<HypreParMatrix>()->Mult(-1.0, x, 1.0, b);
+}
+
+void ParMixedBilinearForm::ParallelEliminateTestEssentialBC(
+   const Array<int> &bdr_attr_is_ess)
+{
+   Array<int> test_tdof_list;
+   test_pfes->GetEssentialTrueDofs(bdr_attr_is_ess, test_tdof_list);
+
+   ParallelEliminateTestTDofs(test_tdof_list);
+}
+
+void ParMixedBilinearForm::ParallelEliminateTestTDofs(
+   const Array<int> &test_tdof_list)
+{
+   p_mat.As<HypreParMatrix>()->EliminateRows(test_tdof_list);
+}
+
 void ParMixedBilinearForm::FormRectangularSystemMatrix(
   const Array<int>
   &trial_tdof_list,
@@ -690,10 +904,8 @@ void ParMixedBilinearForm::FormRectangularSystemMatrix(
      mat = NULL;
      delete mat_e;
      mat_e = NULL;
-      HypreParMatrix *temp =
-         p_mat.As<HypreParMatrix>()->EliminateCols(trial_tdof_list);
-      p_mat.As<HypreParMatrix>()->EliminateRows(test_tdof_list);
-      p_mat_e.Reset(temp, true);
+      ParallelEliminateTrialTDofs(trial_tdof_list);
+      ParallelEliminateTestTDofs(test_tdof_list);
   }

   A = p_mat;
@@ -723,7 +935,7 @@ void ParMixedBilinearForm::FormRectangularLinearSystem(
   test_P->MultTranspose(b, B);
   trial_R->Mult(x, X);

-   p_mat_e.As<HypreParMatrix>()->Mult(-1.0, X, 1.0, B);
+   ParallelEliminateTrialTDofsInRHS(trial_tdof_list, X, B);
   B.SetSubVector(test_tdof_list, 0.0);
 }

@@ -73,7 +73,7 @@ public:
   /** When set to true and the ParBilinearForm has interior face integrators,
       the local SparseMatrix will include the rows (in addition to the columns)
       corresponding to face-neighbor dofs. The default behavior is to disregard
-       those rows. Must be called before the first Assemble call. */
+       those rows. Must be called before the first Assemble() call. */
   void KeepNbrBlock(bool knb = true) { keep_nbr_block = knb; }

   /** @brief Set the operator type id for the parallel matrix/operator when
@@ -101,6 +101,14 @@ public:
       diagonal for this case. */
   void AssembleDiagonal(Vector &diag) const override;

+   /// Returns the matrix assembled on the true dofs, i.e. P^t A P.
+   /** The returned matrix is the internal one, owned by the form. It is not
+       reassembled if it has been already constructed. If FormSystemMatrix()
+       has been called before, it is the system matrix with eliminated
+       essential DOFs, otherwise the parallel matrix is assembled here without
+       the elimination process. */
+   HypreParMatrix *ParallelAssembleInternalMatrix();
+
   /// Returns the matrix assembled on the true dofs, i.e. P^t A P.
   /** The returned matrix has to be deleted by the caller. */
   HypreParMatrix *ParallelAssemble() { return ParallelAssemble(mat); }
@@ -146,6 +154,13 @@ public:
                                     const HypreParVector &X,
                                     HypreParVector &B) const;

+   /// Eliminate essential boundary DOFs from the parallel system matrix.
+   /** The array @a bdr_attr_is_ess marks boundary attributes that constitute
+       the essential part of the boundary. */
+   void ParallelEliminateEssentialBC(const Array<int> &bdr_attr_is_ess,
+                                     const HypreParVector &X,
+                                     HypreParVector &B);
+
   /// Eliminate essential boundary DOFs from a parallel assembled matrix @a A.
   /** The array @a bdr_attr_is_ess marks boundary attributes that constitute
       the essential part of the boundary. The eliminated part is stored in a
@@ -157,6 +172,12 @@ public:
   HypreParMatrix *ParallelEliminateEssentialBC(const Array<int> &bdr_attr_is_ess,
                                                HypreParMatrix &A) const;

+   /// Eliminate essential boundary DOFs from the parallel system matrix.
+   /** The array @a bdr_attr_is_ess marks boundary attributes that constitute
+       the essential part of the boundary. This method relies on
+       ParallelEliminateTDofs(const Array<int> &), see it for details. */
+   void ParallelEliminateEssentialBC(const Array<int> &bdr_attr_is_ess);
+
   /// Eliminate essential true DOFs from a parallel assembled matrix @a A.
   /** Given a list of essential true dofs and the parallel assembled matrix
       @a A, eliminate the true dofs from the matrix, storing the eliminated
@@ -169,6 +190,28 @@ public:
                                          HypreParMatrix &A) const
   { return A.EliminateRowsCols(tdofs_list); }

+   /// Eliminate essential true DOFs from the parallel system matrix.
+   /** Given a list of essential true dofs, eliminate the true dofs from
+       the parallel assembled system matrix, storing the eliminated part
+       internally. This method works in conjunction with
+       ParallelEliminateTDofsInRHS() and allows elimination of boundary
+       conditions in multiple right-hand sides. */
+   void ParallelEliminateTDofs(const Array<int> &tdofs_list);
+
+   /** @brief Use the stored eliminated part of the parallel system matrix for
+       elimination of boundary conditions in the r.h.s. */
+   /** Given a list of essential true dofs, eliminate the true dofs from the
+       right-hand side @a b using the solution vector @a x and the previously
+       stored eliminated part of the parallel assembled system matrix produced
+       by ParallelEliminateTDofs(const Array<int> &). */
+   void ParallelEliminateTDofsInRHS(const Array<int> &tdofs, const Vector &x,
+                                    Vector &b);
+
+   /// @deprecated Use ParallelEliminateTDofsInRHS() instead.
+   MFEM_DEPRECATED void EliminateVDofsInRHS(const Array<int> &vdofs,
+                                            const Vector &x, Vector &b)
+   { ParallelEliminateTDofsInRHS(vdofs, x, b); }
+
   /** @brief Compute @a y += @a a (P^t A P) @a x, where @a x and @a y are
       vectors on the true dofs. */
   void TrueAddMult(const Vector &x, Vector &y, const real_t a = 1.0) const;
@@ -238,8 +281,6 @@ public:

   void Update(FiniteElementSpace *nfes = NULL) override;

-   void EliminateVDofsInRHS(const Array<int> &vdofs, const Vector &x, Vector &b);
-
   virtual ~ParBilinearForm() { }
 };

@@ -257,6 +298,13 @@ protected:
   /// Matrix and eliminated matrix
   OperatorHandle p_mat, p_mat_e;

+   bool keep_nbr_block;
+
+   // Allocate mat - called when (mat == NULL && fbfi.Size() > 0)
+   void pAllocMat();
+
+   void AssembleSharedFaces(int skip_zeros = 1);
+
 private:
   /// Copy construction is not supported; body is undefined.
   ParMixedBilinearForm(const ParMixedBilinearForm &);
@@ -276,6 +324,7 @@ public:
   {
      trial_pfes = trial_fes;
      test_pfes  = test_fes;
+      keep_nbr_block = false;
   }

   /** @brief Create a ParMixedBilinearForm on the given FiniteElementSpace%s
@@ -295,15 +344,89 @@ public:
   {
      trial_pfes = trial_fes;
      test_pfes  = test_fes;
+      keep_nbr_block = false;
   }

+   /** When set to true and the ParMixedBilinearForm has interior face
+       integrators, the local SparseMatrix will include the rows (in addition
+       to the columns) corresponding to face-neighbor dofs. The default
+       behavior is to disregard those rows. Must be called before the first
+       Assemble() call. */
+   void KeepNbrBlock(bool knb = true) { keep_nbr_block = knb; }
+
+   /// Assemble the local matrix
+   void Assemble(int skip_zeros = 1);
+
   /// Returns the matrix assembled on the true dofs, i.e. P_test^t A P_trial.
-   HypreParMatrix *ParallelAssemble();
+   /** The returned matrix is the internal one, owned by the form. It is not
+       reassembled if it has been already constructed. If
+       FormRectangularSystemMatrix() has been called before, it is the system
+       matrix with eliminated essential DOFs, otherwise the parallel matrix is
+       assembled here without the elimination process. */
+   HypreParMatrix *ParallelAssembleInternalMatrix();
+
+   /// Returns the matrix assembled on the true dofs, i.e. P_test^t A P_trial.
+   /** The returned matrix has to be deleted by the caller. */
+   HypreParMatrix *ParallelAssemble() { return ParallelAssemble(mat); }
+
+   /** @brief Returns the eliminated matrix assembled on the true dofs, i.e.
+       P_test^t A_local P_trial. */
+   /** The returned matrix has to be deleted by the caller. */
+   HypreParMatrix *ParallelAssembleElim() { return ParallelAssemble(mat_e); }
+
+   /** @brief Return the matrix @a m assembled on the true dofs, i.e. P_test^t
+       A_local P_trial. */
+   /** The returned matrix has to be deleted by the caller. */
+   HypreParMatrix *ParallelAssemble(SparseMatrix *m);

   /** @brief Returns the matrix assembled on the true dofs, i.e.
       @a A = P_test^t A_local P_trial, in the format (type id) specified by
       @a A. */
-   void ParallelAssemble(OperatorHandle &A);
+   void ParallelAssemble(OperatorHandle &A) { ParallelAssemble(A, mat); }
+
+   /** Returns the eliminated matrix assembled on the true dofs, i.e.
+       @a A_elim = P^t A_elim_local P in the format (type id) specified by @a A.
+    */
+   void ParallelAssembleElim(OperatorHandle &A_elim)
+   { ParallelAssemble(A_elim, mat_e); }
+
+   /** Returns the matrix @a A_local assembled on the true dofs, i.e.
+       @a A = P_test^t A_local P_trial in the format (type id) specified by
+       @a A. */
+   void ParallelAssemble(OperatorHandle &A, SparseMatrix *A_local);
+
+   /// Eliminate essential boundary trial DOFs from the parallel system matrix.
+   /** The array @a bdr_attr_is_ess marks boundary attributes that constitute
+       the essential part of the boundary. This method relies on
+       ParallelEliminateTrialTDofs(const Array<int> &), see it for details. */
+   void ParallelEliminateTrialEssentialBC(const Array<int> &bdr_attr_is_ess);
+
+   /// Eliminate essential trial true DOFs from the parallel system matrix.
+   /** Given a list of essential trial true dofs, eliminate the trial true dofs
+       from the parallel assembled system matrix, storing the eliminated part
+       internally. This method works in conjunction with
+       ParallelEliminateTrialTDofsInRHS() and allows elimination of boundary
+       conditions in multiple right-hand sides. */
+   void ParallelEliminateTrialTDofs(const Array<int> &trial_tdof_list);
+
+   /** @brief Use the stored eliminated part of the parallel system matrix for
+       elimination of boundary conditions in the r.h.s. */
+   /** Given a list of essential trial true dofs, eliminate the trial true dofs
+       from the right-hand side @a B using the solution vector @a X and the
+       previously stored eliminated part of the parallel assembled system
+       matrix produced by ParallelEliminateTrialTDofs(const Array<int> &). */
+   void ParallelEliminateTrialTDofsInRHS(const Array<int> &trial_tdof_list,
+                                         const Vector &X, Vector &B);
+
+   /// Eliminate essential boundary test DOFs from the parallel system matrix.
+   /** The array @a bdr_attr_is_ess marks boundary attributes that constitute
+       the essential part of the boundary. */
+   void ParallelEliminateTestEssentialBC(const Array<int> &bdr_attr_is_ess);
+
+   /// Eliminate essential test true DOFs from the parallel system matrix.
+   /** Given a list of essential test true dofs, eliminate the test true dofs
+       from the parallel assembled system matrix. */
+   void ParallelEliminateTestTDofs(const Array<int> &test_tdof_list);

   using MixedBilinearForm::FormRectangularSystemMatrix;
   using MixedBilinearForm::FormRectangularLinearSystem;
@@ -105,6 +105,59 @@ const SparseMatrix &ParNonlinearForm::GetLocalGradient(const Vector &x) const
   return *Grad;
 }

+void ParNonlinearForm::GradientSharedFaces(const Vector &x,
+                                           int skip_zeros) const
+{
+   ParFiniteElementSpace *pfes = ParFESpace();
+   ParMesh *pmesh = pfes->GetParMesh();
+   FaceElementTransformations *T;
+   Array<int> vdofs1, vdofs2, vdofs_all;
+   DenseMatrix elemmat;
+   Vector el_x, nbr_x, face_x;
+   const Vector &px = Prolongate(x);
+
+   ParGridFunction pgf(pfes, const_cast<Vector&>(px), 0);
+   pgf.ExchangeFaceNbrData();
+
+   int nfaces = pmesh->GetNSharedFaces();
+   for (int i = 0; i < nfaces; i++)
+   {
+      T = pmesh->GetSharedFaceTransformations(i);
+      int Elem2NbrNo = T->Elem2No - pmesh->GetNE();
+
+      pfes->GetElementVDofs(T->Elem1No, vdofs1);
+      pfes->GetFaceNbrElementVDofs(Elem2NbrNo, vdofs2);
+      face_x.SetSize(vdofs1.Size() + vdofs2.Size());
+
+      el_x.MakeRef(face_x, 0, vdofs1.Size());
+      pgf.GetSubVector(vdofs1, el_x);
+
+      nbr_x.MakeRef(face_x, vdofs1.Size(), vdofs2.Size());
+      pgf.FaceNbrData().GetSubVector(vdofs2, nbr_x);
+
+      vdofs1.Copy(vdofs_all);
+      for (int j = 0; j < vdofs2.Size(); j++)
+      {
+         if (vdofs2[j] >= 0)
+         {
+            vdofs2[j] += height;
+         }
+         else
+         {
+            vdofs2[j] -= height;
+         }
+      }
+      vdofs_all.Append(vdofs2);
+      for (int k = 0; k < fnfi.Size(); k++)
+      {
+         fnfi[k]->AssembleFaceGrad(*pfes->GetFE(T->Elem1No),
+                                   *pfes->GetFaceNbrFE(Elem2NbrNo),
+                                   *T, face_x, elemmat);
+         Grad->AddSubMatrix(vdofs1, vdofs_all, elemmat, skip_zeros);
+      }
+   }
+}
+
 Operator &ParNonlinearForm::GetGradient(const Vector &x) const
 {
   if (NonlinearForm::ext) { return NonlinearForm::GetGradient(x); }
@@ -112,19 +165,61 @@ Operator &ParNonlinearForm::GetGradient(const Vector &x) const
   ParFiniteElementSpace *pfes = ParFESpace();

   pGrad.Clear();
+   OperatorHandle dA(pGrad.Type()), Ph(pGrad.Type()), hdA;

-   NonlinearForm::GetGradient(x); // (re)assemble Grad, no b.c.
-
-   OperatorHandle dA(pGrad.Type()), Ph(pGrad.Type());
-
-   if (fnfi.Size() == 0)
+   if (fnfi.Size())
   {
-      dA.MakeSquareBlockDiag(pfes->GetComm(), pfes->GlobalVSize(),
-                             pfes->GetDofOffsets(), Grad);
+      const int skip_zeros = 0;
+
+      pfes->ExchangeFaceNbrData();
+      if (Grad == NULL)
+      {
+         int nbr_size = pfes->GetFaceNbrVSize();
+         Grad = new SparseMatrix(pfes->GetVSize(), pfes->GetVSize() + nbr_size);
+      }
+
+      NonlinearForm::GetGradient(x, false); // (re)assemble Grad, no b.c.
+
+      GradientSharedFaces(x, skip_zeros);
+
+      Grad->Finalize(skip_zeros);
+
+      // handle the case when 'a' contains off-diagonal
+      int lvsize = pfes->GetVSize();
+      const HYPRE_BigInt *face_nbr_glob_ldof = pfes->GetFaceNbrGlobalDofMap();
+      HYPRE_BigInt ldof_offset = pfes->GetMyDofOffset();
+
+      Array<HYPRE_BigInt> glob_J(Grad->NumNonZeroElems());
+      int *J = Grad->GetJ();
+      for (int i = 0; i < glob_J.Size(); i++)
+      {
+         if (J[i] < lvsize)
+         {
+            glob_J[i] = J[i] + ldof_offset;
+         }
+         else
+         {
+            glob_J[i] = face_nbr_glob_ldof[J[i] - lvsize];
+         }
+      }
+
+      // TODO - construct dA directly in the A format
+      hdA.Reset(
+         new HypreParMatrix(pfes->GetComm(), lvsize, pfes->GlobalVSize(),
+                            pfes->GlobalVSize(), Grad->GetI(), glob_J,
+                            Grad->GetData(), pfes->GetDofOffsets(),
+                            pfes->GetDofOffsets()));
+      // - hdA owns the new HypreParMatrix
+      // - the above constructor copies all input arrays
+      glob_J.DeleteAll();
+      dA.ConvertFrom(hdA);
   }
   else
   {
-      MFEM_ABORT("TODO: assemble contributions from shared face terms");
+      NonlinearForm::GetGradient(x); // (re)assemble Grad, no b.c.
+
+      dA.MakeSquareBlockDiag(pfes->GetComm(), pfes->GlobalVSize(),
+                             pfes->GetDofOffsets(), Grad);
   }

   // RAP the local gradient dA.
@@ -271,7 +366,70 @@ void ParBlockNonlinearForm::Mult(const Vector &x, Vector &y) const

   if (fnfi.Size() > 0)
   {
-      MFEM_ABORT("TODO: assemble contributions from shared face terms");
+      // Terms over shared interior faces in parallel.
+      ParMesh *pmesh = ParFESpace(0)->GetParMesh();
+      FaceElementTransformations *tr;
+
+      Array<Array<int> *>vdofs(fes.Size());
+      Array<Array<int> *>vdofs2(fes.Size());
+      Array<Vector *> el_x(fes.Size());
+      Array<const Vector *> el_x_const(fes.Size());
+      Array<Vector *> el_y(fes.Size());
+      Array<const FiniteElement *> fe(fes.Size());
+      Array<const FiniteElement *> fe2(fes.Size());
+      Array<ParGridFunction *> pgfs(fes.Size());
+      for (int s=0; s<fes.Size(); ++s)
+      {
+         el_x_const[s] = el_x[s] = new Vector();
+         el_y[s] = new Vector();
+         vdofs[s] = new Array<int>;
+         vdofs2[s] = new Array<int>;
+         pgfs[s] = new ParGridFunction(const_cast<ParFiniteElementSpace*>(ParFESpace(s)),
+                                       xs.GetBlock(s));
+         pgfs[s]->ExchangeFaceNbrData();
+      }
+
+      const int n_shared_faces = pmesh->GetNSharedFaces();
+      for (int i = 0; i < n_shared_faces; i++)
+      {
+         tr = pmesh->GetSharedFaceTransformations(i, true);
+         int Elem2NbrNo = tr->Elem2No - pmesh->GetNE();
+
+         for (int s=0; s<fes.Size(); ++s)
+         {
+            const ParFiniteElementSpace *pfes = ParFESpace(s);
+            fe[s] = pfes->GetFE(tr->Elem1No);
+            fe2[s] = pfes->GetFaceNbrFE(Elem2NbrNo);
+
+            pfes->GetElementVDofs(tr->Elem1No, *(vdofs[s]));
+            pfes->GetFaceNbrElementVDofs(Elem2NbrNo, *(vdofs2[s]));
+
+            el_x[s]->SetSize(vdofs[s]->Size() + vdofs2[s]->Size());
+            xs.GetBlock(s).GetSubVector(*(vdofs[s]), el_x[s]->GetData());
+            pgfs[s]->FaceNbrData().GetSubVector(*(vdofs2[s]),
+                                                el_x[s]->GetData() + vdofs[s]->Size());
+         }
+
+         for (int k = 0; k < fnfi.Size(); ++k)
+         {
+            fnfi[k]->AssembleFaceVector(fe, fe2, *tr, el_x_const, el_y);
+
+            for (int s=0; s<fes.Size(); ++s)
+            {
+               if (el_y[s]->Size() == 0) { continue; }
+               ys.GetBlock(s).AddElementVector(*(vdofs[s]), *el_y[s]);
+            }
+         }
+      }
+
+      for (int s=0; s<fes.Size(); ++s)
+      {
+         delete pgfs[s];
+         delete vdofs2[s];
+         delete vdofs[s];
+         delete el_y[s];
+         delete el_x[s];
+      }
   }

   for (int s=0; s<fes.Size(); ++s)
@@ -328,6 +486,106 @@ void ParBlockNonlinearForm::SetGradientType(Operator::Type tid)
   }
 }

+void ParBlockNonlinearForm::GradientSharedFaces(const BlockVector &xs,
+                                                int skip_zeros) const
+{
+   // Terms over shared interior faces in parallel.
+   ParMesh *pmesh = ParFESpace(0)->GetParMesh();
+   FaceElementTransformations *tr;
+
+   Array<Array<int> *>vdofs(fes.Size());
+   Array<Array<int> *>vdofs2(fes.Size());
+   Array<Array<int> *>vdofs_all(fes.Size());
+   Array<Vector *> el_x(fes.Size());
+   Array<const Vector *> el_x_const(fes.Size());
+   Array2D<DenseMatrix *> elmats(fes.Size(), fes.Size());
+   Array<const FiniteElement *> fe(fes.Size());
+   Array<const FiniteElement *> fe2(fes.Size());
+   Array<ParGridFunction *> pgfs(fes.Size());
+
+   for (int s1=0; s1<fes.Size(); ++s1)
+   {
+      el_x_const[s1] = el_x[s1] = new Vector();
+      vdofs[s1] = new Array<int>;
+      vdofs2[s1] = new Array<int>;
+      vdofs_all[s1] = new Array<int>;
+      pgfs[s1] = new ParGridFunction(
+         const_cast<ParFiniteElementSpace*>(ParFESpace(s1)),
+         const_cast<Vector&>(xs.GetBlock(s1)));
+      pgfs[s1]->ExchangeFaceNbrData();
+      for (int s2=0; s2<fes.Size(); ++s2)
+      {
+         elmats(s1,s2) = new DenseMatrix();
+      }
+   }
+
+   const int n_shared_faces = pmesh->GetNSharedFaces();
+   for (int i = 0; i < n_shared_faces; i++)
+   {
+      tr = pmesh->GetSharedFaceTransformations(i, true);
+      int Elem2NbrNo = tr->Elem2No - pmesh->GetNE();
+
+      for (int s=0; s<fes.Size(); ++s)
+      {
+         const ParFiniteElementSpace *pfes = ParFESpace(s);
+         fe[s] = pfes->GetFE(tr->Elem1No);
+         fe2[s] = pfes->GetFaceNbrFE(Elem2NbrNo);
+
+         pfes->GetElementVDofs(tr->Elem1No, *(vdofs[s]));
+         pfes->GetFaceNbrElementVDofs(Elem2NbrNo, *(vdofs2[s]));
+
+         el_x[s]->SetSize(vdofs[s]->Size() + vdofs2[s]->Size());
+         xs.GetBlock(s).GetSubVector(*(vdofs[s]), el_x[s]->GetData());
+         pgfs[s]->FaceNbrData().GetSubVector(*(vdofs2[s]),
+                                             el_x[s]->GetData() + vdofs[s]->Size());
+
+         vdofs[s]->Copy(*vdofs_all[s]);
+
+         const int lvsize = pfes->GetVSize();
+         for (int j = 0; j < vdofs2[s]->Size(); j++)
+         {
+            if ((*vdofs2[s])[j] >= 0)
+            {
+               (*vdofs2[s])[j] += lvsize;
+            }
+            else
+            {
+               (*vdofs2[s])[j] -= lvsize;
+            }
+         }
+         vdofs_all[s]->Append(*(vdofs2[s]));
+      }
+
+      for (int k = 0; k < fnfi.Size(); ++k)
+      {
+         fnfi[k]->AssembleFaceGrad(fe, fe2, *tr, el_x_const, elmats);
+
+         for (int s1=0; s1<fes.Size(); ++s1)
+         {
+            for (int s2=0; s2<fes.Size(); ++s2)
+            {
+               if (elmats(s1,s2)->Height() == 0) { continue; }
+               Grads(s1,s2)->AddSubMatrix(*vdofs[s1], *vdofs_all[s2],
+                                          *elmats(s1,s2), skip_zeros);
+            }
+         }
+      }
+   }
+
+   for (int s1=0; s1<fes.Size(); ++s1)
+   {
+      delete pgfs[s1];
+      delete vdofs_all[s1];
+      delete vdofs2[s1];
+      delete vdofs[s1];
+      delete el_x[s1];
+      for (int s2=0; s2<fes.Size(); ++s2)
+      {
+         delete elmats(s1,s2);
+      }
+   }
+}
+
 BlockOperator & ParBlockNonlinearForm::GetGradient(const Vector &x) const
 {
   if (pBlockGrad == NULL)
@@ -347,49 +605,155 @@ BlockOperator & ParBlockNonlinearForm::GetGradient(const Vector &x) const
      }
   }

-   GetLocalGradient(x); // gradients are stored in 'Grads'
+   // xs_true is not modified, so const_cast is okay
+   xs_true.Update(const_cast<Vector &>(x), block_trueOffsets);
+   xs.Update(block_offsets);
+
+   for (int s=0; s<fes.Size(); ++s)
+   {
+      fes[s]->GetProlongationMatrix()->Mult(
+         xs_true.GetBlock(s), xs.GetBlock(s));
+   }

   if (fnfi.Size() > 0)
   {
-      MFEM_ABORT("TODO: assemble contributions from shared face terms");
-   }
+      const int skip_zeros = 0;

-   for (int s1=0; s1<fes.Size(); ++s1)
-   {
-      for (int s2=0; s2<fes.Size(); ++s2)
+      for (int s=0; s<fes.Size(); ++s)
      {
-         OperatorHandle dA(phBlockGrad(s1,s2)->Type()),
-                        Ph(phBlockGrad(s1,s2)->Type()),
-                        Rh(phBlockGrad(s1,s2)->Type());
+         const_cast<ParFiniteElementSpace*>(pfes[s])->ExchangeFaceNbrData();
+      }

-         if (s1 == s2)
+      for (int s1=0; s1<fes.Size(); ++s1)
+      {
+         for (int s2=0; s2<fes.Size(); ++s2)
         {
-            dA.MakeSquareBlockDiag(pfes[s1]->GetComm(), pfes[s1]->GlobalVSize(),
-                                   pfes[s1]->GetDofOffsets(), Grads(s1,s1));
-            Ph.ConvertFrom(pfes[s1]->Dof_TrueDof_Matrix());
-            phBlockGrad(s1,s1)->MakePtAP(dA, Ph);
-
-            OperatorHandle Ae;
-            Ae.EliminateRowsCols(*phBlockGrad(s1,s1), *ess_tdofs[s1]);
+            if (Grads(s1,s2) == NULL)
+            {
+               int nbr_size = pfes[s2]->GetFaceNbrVSize();
+               Grads(s1,s2) = new SparseMatrix(pfes[s1]->GetVSize(),
+                                               pfes[s2]->GetVSize() + nbr_size);
+            }
         }
-         else
+      }
+
+      // (re)assemble Grad without b.c. into 'Grads'
+      BlockNonlinearForm::ComputeGradientBlocked(xs, false);
+
+      GradientSharedFaces(xs, skip_zeros);
+
+      // finalize the gradients
+      for (int s1=0; s1<fes.Size(); ++s1)
+         for (int s2=0; s2<fes.Size(); ++s2)
         {
-            dA.MakeRectangularBlockDiag(pfes[s1]->GetComm(),
-                                        pfes[s1]->GlobalVSize(),
-                                        pfes[s2]->GlobalVSize(),
-                                        pfes[s1]->GetDofOffsets(),
-                                        pfes[s2]->GetDofOffsets(),
-                                        Grads(s1,s2));
-            Rh.ConvertFrom(pfes[s1]->Dof_TrueDof_Matrix());
-            Ph.ConvertFrom(pfes[s2]->Dof_TrueDof_Matrix());
-
-            phBlockGrad(s1,s2)->MakeRAP(Rh, dA, Ph);
-
-            phBlockGrad(s1,s2)->EliminateRows(*ess_tdofs[s1]);
-            phBlockGrad(s1,s2)->EliminateCols(*ess_tdofs[s2]);
+            Grads(s1,s2)->Finalize(skip_zeros);
         }

-         pBlockGrad->SetBlock(s1, s2, phBlockGrad(s1,s2)->Ptr());
+      for (int s1=0; s1<fes.Size(); ++s1)
+      {
+         for (int s2=0; s2<fes.Size(); ++s2)
+         {
+            OperatorHandle hdA;
+            OperatorHandle dA(phBlockGrad(s1,s2)->Type()),
+                           Ph(phBlockGrad(s1,s2)->Type()),
+                           Rh(phBlockGrad(s1,s2)->Type());
+
+            // handle the case when 'a' contains off-diagonal
+            int lvsize = pfes[s2]->GetVSize();
+            const HYPRE_BigInt *face_nbr_glob_ldof =
+               const_cast<ParFiniteElementSpace*>(pfes[s2])->GetFaceNbrGlobalDofMap();
+            HYPRE_BigInt ldof_offset = pfes[s2]->GetMyDofOffset();
+
+            Array<HYPRE_BigInt> glob_J(Grads(s1,s2)->NumNonZeroElems());
+            int *J = Grads(s1,s2)->GetJ();
+            for (int i = 0; i < glob_J.Size(); i++)
+            {
+               if (J[i] < lvsize)
+               {
+                  glob_J[i] = J[i] + ldof_offset;
+               }
+               else
+               {
+                  glob_J[i] = face_nbr_glob_ldof[J[i] - lvsize];
+               }
+            }
+
+            // TODO - construct dA directly in the A format
+            hdA.Reset(
+               new HypreParMatrix(pfes[s2]->GetComm(), pfes[s1]->GetVSize(),
+                                  pfes[s1]->GlobalVSize(), pfes[s2]->GlobalVSize(),
+                                  Grads(s1,s2)->GetI(), glob_J, Grads(s1,s2)->GetData(),
+                                  pfes[s1]->GetDofOffsets(), pfes[s2]->GetDofOffsets()));
+            // - hdA owns the new HypreParMatrix
+            // - the above constructor copies all input arrays
+            glob_J.DeleteAll();
+            dA.ConvertFrom(hdA);
+
+            if (s1 == s2)
+            {
+               Ph.ConvertFrom(pfes[s1]->Dof_TrueDof_Matrix());
+               phBlockGrad(s1,s1)->MakePtAP(dA, Ph);
+
+               OperatorHandle Ae;
+               Ae.EliminateRowsCols(*phBlockGrad(s1,s1), *ess_tdofs[s1]);
+            }
+            else
+            {
+               Rh.ConvertFrom(pfes[s1]->Dof_TrueDof_Matrix());
+               Ph.ConvertFrom(pfes[s2]->Dof_TrueDof_Matrix());
+
+               phBlockGrad(s1,s2)->MakeRAP(Rh, dA, Ph);
+
+               phBlockGrad(s1,s2)->EliminateRows(*ess_tdofs[s1]);
+               phBlockGrad(s1,s2)->EliminateCols(*ess_tdofs[s2]);
+            }
+
+            pBlockGrad->SetBlock(s1, s2, phBlockGrad(s1,s2)->Ptr());
+         }
+      }
+   }
+   else
+   {
+      // (re)assemble Grad without b.c. into 'Grads'
+      BlockNonlinearForm::ComputeGradientBlocked(xs);
+
+      for (int s1=0; s1<fes.Size(); ++s1)
+      {
+         for (int s2=0; s2<fes.Size(); ++s2)
+         {
+            OperatorHandle dA(phBlockGrad(s1,s2)->Type()),
+                           Ph(phBlockGrad(s1,s2)->Type()),
+                           Rh(phBlockGrad(s1,s2)->Type());
+
+            if (s1 == s2)
+            {
+               dA.MakeSquareBlockDiag(pfes[s1]->GetComm(), pfes[s1]->GlobalVSize(),
+                                      pfes[s1]->GetDofOffsets(), Grads(s1,s1));
+               Ph.ConvertFrom(pfes[s1]->Dof_TrueDof_Matrix());
+               phBlockGrad(s1,s1)->MakePtAP(dA, Ph);
+
+               OperatorHandle Ae;
+               Ae.EliminateRowsCols(*phBlockGrad(s1,s1), *ess_tdofs[s1]);
+            }
+            else
+            {
+               dA.MakeRectangularBlockDiag(pfes[s1]->GetComm(),
+                                           pfes[s1]->GlobalVSize(),
+                                           pfes[s2]->GlobalVSize(),
+                                           pfes[s1]->GetDofOffsets(),
+                                           pfes[s2]->GetDofOffsets(),
+                                           Grads(s1,s2));
+               Rh.ConvertFrom(pfes[s1]->Dof_TrueDof_Matrix());
+               Ph.ConvertFrom(pfes[s2]->Dof_TrueDof_Matrix());
+
+               phBlockGrad(s1,s2)->MakeRAP(Rh, dA, Ph);
+
+               phBlockGrad(s1,s2)->EliminateRows(*ess_tdofs[s1]);
+               phBlockGrad(s1,s2)->EliminateCols(*ess_tdofs[s2]);
+            }
+
+            pBlockGrad->SetBlock(s1, s2, phBlockGrad(s1,s2)->Ptr());
+         }
      }
   }

@@ -29,6 +29,8 @@ protected:
   mutable ParGridFunction X, Y;
   mutable OperatorHandle pGrad;

+   void GradientSharedFaces(const Vector &x, int skip_zeros = 1) const;
+
 public:
   ParNonlinearForm(ParFiniteElementSpace *pf);

@@ -81,6 +83,8 @@ protected:
   mutable Array2D<OperatorHandle *> phBlockGrad;
   mutable BlockOperator *pBlockGrad;

+   void GradientSharedFaces(const BlockVector &xs, int skip_zeros) const;
+
 public:
   /// Computes the energy of the system
   real_t GetEnergy(const Vector &x) const override;
@@ -1,214 +0,0 @@
-// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
-// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
-// LICENSE and NOTICE for details. LLNL-CODE-806117.
-//
-// This file is part of the MFEM library. For more information and source code
-// availability visit https://mfem.org.
-//
-// MFEM is free software; you can redistribute it and/or modify it under the
-// terms of the BSD-3 license. We welcome feedback and contributions, see file
-// CONTRIBUTING.md for details.
-
-
-// Abstract array data type
-
-#include "array.hpp"
-#include "../general/forall.hpp"
-#include <fstream>
-#include <type_traits>
-
-namespace mfem
-{
-
-template <class T>
-void Array<T>::Print(std::ostream &os, int width) const
-{
-   for (int i = 0; i < size; i++)
-   {
-      os << data[i];
-      if ( !((i+1) % width) || i+1 == size )
-      {
-         os << '\n';
-      }
-      else
-      {
-         os << " ";
-      }
-   }
-}
-
-template <class T>
-void Array<T>::Save(std::ostream &os, int fmt) const
-{
-   if (fmt == 0)
-   {
-      os << size << '\n';
-   }
-   for (int i = 0; i < size; i++)
-   {
-      os << operator[](i) << '\n';
-   }
-}
-
-template <class T>
-void Array<T>::Load(std::istream &in, int fmt)
-{
-   if (fmt == 0)
-   {
-      int new_size;
-      in >> new_size;
-      SetSize(new_size);
-   }
-   for (int i = 0; i < size; i++)
-   {
-      in >> operator[](i);
-   }
-}
-
-template <class T>
-T Array<T>::Max() const
-{
-   MFEM_ASSERT(size > 0, "Array is empty with size " << size);
-
-   T max = operator[](0);
-   for (int i = 1; i < size; i++)
-   {
-      if (max < operator[](i))
-      {
-         max = operator[](i);
-      }
-   }
-
-   return max;
-}
-
-template <class T>
-T Array<T>::Min() const
-{
-   MFEM_ASSERT(size > 0, "Array is empty with size " << size);
-
-   T min = operator[](0);
-   for (int i = 1; i < size; i++)
-   {
-      if (operator[](i) < min)
-      {
-         min = operator[](i);
-      }
-   }
-
-   return min;
-}
-
-// Partial Sum
-template <class T>
-void Array<T>::PartialSum()
-{
-   T sum = static_cast<T>(0);
-   for (int i = 0; i < size; i++)
-   {
-      sum+=operator[](i);
-      operator[](i) = sum;
-   }
-}
-
-template <class T>
-void Array<T>::Abs()
-{
-   static_assert(std::is_arithmetic<T>::value, "Use with arithmetic types!");
-   const bool useDevice = UseDevice();
-   const int N = size;
-   auto y = ReadWrite(useDevice);
-   mfem::forall_switch(useDevice, N, [=] MFEM_HOST_DEVICE (int i)
-   {
-      y[i] = std::abs(y[i]);
-   });
-}
-
-// Sum
-template <class T>
-T Array<T>::Sum() const
-{
-   T sum = static_cast<T>(0);
-   for (int i = 0; i < size; i++)
-   {
-      sum+=operator[](i);
-   }
-
-   return sum;
-}
-
-template <class T>
-int Array<T>::IsSorted() const
-{
-   T val_prev = operator[](0), val;
-   for (int i = 1; i < size; i++)
-   {
-      val=operator[](i);
-      if (val < val_prev)
-      {
-         return 0;
-      }
-      val_prev = val;
-   }
-
-   return 1;
-}
-
-template <class T>
-bool Array<T>::IsConstant() const
-{
-   if (size < 2) { return true; }
-   const T v0 = data[0];
-   for (int i = 1; i < size; i++)
-   {
-      if (data[i] != v0)
-      {
-         return false;
-      }
-   }
-
-   return true;
-}
-
-template <class T>
-void Array2D<T>::Load(const char *filename, int fmt)
-{
-   std::ifstream in;
-   in.open(filename, std::ifstream::in);
-   MFEM_VERIFY(in.is_open(), "File " << filename << " does not exist.");
-   Load(in, fmt);
-   in.close();
-}
-
-template <class T>
-void Array2D<T>::Print(std::ostream &os, int width_)
-{
-   int height = this->NumRows();
-   int width  = this->NumCols();
-
-   for (int i = 0; i < height; i++)
-   {
-      os << "[row " << i << "]\n";
-      for (int j = 0; j < width; j++)
-      {
-         os << (*this)(i,j);
-         if ( (j+1) == width_ || (j+1) % width_ == 0 )
-         {
-            os << '\n';
-         }
-         else
-         {
-            os << ' ';
-         }
-      }
-   }
-}
-
-template class Array<char>;
-template class Array<int>;
-template class Array<long long>;
-template class Array<real_t>;
-template class Array2D<int>;
-template class Array2D<real_t>;
-
-} // namespace mfem
@@ -16,9 +16,13 @@
 #include "mem_manager.hpp"
 #include "device.hpp"
 #include "error.hpp"
+#include "forall.hpp"
 #include "globals.hpp"
+#include "reducers.hpp"
+#include "scan.hpp"

 #include <iostream>
+#include <fstream>
 #include <cstdlib>
 #include <cstring>
 #include <algorithm>
@@ -135,6 +139,8 @@ public:
   /// Return the device flag of the Memory object used by the Array
   bool UseDevice() const { return data.UseDevice(); }

+   void UseDevice(bool use_dev) { data.UseDevice(use_dev); }
+
   /// Return true if the data will be deleted by the Array
   inline bool OwnsData() const { return data.OwnsHostPtr(); }

@@ -275,11 +281,11 @@ public:

   /** @brief Find the maximal element in the array, using the comparison
       operator `<` for class T. */
-   T Max() const;
+   inline T Max() const;

   /** @brief Find the minimal element in the array, using the comparison
       operator `<` for class T. */
-   T Min() const;
+   inline T Min() const;

   /// Sorts the array in ascending order. This requires operator< to be defined for T.
   void Sort() { std::sort((T*)data, data + size); }
@@ -297,22 +303,22 @@ public:
   }

   /// Return 1 if the array is sorted from lowest to highest.  Otherwise return 0.
-   int IsSorted() const;
+   inline int IsSorted() const;

   /// Does the Array have Size zero.
   bool IsEmpty() const { return Size() == 0; }

   /// Return true if all entries of the array are the same.
-   bool IsConstant() const;
+   inline bool IsConstant() const;

   /// Fill the entries of the array with the cumulative sum of the entries.
-   void PartialSum();
+   inline void PartialSum();

   /// Replace each entry of the array with its absolute value.
-   void Abs();
+   inline void Abs();

   /// Return the sum of all the array entries using the '+'' operator for class 'T'.
-   T Sum() const;
+   inline T Sum() const;

   /// Set all entries of the array to the provided constant.
   inline void operator=(const T &a);
@@ -797,8 +803,14 @@ template <typename T> template <typename CT>
 inline Array<T> &Array<T>::operator=(const Array<CT> &src)
 {
   SetSize(src.Size());
-   for (int i = 0; i < size; i++) { (*this)[i] = T(src[i]); }
-   return *this;
+
+   const bool use_dev = UseDevice() || src.UseDevice();
+   const auto x = src.Read(use_dev);
+   auto y = Write(use_dev);
+   mfem::forall_switch(use_dev, size, [=] MFEM_HOST_DEVICE (int i)
+   {
+      y[i] = x[i];
+   });
 }

 template <class T>
@@ -1014,19 +1026,24 @@ template <class T>
 inline void Array<T>::GetSubArray(int offset, int sa_size, Array<T> &sa) const
 {
   sa.SetSize(sa_size);
-   for (int i = 0; i < sa_size; i++)
+   const bool use_dev = UseDevice() || sa.UseDevice();
+   const auto x = Read(use_dev);
+   auto y = sa.Write(use_dev);
+   mfem::forall_switch(use_dev, sa_size, [=] MFEM_HOST_DEVICE (int i)
   {
-      sa[i] = (*this)[offset+i];
-   }
+      y[i] = x[offset + i];
+   });
 }

 template <class T>
 inline void Array<T>::operator=(const T &a)
 {
-   for (int i = 0; i < size; i++)
+   const bool use_dev = UseDevice();
+   auto x = Write(use_dev);
+   mfem::forall_switch(use_dev, size, [=] MFEM_HOST_DEVICE (int i)
   {
-      data[i] = a;
-   }
+      x[i] = a;
+   });
 }

 template <class T>
@@ -1035,6 +1052,153 @@ inline void Array<T>::Assign(const T *p)
   data.CopyFromHost(p, Size());
 }

+template <class T>
+inline void Array<T>::Print(std::ostream &os, int width) const
+{
+   for (int i = 0; i < size; i++)
+   {
+      os << data[i];
+      if ( !((i+1) % width) || i+1 == size )
+      {
+         os << '\n';
+      }
+      else
+      {
+         os << " ";
+      }
+   }
+}
+
+template <class T>
+inline void Array<T>::Save(std::ostream &os, int fmt) const
+{
+   if (fmt == 0)
+   {
+      os << size << '\n';
+   }
+   for (int i = 0; i < size; i++)
+   {
+      os << operator[](i) << '\n';
+   }
+}
+
+template <class T>
+void Array<T>::Load(std::istream &in, int fmt)
+{
+   if (fmt == 0)
+   {
+      int new_size;
+      in >> new_size;
+      SetSize(new_size);
+   }
+   for (int i = 0; i < size; i++)
+   {
+      in >> operator[](i);
+   }
+}
+
+template <class T>
+inline T Array<T>::Max() const
+{
+   MFEM_ASSERT(size > 0, "Array is empty with size " << size);
+
+   T max = operator[](0);
+   for (int i = 1; i < size; i++)
+   {
+      if (max < operator[](i))
+      {
+         max = operator[](i);
+      }
+   }
+
+   return max;
+}
+
+template <class T>
+inline T Array<T>::Min() const
+{
+   MFEM_ASSERT(size > 0, "Array is empty with size " << size);
+
+   T min = operator[](0);
+   for (int i = 1; i < size; i++)
+   {
+      if (operator[](i) < min)
+      {
+         min = operator[](i);
+      }
+   }
+
+   return min;
+}
+
+// Partial Sum
+template <class T>
+inline void Array<T>::PartialSum()
+{
+   auto data_ptr = ReadWrite(UseDevice());
+   InclusiveScan(UseDevice(), data_ptr, data_ptr, size);
+}
+
+template <class T>
+inline void Array<T>::Abs()
+{
+   static_assert(std::is_arithmetic<T>::value, "Use with arithmetic types!");
+   const bool useDevice = UseDevice();
+   const int N = size;
+   auto y = ReadWrite(useDevice);
+   mfem::forall_switch(useDevice, N, [=] MFEM_HOST_DEVICE (int i)
+   {
+      y[i] = std::abs(y[i]);
+   });
+}
+
+// Sum
+template <class T>
+inline T Array<T>::Sum() const
+{
+   T sum = static_cast<T>(0);
+   if (size > 0)
+   {
+      const auto m_data = Read(UseDevice());
+      reduce(size, sum, [=] MFEM_HOST_DEVICE(int i, T &r) { r += m_data[i]; },
+      /*  */ SumReducer<T> {}, UseDevice());
+   }
+   return sum;
+}
+
+template <class T>
+inline int Array<T>::IsSorted() const
+{
+   T val_prev = operator[](0), val;
+   for (int i = 1; i < size; i++)
+   {
+      val=operator[](i);
+      if (val < val_prev)
+      {
+         return 0;
+      }
+      val_prev = val;
+   }
+
+   return 1;
+}
+
+template <class T>
+inline bool Array<T>::IsConstant() const
+{
+   if (size < 2) { return true; }
+   const T v0 = data[0];
+   for (int i = 1; i < size; i++)
+   {
+      if (data[i] != v0)
+      {
+         return false;
+      }
+   }
+
+   return true;
+}
+

 template <class T>
 inline const T &Array2D<T>::operator()(int i, int j) const
@@ -1074,6 +1238,40 @@ inline T *Array2D<T>::operator[](int i)
   return &array1d[i*N];
 }

+template <class T>
+void Array2D<T>::Load(const char *filename, int fmt)
+{
+   std::ifstream in;
+   in.open(filename, std::ifstream::in);
+   MFEM_VERIFY(in.is_open(), "File " << filename << " does not exist.");
+   Load(in, fmt);
+   in.close();
+}
+
+template <class T>
+void Array2D<T>::Print(std::ostream &os, int width_)
+{
+   int height = this->NumRows();
+   int width  = this->NumCols();
+
+   for (int i = 0; i < height; i++)
+   {
+      os << "[row " << i << "]\n";
+      for (int j = 0; j < width; j++)
+      {
+         os << (*this)(i,j);
+         if ( (j+1) == width_ || (j+1) % width_ == 0 )
+         {
+            os << '\n';
+         }
+         else
+         {
+            os << ' ';
+         }
+      }
+   }
+}
+

 template <class T>
 inline void Swap(Array2D<T> &a, Array2D<T> &b)
@@ -12,7 +12,6 @@
 #ifndef MFEM_REDUCERS_HPP
 #define MFEM_REDUCERS_HPP

-#include "array.hpp"
 #include "forall.hpp"

 #include <cmath>
@@ -514,6 +513,33 @@ template<class B, class R> struct reduction_kernel
      }
   }
 };
+
+template <class T>
+class ReductionWorkspace
+{
+   Memory<T> workspace;
+
+   static ReductionWorkspace &Instance()
+   {
+      static ReductionWorkspace instance;
+      return instance;
+   }
+
+   ~ReductionWorkspace() { workspace.Delete(); }
+
+public:
+   static T *Get(int num_blocks)
+   {
+      ReductionWorkspace &instance = Instance();
+      if (instance.workspace.Capacity() < num_blocks)
+      {
+         instance.workspace.Delete();
+         instance.workspace.New(num_blocks, MemoryType::HOST_PINNED);
+      }
+      return instance.workspace;
+   }
+};
+
 }

 /**
@@ -529,8 +555,7 @@ template<class B, class R> struct reduction_kernel
 @tparam T value_type to operate on
 */
 template <class T, class B, class R>
-void reduce(int N, T &res, B &&body, const R &reducer, bool use_dev,
-            Array<T> &workspace)
+void reduce(int N, T &res, B &&body, const R &reducer, bool use_dev)
 {
   if (N == 0)
   {
@@ -567,13 +592,7 @@ void reduce(int N, T &res, B &&body, const R &reducer, bool use_dev,

      red_type red{nullptr, std::forward<B>(body), reducer, N, items_per_thread};
      // allocate res to fit block_size entries
-      auto mt = workspace.GetMemory().GetMemoryType();
-      if (mt != MemoryType::HOST_PINNED && mt != MemoryType::MANAGED)
-      {
-         mt = MemoryType::HOST_PINNED;
-      }
-      workspace.SetSize(nblocks, mt);
-      auto work = workspace.HostWrite();
+      auto work = internal::ReductionWorkspace<T>::Get(nblocks);
      red.work = work;
      forall_2D(nblocks, block_size, 1, std::move(red));
      // wait for results
@@ -28,8 +28,37 @@

 namespace mfem
 {
-/// Equivalent to InclusiveScan(use_dev, d_in, d_out, num_items, workspace,
-/// std::plus<>{})
+
+namespace internal
+{
+class ScanWorkspace
+{
+   Memory<std::byte> workspace;
+   static ScanWorkspace &Instance()
+   {
+      static ScanWorkspace instance;
+      return instance;
+   }
+   ~ScanWorkspace() { workspace.Delete(); }
+public:
+   static std::byte *Get(int num_bytes)
+   {
+      ScanWorkspace &instance = Instance();
+      if (Size() < num_bytes)
+      {
+         instance.workspace.Delete();
+         instance.workspace.New(num_bytes);
+      }
+      return instance.workspace.Write(MemoryClass::DEVICE, Size());
+   }
+   static int Size()
+   {
+      return Instance().workspace.Capacity();
+   }
+};
+}
+
+/// Equivalent to InclusiveScan(use_dev, d_in, d_out, num_items, std::plus<>{})
 template <class InputIt, class OutputIt>
 void InclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items)
 {
@@ -37,12 +66,12 @@ void InclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items)
 #if defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP)
   if (use_dev && mfem::Device::Allows(Backend::CUDA_MASK | Backend::HIP_MASK))
   {
-      static Array<std::byte> workspace;
-      size_t bytes = workspace.Size();
-      if (bytes)
+      using internal::ScanWorkspace;
+      size_t bytes = ScanWorkspace::Size();
+      if (bytes > 0)
      {
         auto err = MFEM_CUB_NAMESPACE::DeviceScan::InclusiveSum(
-                       workspace.Write(), bytes, d_in, d_out, num_items);
+                       ScanWorkspace::Get(bytes), bytes, d_in, d_out, num_items);
 #if defined(MFEM_USE_CUDA)
         if (err == cudaSuccess)
         {
@@ -57,11 +86,12 @@ void InclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items)
      }
      // try allocating a larger buffer
      bytes = 0;
+      // get size of buffer
      MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::InclusiveSum(
                        nullptr, bytes, d_in, d_out, num_items));
-      workspace.SetSize(bytes);
+      // resize buffer (in ScanWorkspace::Get) and try again
      MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::InclusiveSum(
-                        workspace.Write(), bytes, d_in, d_out, num_items));
+                        ScanWorkspace::Get(bytes), bytes, d_in, d_out, num_items));
      return;
   }
 #endif
@@ -101,12 +131,13 @@ void InclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
 #if defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP)
   if (use_dev && mfem::Device::Allows(Backend::CUDA_MASK | Backend::HIP_MASK))
   {
-      static Array<std::byte> workspace;
-      size_t bytes = workspace.Size();
-      if (bytes)
+      using internal::ScanWorkspace;
+      size_t bytes = ScanWorkspace::Size();
+      if (bytes > 0)
      {
         auto err = MFEM_CUB_NAMESPACE::DeviceScan::InclusiveScan(
-                       workspace.Write(), bytes, d_in, d_out, scan_op, num_items);
+                       ScanWorkspace::Get(bytes), bytes, d_in, d_out, scan_op,
+                       num_items);
 #if defined(MFEM_USE_CUDA)
         if (err == cudaSuccess)
         {
@@ -123,9 +154,9 @@ void InclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
      bytes = 0;
      MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::InclusiveScan(
                        nullptr, bytes, d_in, d_out, scan_op, num_items));
-      workspace.SetSize(bytes);
      MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::InclusiveScan(
-                        workspace.Write(), bytes, d_in, d_out, scan_op, num_items));
+                        ScanWorkspace::Get(bytes), bytes, d_in, d_out, scan_op,
+                        num_items));
      return;
   }
 #endif
@@ -164,13 +195,13 @@ void ExclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
 #if defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP)
   if (use_dev && mfem::Device::Allows(Backend::CUDA_MASK | Backend::HIP_MASK))
   {
-      static Array<std::byte> workspace;
-      size_t bytes = workspace.Size();
+      using internal::ScanWorkspace;
+      size_t bytes = ScanWorkspace::Size();
      if (bytes)
      {
         auto err = MFEM_CUB_NAMESPACE::DeviceScan::ExclusiveScan(
-                       workspace.Write(), bytes, d_in, d_out, scan_op, init_value,
-                       num_items);
+                       ScanWorkspace::Get(bytes), bytes, d_in, d_out, scan_op,
+                       init_value, num_items);
 #if defined(MFEM_USE_CUDA)
         if (err == cudaSuccess)
         {
@@ -187,10 +218,9 @@ void ExclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
      bytes = 0;
      MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::ExclusiveScan(
                        nullptr, bytes, d_in, d_out, scan_op, init_value, num_items));
-      workspace.SetSize(bytes);
      MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::ExclusiveScan(
-                        workspace.Write(), bytes, d_in, d_out, scan_op, init_value,
-                        num_items));
+                        ScanWorkspace::Get(bytes), bytes, d_in, d_out, scan_op,
+                        init_value, num_items));
      return;
   }
 #endif
@@ -213,7 +243,7 @@ void ExclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
 }

 /// Equivalent to ExclusiveScan(use_dev, d_in, d_out, num_items, init_value,
-/// workspace, std::plus<>{})
+/// std::plus<>{})
 template <class InputIt, class OutputIt, class T>
 void ExclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
                   T init_value)
@@ -167,7 +167,7 @@ void MagmaBatchedLinAlg::Invert(DenseTensor &A) const
   magma_int_t status;

   status = MFEM_MAGMA_PREFIX(getrf_batched)(
-               n, n, d_A_ptrs, n, d_P_ptrs, info_array.Write(), n_mat,
+               n, n, d_LU_ptrs, n, d_P_ptrs, info_array.Write(), n_mat,
               Magma::Queue());
   MFEM_VERIFY(status == MAGMA_SUCCESS, "");

@@ -561,7 +561,8 @@ void CopyMemory(Memory<T> &src, Memory<T> &dst, MemoryClass dst_mc,
    this function. In particular, @a dst should be empty or deleted before
    calling this function. */
 template <typename SrcT, typename DstT>
-void CopyConvertMemory(Memory<SrcT> &src, MemoryClass dst_mc, Memory<DstT> &dst)
+void CopyConvertMemory(const Memory<SrcT> &src, MemoryClass dst_mc,
+                       Memory<DstT> &dst)
 {
   auto capacity = src.Capacity();
   dst.New(capacity, GetMemoryType(dst_mc));
@@ -842,8 +843,8 @@ static int GetPartitioningArraySize(MPI_Comm comm)
 ///
 /// Both @a row and @a col are partitioning arrays, whose length is returned by
 /// GetPartitioningArraySize(), see @ref hypre_partitioning_descr.
-static bool RowAndColStartsAreEqual(MPI_Comm comm, HYPRE_BigInt *rows,
-                                    HYPRE_BigInt *cols)
+static bool RowAndColStartsAreEqual(MPI_Comm comm, const HYPRE_BigInt *rows,
+                                    const HYPRE_BigInt *cols)
 {
   const int part_size = GetPartitioningArraySize(comm);
   bool are_equal = true;
@@ -1131,7 +1132,7 @@ HypreParMatrix::HypreParMatrix(
 HypreParMatrix::HypreParMatrix(MPI_Comm comm,
                               HYPRE_BigInt *row_starts,
                               HYPRE_BigInt *col_starts,
-                               SparseMatrix *sm_a)
+                               const SparseMatrix *sm_a)
 {
   MFEM_ASSERT(sm_a != NULL, "invalid input");
   MFEM_VERIFY(!HYPRE_AssumedPartitionCheck(),
@@ -1145,7 +1146,7 @@ HypreParMatrix::HypreParMatrix(MPI_Comm comm,

   hypre_CSRMatrixSetDataOwner(csr_a,0);
   MemoryIJData mem_a;
-   CopyCSR(sm_a, mem_a, csr_a, false);
+   CopyCSR(const_cast<SparseMatrix*>(sm_a), mem_a, csr_a, false);
   hypre_CSRMatrixSetRownnz(csr_a);

   // NOTE: this call creates a matrix on host even when device support is
@@ -1307,10 +1308,11 @@ HypreParMatrix::HypreParMatrix(MPI_Comm comm, int id, int np,
 HypreParMatrix::HypreParMatrix(MPI_Comm comm, int nrows,
                               HYPRE_BigInt glob_nrows,
                               HYPRE_BigInt glob_ncols,
-                               int *I, HYPRE_BigInt *J,
-                               real_t *data,
-                               HYPRE_BigInt *rows,
-                               HYPRE_BigInt *cols)
+                               const int *I,
+                               const HYPRE_BigInt *J,
+                               const real_t *data,
+                               const HYPRE_BigInt *rows,
+                               const HYPRE_BigInt *cols)
 {
   Init();

@@ -565,7 +565,7 @@ public:
       partitioning arrays @a row_starts and @a col_starts. */
   HypreParMatrix(MPI_Comm comm, HYPRE_BigInt *row_starts,
                  HYPRE_BigInt *col_starts,
-                  SparseMatrix *a); // constructor with 4 arguments, v2
+                  const SparseMatrix *a); // constructor with 4 arguments, v2

   /// Creates boolean block-diagonal rectangular parallel matrix.
   /** The new HypreParMatrix does not take ownership of any of the input
@@ -594,9 +594,9 @@ public:
       arrays (so they can be deleted). See @ref hypre_partitioning_descr "here"
       for a description of the partitioning arrays @a rows and @a cols. */
   HypreParMatrix(MPI_Comm comm, int nrows, HYPRE_BigInt glob_nrows,
-                  HYPRE_BigInt glob_ncols, int *I, HYPRE_BigInt *J,
-                  real_t *data, HYPRE_BigInt *rows,
-                  HYPRE_BigInt *cols); // constructor with 9 arguments
+                  HYPRE_BigInt glob_ncols, const int *I, const HYPRE_BigInt *J,
+                  const real_t *data, const HYPRE_BigInt *rows,
+                  const HYPRE_BigInt *cols); // constructor with 9 arguments

   /** @brief Copy constructor for a ParCSR matrix which creates a deep copy of
       structure and data from @a P. */
@@ -46,6 +46,12 @@
 #define MFEM_GPUSPARSE_ALG HIPSPARSE_CSRMV_ALG1
 #endif // defined(MFEM_USE_CUDA)

+#if defined(MFEM_USE_SINGLE)
+#define MFEM_CUDA_or_HIP_REAL_T MFEM_CUDA_or_HIP(_R_32F)
+#elif defined(MFEM_USE_DOUBLE)
+#define MFEM_CUDA_or_HIP_REAL_T MFEM_CUDA_or_HIP(_R_64F)
+#endif
+
 namespace mfem
 {

@@ -57,8 +63,10 @@ int SparseMatrix::SparseMatrixCount = 0;
 /// @cond Suppress_Doxygen_warnings
 MFEM_cu_or_hip(sparseHandle_t) SparseMatrix::handle = nullptr;
 /// @endcond
+#ifndef MFEM_CUDA_1897_WORKAROUND
 size_t SparseMatrix::bufferSize = 0;
 void * SparseMatrix::dBuffer = nullptr;
+#endif
 #endif // MFEM_USE_CUDA_OR_HIP

 void SparseMatrix::InitGPUSparse()
@@ -464,109 +472,67 @@ void SparseMatrix::SortColumnIndices()
   }

 #ifdef MFEM_USE_CUDA_OR_HIP
-   if ( Device::Allows( Backend::CUDA_MASK ))
+   if (Device::Allows(Backend::CUDA_MASK) || Device::Allows(Backend::HIP_MASK))
   {
-#if defined(MFEM_USE_CUDA)
-      size_t pBufferSizeInBytes = 0;
-      void *pBuffer = NULL;
-
-      const int n = Height();
-      const int m = Width();
+      const int m = Height();
+      const int n = Width();
      const int nnzA = J.Capacity();
-      real_t * d_a_sorted = ReadWriteData();
-      const int * d_ia = ReadI();
-      int * d_ja_sorted = ReadWriteJ();
-      csru2csrInfo_t sortInfoA;
+      const int *d_ia = ReadI();
+      int *d_ja = ReadWriteJ();

-      cusparseMatDescr_t matA_descr;
-      cusparseCreateMatDescr( &matA_descr );
-      cusparseSetMatIndexBase( matA_descr, CUSPARSE_INDEX_BASE_ZERO );
-      cusparseSetMatType( matA_descr, CUSPARSE_MATRIX_TYPE_GENERAL );
+      // Get size of temporary buffer needed to sort the column indices,
+      // allocate the temporary buffer.
+      size_t pBufferSizeInBytes;
+      MFEM_cu_or_hip(sparseXcsrsort_bufferSizeExt)(handle, m, n, nnzA, d_ia,
+                                                   d_ja, &pBufferSizeInBytes);
+      void *pBuffer = MFEM_Cu_or_Hip(MemAlloc)(&pBuffer, pBufferSizeInBytes);

-      cusparseCreateCsru2csrInfo( &sortInfoA );
+      // Create matrix descriptor, will have default values
+      // CUSPARSE_INDEX_BASE_ZERO and CUSPARSE_MATRIX_TYPE_GENERAL.
+      MFEM_cu_or_hip(sparseMatDescr_t) matA_descr;
+      MFEM_cu_or_hip(sparseCreateMatDescr)(&matA_descr);

-#ifdef MFEM_USE_SINGLE
-      cusparseScsru2csr_bufferSizeExt( handle, n, m, nnzA, d_a_sorted, d_ia,
-                                       d_ja_sorted, sortInfoA,
-                                       &pBufferSizeInBytes);
-#elif defined MFEM_USE_DOUBLE
-      cusparseDcsru2csr_bufferSizeExt( handle, n, m, nnzA, d_a_sorted, d_ia,
-                                       d_ja_sorted, sortInfoA,
-                                       &pBufferSizeInBytes);
-#else
-      MFEM_ABORT("Floating point type undefined");
-#endif
+      // Initialize permutation to identity
+      Array<int> P(nnzA);
+      int *d_P = P.Write();
+      mfem::forall(nnzA, [=] MFEM_HOST_DEVICE (int i) { d_P[i] = i; });

-      CuMemAlloc( &pBuffer, pBufferSizeInBytes );
+      // Sort the column indices. The array d_ja will now be sorted. The
+      // permutation required to sort the values will be returned in d_P.
+      MFEM_cu_or_hip(sparseXcsrsort)(handle, m, n, nnzA, matA_descr, d_ia, d_ja,
+                                     d_P, pBuffer);

-#ifdef MFEM_USE_SINGLE
-      cusparseScsru2csr( handle, n, m, nnzA, matA_descr, d_a_sorted, d_ia,
-                         d_ja_sorted, sortInfoA, pBuffer);
-#elif defined MFEM_USE_DOUBLE
-      cusparseDcsru2csr( handle, n, m, nnzA, matA_descr, d_a_sorted, d_ia,
-                         d_ja_sorted, sortInfoA, pBuffer);
-#else
-      MFEM_ABORT("Floating point type undefined");
-#endif
+      // Create a copy of the unsorted matrix values.
+      real_t *d_a = ReadWriteData();
+      void *d_a_unsorted = MFEM_Cu_or_Hip(MemAlloc)(&d_a_unsorted,
+                                                    nnzA * sizeof(real_t));
+      MFEM_Cu_or_Hip(MemcpyDtoD)(d_a_unsorted, d_a, nnzA * sizeof(real_t));

-      // The above call is (at least in some cases) asynchronous, so we need to
-      // wait for it to finish before we can free device temporaries.
+      // Create the (input) dense vector with the unsorted values.
+      MFEM_cu_or_hip(sparseDnVecDescr_t) d_a_dense;
+      MFEM_cu_or_hip(sparseCreateDnVec)(&d_a_dense, nnzA, d_a_unsorted,
+                                        MFEM_CUDA_or_HIP_REAL_T);
+
+      // Create the (output) sparse vector that will have the sorted values.
+      MFEM_cu_or_hip(sparseSpVecDescr_t) d_a_sparse;
+      MFEM_cu_or_hip(sparseCreateSpVec)(&d_a_sparse, nnzA, nnzA, d_P, d_a,
+                                        MFEM_CU_or_HIP(SPARSE_INDEX_32I),
+                                        MFEM_CU_or_HIP(SPARSE_INDEX_BASE_ZERO),
+                                        MFEM_CUDA_or_HIP_REAL_T);
+
+      // Sort the matrix values using the permutation vector.
+      MFEM_cu_or_hip(sparseGather)(handle, d_a_dense, d_a_sparse);
+
+      // The above calls may be asynchronous, so we need to wait for them to
+      // finish before we can free memory.
      MFEM_STREAM_SYNC;

-      cusparseDestroyCsru2csrInfo( sortInfoA );
-      cusparseDestroyMatDescr( matA_descr );
+      MFEM_cu_or_hip(sparseDestroyDnVec)(d_a_dense);
+      MFEM_cu_or_hip(sparseDestroySpVec)(d_a_sparse);
+      MFEM_cu_or_hip(sparseDestroyMatDescr)(matA_descr);

-      CuMemFree( pBuffer );
-#endif
-   }
-   else if ( Device::Allows( Backend::HIP_MASK ))
-   {
-#if defined(MFEM_USE_HIP)
-      size_t pBufferSizeInBytes = 0;
-      void *pBuffer = NULL;
-      int *P = NULL;
-
-      const int n = Height();
-      const int m = Width();
-      const int nnzA = J.Capacity();
-      real_t * d_a_sorted = ReadWriteData();
-      const int * d_ia = ReadI();
-      int * d_ja_sorted = ReadWriteJ();
-
-      hipsparseMatDescr_t descrA;
-      hipsparseCreateMatDescr( &descrA );
-      // FIXME: There is not in-place version of csr sort in hipSPARSE currently, so we make
-      //        a temporary copy of the data for gthr, sort that, and then copy the sorted values
-      //        back to the array being returned. Where there is an in-place version available,
-      //        we should use it.
-      Array< real_t > a_tmp( nnzA );
-      real_t *d_a_tmp = a_tmp.Write();
-
-      hipsparseXcsrsort_bufferSizeExt(handle, n, m, nnzA, d_ia, d_ja_sorted,
-                                      &pBufferSizeInBytes);
-
-      HipMemAlloc( &pBuffer, pBufferSizeInBytes );
-      HipMemAlloc( (void**)&P, nnzA * sizeof(int) );
-
-      hipsparseCreateIdentityPermutation(handle, nnzA, P);
-      hipsparseXcsrsort(handle, n, m, nnzA, descrA, d_ia, d_ja_sorted, P, pBuffer);
-
-#if defined(MFEM_USE_SINGLE)
-      hipsparseSgthr(handle, nnzA, d_a_sorted, d_a_tmp, P,
-                     HIPSPARSE_INDEX_BASE_ZERO);
-#elif defined(MFEM_USE_DOUBLE)
-      hipsparseDgthr(handle, nnzA, d_a_sorted, d_a_tmp, P,
-                     HIPSPARSE_INDEX_BASE_ZERO);
-#else
-      MFEM_ABORT("Unsupported floating point type!");
-#endif
-
-      A.CopyFrom( a_tmp.GetMemory(), nnzA );
-      hipsparseDestroyMatDescr( descrA );
-
-      HipMemFree( pBuffer );
-      HipMemFree( P );
-#endif
+      MFEM_Cu_or_Hip(MemFree)(d_a_unsorted);
+      MFEM_Cu_or_Hip(MemFree)(pBuffer);
   }
   else
 #endif // MFEM_USE_CUDA_OR_HIP
@@ -821,27 +787,15 @@ void SparseMatrix::AddMult(const Vector &x, Vector &y, const real_t a) const
            MFEM_CU_or_HIP(SPARSE_INDEX_32I),
            MFEM_CU_or_HIP(SPARSE_INDEX_32I),
            MFEM_CU_or_HIP(SPARSE_INDEX_BASE_ZERO),
-#ifdef MFEM_USE_SINGLE
-            MFEM_CUDA_or_HIP(_R_32F));
-#else
-            MFEM_CUDA_or_HIP(_R_64F));
-#endif
+            MFEM_CUDA_or_HIP_REAL_T);

         // Create handles for input/output vectors
         MFEM_cu_or_hip(sparseCreateDnVec)(&vecX_descr,
                                           x.Size(),
                                           const_cast<real_t *>(d_x),
-#ifdef MFEM_USE_SINGLE
-                                           MFEM_CUDA_or_HIP(_R_32F));
-#else
-                                           MFEM_CUDA_or_HIP(_R_64F));
-#endif
+                                           MFEM_CUDA_or_HIP_REAL_T);
         MFEM_cu_or_hip(sparseCreateDnVec)(&vecY_descr, y.Size(), d_y,
-#ifdef MFEM_USE_SINGLE
-                                           MFEM_CUDA_or_HIP(_R_32F));
-#else
-                                           MFEM_CUDA_or_HIP(_R_64F));
-#endif
+                                           MFEM_CUDA_or_HIP_REAL_T);
 #else
         cusparseCreateMatDescr(&matA_descr);
         cusparseSetMatIndexBase(matA_descr, CUSPARSE_INDEX_BASE_ZERO);
@@ -860,11 +814,7 @@ void SparseMatrix::AddMult(const Vector &x, Vector &y, const real_t a) const
         vecX_descr,
         &beta,
         vecY_descr,
-#ifdef MFEM_USE_SINGLE
-         MFEM_CUDA_or_HIP(_R_32F),
-#else
-         MFEM_CUDA_or_HIP(_R_64F),
-#endif
+         MFEM_CUDA_or_HIP_REAL_T,
         MFEM_GPUSPARSE_ALG,
         &newBufferSize);

@@ -891,11 +841,7 @@ void SparseMatrix::AddMult(const Vector &x, Vector &y, const real_t a) const
         vecX_descr,
         &beta,
         vecY_descr,
-#ifdef MFEM_USE_SINGLE
-         MFEM_CUDA_or_HIP(_R_32F),
-#else
-         MFEM_CUDA_or_HIP(_R_64F),
-#endif
+         MFEM_CUDA_or_HIP_REAL_T,
         MFEM_GPUSPARSE_ALG,
         dBuffer);
 #else
@@ -4372,6 +4318,14 @@ SparseMatrix::~SparseMatrix()
 #ifdef MFEM_USE_CUDA_OR_HIP
   if (Device::Allows(Backend::CUDA_MASK | Backend::HIP_MASK))
   {
+#ifdef MFEM_CUDA_1897_WORKAROUND
+      if (dBuffer)
+      {
+         MFEM_Cu_or_Hip(MemFree)(dBuffer);
+         dBuffer = nullptr;
+         bufferSize = 0;
+      }
+#endif
      if (SparseMatrixCount==1)
      {
         if (handle)
@@ -4379,12 +4333,14 @@ SparseMatrix::~SparseMatrix()
            MFEM_cu_or_hip(sparseDestroy)(handle);
            handle = nullptr;
         }
+#ifndef MFEM_CUDA_1897_WORKAROUND
         if (dBuffer)
         {
            MFEM_Cu_or_Hip(MemFree)(dBuffer);
            dBuffer = nullptr;
            bufferSize = 0;
         }
+#endif
      }
      SparseMatrixCount--;
   }
@@ -98,9 +98,17 @@ protected:
 #ifdef MFEM_USE_CUDA_OR_HIP
   // common for hipSPARSE and cuSPARSE
   static int SparseMatrixCount;
+   mutable bool initBuffers = false;
+
+#if defined(MFEM_USE_CUDA) && CUDA_VERSION >= 12300 && CUDA_VERSION < 12602
+   // Workaround for bug CUSPARSE-1897
+#define MFEM_CUDA_1897_WORKAROUND
+   mutable size_t bufferSize = 0;
+   mutable void *dBuffer = nullptr;
+#else
   static size_t bufferSize;
   static void *dBuffer;
-   mutable bool initBuffers = false;
+#endif

 #if defined(MFEM_USE_CUDA)
   cusparseStatus_t status;
@@ -92,18 +92,6 @@ struct LpReducer
   }
 };

-static Array<real_t>& vector_workspace()
-{
-   static Array<real_t> instance;
-   return instance;
-}
-
-static Array<DevicePair<real_t, real_t>> &Lpvector_workspace()
-{
-   static Array<DevicePair<real_t, real_t>> instance;
-   return instance;
-}
-
 Vector::Vector(const Vector &v)
 {
   const int s = v.Size();
@@ -991,7 +979,7 @@ real_t Vector::Norml2() const
         }
      }
   },
-   L2Reducer{}, UseDevice(), Lpvector_workspace());
+   L2Reducer{}, UseDevice());
   // final answer
   return res.second * sqrt(res.first);
 }
@@ -1006,7 +994,7 @@ real_t Vector::Normlinf() const
   {
      r = fmax(r, fabs(m_data[i]));
   },
-   MaxReducer<real_t> {}, UseDevice(), vector_workspace());
+   MaxReducer<real_t> {}, UseDevice());
   return res;
 }

@@ -1020,7 +1008,7 @@ real_t Vector::Norml1() const
   {
      r += fabs(m_data[i]);
   },
-   SumReducer<real_t> {}, UseDevice(), vector_workspace());
+   SumReducer<real_t> {}, UseDevice());
   return res;
 }

@@ -1063,7 +1051,7 @@ real_t Vector::Normlp(real_t p) const
            }
         }
      },
-      LpReducer{p}, UseDevice(), Lpvector_workspace());
+      LpReducer{p}, UseDevice());
      // final answer
      return res.second * pow(res.first, 1.0 / p);
   } // end if p < infinity()
@@ -1096,7 +1084,7 @@ real_t Vector::operator*(const Vector &v) const
      {
         r += m_data[i] * v_data[i];
      },
-      SumReducer<real_t> {}, use_dev, vector_workspace());
+      SumReducer<real_t> {}, use_dev);
      return res;
   };

@@ -1167,7 +1155,7 @@ real_t Vector::Min() const
      {
         r = fmin(r, m_data[i]);
      },
-      MinReducer<real_t> {}, use_dev, vector_workspace());
+      MinReducer<real_t> {}, use_dev);
      return res;
   };

@@ -1213,7 +1201,7 @@ real_t Vector::Max() const
      {
         r = fmax(r, m_data[i]);
      },
-      MaxReducer<real_t> {}, use_dev, vector_workspace());
+      MaxReducer<real_t> {}, use_dev);
      return res;
   };

@@ -1248,7 +1236,7 @@ real_t Vector::Sum() const
   {
      r += m_data[i];
   },
-   SumReducer<real_t> {}, UseDevice(), vector_workspace());
+   SumReducer<real_t> {}, UseDevice());
   return res;
 }

@@ -461,7 +461,7 @@ void ConductionOperator::Mult(const Vector &u, Vector &du_dt) const

   Kmat.Mult(u, z);
   z.Neg(); // z = -z
-   K->EliminateVDofsInRHS(ess_tdof_list, u, z);
+   K->ParallelEliminateTDofsInRHS(ess_tdof_list, u, z);

   M_solver.Mult(z, du_dt);
   du_dt.Print();
@@ -483,7 +483,7 @@ void ConductionOperator::ImplicitSolve(const real_t dt,
   MFEM_VERIFY(dt == current_dt, ""); // SDIRK methods use the same dt
   Kmat.Mult(u, z);
   z.Neg();
-   K->EliminateVDofsInRHS(ess_tdof_list, u, z);
+   K->ParallelEliminateTDofsInRHS(ess_tdof_list, u, z);

   T_solver.Mult(z, du_dt);
   du_dt.SetSubVector(ess_tdof_list, 0.0);
@@ -87,6 +87,8 @@
 //
 //   Problem 4: level set: Union of doughnut and swiss cheese shapes
 //     mpirun -np 4 distance -m ../../data/inline-hex.mesh -rs 3 -o 2 -t 1.0 -p 4
+//   Problem 5: point source in mfem mesh.
+//     mpirun -np 4 distance -m ../../data/mfem.mesh -p 5 -rs 3 -t 300.0

 #include <fstream>
 #include <iostream>
@@ -233,7 +235,8 @@ int main(int argc, char *argv[])
                  "1: Circle / sphere level set in 2D / 3D\n\t"
                  "2: 2D sine-looking level set\n\t"
                  "3: Gyroid level set in 2D or 3D\n\t"
-                  "4: Combo of a doughnut and swiss cheese shapes in 3D.");
+                  "4: Combo of a doughnut and swiss cheese shapes in 3D.\n\t"
+                  "5: Point source in MFEM mesh.");
   args.AddOption(&rs_levels, "-rs", "--refine-serial",
                  "Number of times to refine the mesh uniformly in serial.");
   args.AddOption(&order, "-o", "--order",
@@ -299,6 +302,11 @@ int main(int argc, char *argv[])
      ls_coeff = new FunctionCoefficient(doughnut_cheese);
      smooth_steps = 0;
   }
+   else if (problem == 5)
+   {
+      ls_coeff = new DeltaCoefficient(0.0, 0.0, 1000.0);
+      smooth_steps = 0;
+   }
   else { MFEM_ABORT("Unrecognized -problem option."); }

   const real_t dx = AvgElementSize(pmesh);
@@ -306,7 +314,7 @@ int main(int argc, char *argv[])
   if (solver_type == 0)
   {
      auto ds = new HeatDistanceSolver(t_param * dx * dx);
-      if (problem == 0)
+      if (problem == 0 || problem == 5)
      {
         ds->transform = false;
      }
@@ -334,7 +342,7 @@ int main(int argc, char *argv[])
   // Smooth-out Gibbs oscillations from the input level set. The smoothing
   // parameter here is specified to be mesh dependent with length scale dx.
   ParGridFunction filt_gf(&pfes_s);
-   if (problem != 0)
+   if (problem != 0 && problem != 5)
   {
      real_t filter_weight = dx;
      // The normalization-based solver needs a more diffused input.
@@ -22,7 +22,6 @@ using namespace mfem;

 TEST_CASE("Reduce Sum", "[Reduction],[GPU]")
 {
-   Array<int> workspace;
   Array<int> a(1000);
   a.HostReadWrite();
   for (int i = 0; i < a.Size(); ++i)
@@ -36,7 +35,7 @@ TEST_CASE("Reduce Sum", "[Reduction],[GPU]")
      int res = 0;
      mfem::reduce(
      a.Size(), res, [=] MFEM_HOST_DEVICE(int i, int &r) { r += dptr[i]; },
-      SumReducer<int> {}, use_dev, workspace);
+      SumReducer<int> {}, use_dev);
      // correct for even-length summations
      int expected = (AsConst(a)[0] + AsConst(a)[a.Size() - 1]) * a.Size() / 2;
      CAPTURE(use_dev);
@@ -46,7 +45,6 @@ TEST_CASE("Reduce Sum", "[Reduction],[GPU]")

 TEST_CASE("Reduce Mult", "[Reduction],[GPU]")
 {
-   Array<long long> workspace;
   Array<long long> a(64);
   a.HostReadWrite();
   for (int i = 0; i < a.Size(); ++i)
@@ -64,7 +62,7 @@ TEST_CASE("Reduce Mult", "[Reduction],[GPU]")
         mfem::reduce(
            a.Size(), res,
         [=] MFEM_HOST_DEVICE(int i, long long &r) { r *= dptr[i]; },
-         MultReducer<long long> {}, use_dev, workspace);
+         MultReducer<long long> {}, use_dev);
         long long expected = 0;
         CAPTURE(use_dev);
         REQUIRE(res == expected);
@@ -76,7 +74,7 @@ TEST_CASE("Reduce Mult", "[Reduction],[GPU]")
         mfem::reduce(
            a.Size(), res,
         [=] MFEM_HOST_DEVICE(int i, long long &r) { r *= dptr[i]; },
-         MultReducer<long long> {}, use_dev, workspace);
+         MultReducer<long long> {}, use_dev);
         long long expected = 21936950640377856;
         CAPTURE(use_dev);
         REQUIRE(res == expected);
@@ -86,7 +84,6 @@ TEST_CASE("Reduce Mult", "[Reduction],[GPU]")

 TEST_CASE("Reduce BAnd", "[Reduction],[GPU]")
 {
-   Array<unsigned> workspace;
   Array<unsigned> a(10);
   SECTION("{ Bit unset }")
   {
@@ -108,7 +105,7 @@ TEST_CASE("Reduce BAnd", "[Reduction],[GPU]")
         mfem::reduce(
            a.Size(), res,
         [=] MFEM_HOST_DEVICE(int i, unsigned &r) { r &= dptr[i]; },
-         BAndReducer<unsigned> {}, use_dev, workspace);
+         BAndReducer<unsigned> {}, use_dev);
         CAPTURE(use_dev);
         REQUIRE(res == ((~1u) & ~(1u << unset_bit)));
         REQUIRE((res & (1u << unset_bit)) == 0);
@@ -132,7 +129,7 @@ TEST_CASE("Reduce BAnd", "[Reduction],[GPU]")
         mfem::reduce(
            a.Size(), res,
         [=] MFEM_HOST_DEVICE(int i, unsigned &r) { r &= dptr[i]; },
-         BAndReducer<unsigned> {}, use_dev, workspace);
+         BAndReducer<unsigned> {}, use_dev);
         CAPTURE(use_dev);
         REQUIRE(res == (1u << set_bit));
      }
@@ -141,7 +138,6 @@ TEST_CASE("Reduce BAnd", "[Reduction],[GPU]")

 TEST_CASE("Reduce BOr", "[Reduction],[GPU]")
 {
-   Array<unsigned> workspace;
   Array<unsigned> a(0x210);
   a.HostReadWrite();
   for (int i = 0; i < a.Size(); ++i)
@@ -157,7 +153,7 @@ TEST_CASE("Reduce BOr", "[Reduction],[GPU]")
      mfem::reduce(
         a.Size(), res,
      [=] MFEM_HOST_DEVICE(int i, unsigned &r) { r |= dptr[i]; },
-      BOrReducer<unsigned> {}, use_dev, workspace);
+      BOrReducer<unsigned> {}, use_dev);
      CAPTURE(use_dev);
      REQUIRE(res == 0x3ffu);
   }
@@ -165,7 +161,6 @@ TEST_CASE("Reduce BOr", "[Reduction],[GPU]")

 TEST_CASE("Reduce Min", "[Reduction],[GPU]")
 {
-   Array<int> workspace;
   Array<int> a(1000);
   auto hptr = a.HostReadWrite();
   for (int i = 0; i < a.Size(); ++i)
@@ -190,7 +185,7 @@ TEST_CASE("Reduce Min", "[Reduction],[GPU]")
            r = dptr[i];
         }
      },
-      MinReducer<int> {}, use_dev, workspace);
+      MinReducer<int> {}, use_dev);
      CAPTURE(use_dev);
      REQUIRE(res == -10);
   }
@@ -198,7 +193,6 @@ TEST_CASE("Reduce Min", "[Reduction],[GPU]")

 TEST_CASE("Reduce Max", "[Reduction],[GPU]")
 {
-   Array<int> workspace;
   Array<int> a(1000);
   auto hptr = a.HostReadWrite();
   for (int i = 0; i < a.Size(); ++i)
@@ -223,7 +217,7 @@ TEST_CASE("Reduce Max", "[Reduction],[GPU]")
            r = dptr[i];
         }
      },
-      MaxReducer<int> {}, use_dev, workspace);
+      MaxReducer<int> {}, use_dev);
      CAPTURE(use_dev);
      REQUIRE(res == 999 - 10);
   }
@@ -231,7 +225,6 @@ TEST_CASE("Reduce Max", "[Reduction],[GPU]")

 TEST_CASE("Reduce MinMax", "[Reduction],[GPU]")
 {
-   Array<DevicePair<int, int>> workspace;
   Array<int> a(1000);
   auto hptr = a.HostReadWrite();
   for (int i = 0; i < a.Size(); ++i)
@@ -262,7 +255,7 @@ TEST_CASE("Reduce MinMax", "[Reduction],[GPU]")
            r.second = dptr[i];
         }
      },
-      MinMaxReducer<int> {}, use_dev, workspace);
+      MinMaxReducer<int> {}, use_dev);
      CAPTURE(use_dev);
      REQUIRE(res.first == -10);
      REQUIRE(res.second == a.Size() - 11);
@@ -271,7 +264,6 @@ TEST_CASE("Reduce MinMax", "[Reduction],[GPU]")

 TEST_CASE("Reduce ArgMin", "[Reduction],[GPU]")
 {
-   Array<DevicePair<double, int>> workspace;
   Array<double> a(1000);
   auto hptr = a.HostReadWrite();
   for (int i = 0; i < a.Size(); ++i)
@@ -297,7 +289,7 @@ TEST_CASE("Reduce ArgMin", "[Reduction],[GPU]")
            r.second = i;
         }
      },
-      ArgMinReducer<double, int> {}, use_dev, workspace);
+      ArgMinReducer<double, int> {}, use_dev);
      CAPTURE(use_dev);
      REQUIRE(res.first == -10);
      REQUIRE(res.second >= 0);
@@ -308,7 +300,6 @@ TEST_CASE("Reduce ArgMin", "[Reduction],[GPU]")

 TEST_CASE("Reduce ArgMax", "[Reduction],[GPU]")
 {
-   Array<DevicePair<double, int>> workspace;
   Array<double> a(1000);

   auto hptr = a.HostReadWrite();
@@ -337,7 +328,7 @@ TEST_CASE("Reduce ArgMax", "[Reduction],[GPU]")
            r.second = i;
         }
      },
-      ArgMaxReducer<double, int> {}, use_dev, workspace);
+      ArgMaxReducer<double, int> {}, use_dev);
      CAPTURE(use_dev);
      REQUIRE(res.first == a.Size() - 11);
      REQUIRE(res.second >= 0);
@@ -348,7 +339,6 @@ TEST_CASE("Reduce ArgMax", "[Reduction],[GPU]")

 TEST_CASE("Reduce ArgMinMax", "[Reduction],[GPU]")
 {
-   Array<MinMaxLocScalar<double, int>> workspace;
   Array<double> a(1000);
   auto hptr = a.HostReadWrite();
   for (int i = 0; i < a.Size(); ++i)
@@ -383,7 +373,7 @@ TEST_CASE("Reduce ArgMinMax", "[Reduction],[GPU]")
            r.max_loc = i;
         }
      },
-      ArgMinMaxReducer<double, int> {}, use_dev, workspace);
+      ArgMinMaxReducer<double, int> {}, use_dev);
      CAPTURE(use_dev);
      REQUIRE(res.min_val == -10);
      REQUIRE(res.min_loc >= 0);
@@ -373,6 +373,7 @@ TEST_CASE("Batched Linear Algebra",
   const int n_rhs = 2;

   DenseTensor A_batch(n, n, n_mat);
+   DenseTensor A_inv_batch(n, n, n_mat);
   Vector x_batch(n * n_rhs * n_mat), y_batch(n * n_rhs * n_mat);
   std::vector<DenseMatrix> As;
   std::vector<DenseMatrix> xs, ys;
@@ -404,6 +405,7 @@ TEST_CASE("Batched Linear Algebra",
      ys.back() = 0.0;
      AddMult_a(1.5, As.back(), xs.back(), ys.back());
      A_batch(i) = As.back();
+      A_inv_batch(i) = As.back();
   }

   // Test batched matrix-vector products
@@ -463,6 +465,33 @@ TEST_CASE("Batched Linear Algebra",
         }
      }
   }
+
+   // Test batched matrix inverse
+   BatchedLinAlg::Get(backend).Invert(A_inv_batch);
+   A_inv_batch.HostReadWrite();
+   Vector output_col(n);
+   Vector col;
+   for (int i = 0; i < n_mat; ++i)
+   {
+      DenseMatrix Ai_inv(A_inv_batch(i));
+      for (int j = 0; j < n; ++j)
+      {
+         output_col = 0.0;
+         As[i].GetColumnReference(j, col);
+         Ai_inv.Mult(col, output_col);
+         for (int k = 0; k < n; ++k)
+         {
+            if (j == k)
+            {
+               REQUIRE(output_col(k) == MFEM_Approx(1.0));
+            }
+            else
+            {
+               REQUIRE(output_col(k) == MFEM_Approx(0.0));
+            }
+         }
+      }
+   }
 }

 TEST_CASE("DenseTensor copy", "[DenseMatrix][DenseTensor]")
@@ -240,4 +240,54 @@ TEST_CASE("SparseMatrix printing", "[SparseMatrix]")
   }
 }

+TEST_CASE("SparseMatrix cuSPARSE Bug", "[SparseMatrix][GPU]")
+{
+   // This test case ensures that we have a functioning workaround for the bug
+   // CUSPARSE-1897. In versions of cuSPARSE before 12.8, the internal buffer
+   // used for cusparseSpMV must be the same when it is called with the same
+   // matrix.
+   //
+   // By default, MFEM uses one buffer, that is shared by all sparse matrices.
+   // In the code below, a buffer is created for A, then modified for B, then
+   // used again for A. Without the workaround, this fails with cuSPARSE version
+   // earlier than 12.8 (confirmed to fail with 12.4).
+
+   const int n = 100;
+   SparseMatrix A(n, n);
+   Vector d(n);
+   d.Randomize(1);
+   for (int i = 0; i < n; ++i)
+   {
+      A.Set(i, i, d[i]);
+   }
+   A.Finalize();
+
+   Vector x(n);
+   x = 1.0;
+
+   Vector y(n);
+   A.Mult(x, y);
+
+   {
+      SparseMatrix B(20, 20);
+      for (int i = 0; i < 20; ++i)
+      {
+         for (int j = 0; j < 20; ++j)
+         {
+            B.Set(i, j, 1.0);
+         }
+      }
+      B.Finalize();
+      Vector u(20);
+      u = 1.0;
+      Vector v(20);
+      B.Mult(u, v);
+   }
+
+   A.Mult(x, y);
+
+   y -= d;
+   REQUIRE(y.Normlinf() == MFEM_Approx(0.0));
+}
+
 } // namespace mfem
Author	SHA1	Message	Date
Will Pazner	552d6857cb	Add GPU scan support for Array<T>::PartialSum	2025-10-21 09:57:37 -07:00
Will Pazner	a37d46e917	Fix comments in scan.hpp	2025-10-21 09:57:37 -07:00
Will Pazner	4acdb072b6	Add more device support to Array<T>	2025-10-21 09:57:37 -07:00
Will Pazner	9dbb184537	Remove need to explicitly pass workspace array to reducers	2025-10-21 09:57:37 -07:00
Will Pazner	d67762a1c9	Move contents of array.cpp to array.hpp Remove explicit template instantiations	2025-10-20 15:38:54 -07:00
Tzanio Kolev	cd4e583f9f	Merge pull request #4659 from mfem/najlkin/parform-impro Improvements of Par(Mixed)BilinearForm and Par(Block)NonlinearForm	2025-10-18 10:48:44 -07:00
Tzanio Kolev	ee94776558	Merge pull request #5066 from mfem/najlkin/fix-nc-blknlform [BUG] Non-conforming meshes in BlockNonlinearForm	2025-10-18 10:48:25 -07:00
Tzanio Kolev	f9be12c539	Merge pull request #5060 from mfem/hypre-cmake Remove hypre+umpire check in CMake	2025-10-15 16:13:01 -07:00
Tzanio Kolev	af478afd00	Merge branch 'master' into najlkin/fix-nc-blknlform	2025-10-15 16:11:05 -07:00
Tzanio Kolev	519e86f134	Merge pull request #5069 from mfem/stale-action-update Fix duplicated job name in Stale action	2025-10-14 15:47:35 -07:00
Mittal, Ketan	f2578fd479	change action name	2025-10-14 09:01:53 -07:00
Tzanio Kolev	563cd25971	Merge pull request #5016 from mfem/stale-action-update Stale action for PRs in-review	2025-10-12 17:03:19 -07:00
Tzanio Kolev	875c62ca1f	Merge pull request #5047 from mfem/fix-magma-inverse fix magma's Invert	2025-10-12 17:02:42 -07:00
Tzanio Kolev	731b0b7049	Merge pull request #5053 from mfem/mfem-mesh MFEM mesh	2025-10-12 17:00:40 -07:00
Tzanio Kolev	3c77ad25e9	Merge pull request #4524 from adam-sim-dev/cusparseXcsrsort Use cusparseXcsrsort instead of the deprecated cusparseXcsru2csr	2025-10-12 16:58:01 -07:00
Jan Nikl	e13d1a1d53	Fixed support of non-conforming meshes in BlockNonlinearForm.	2025-10-10 13:32:53 -07:00
Victor A. P. Magri	09e73b5cee	Remove check	2025-10-09 17:11:13 -04:00
Victor A. P. Magri	464b441409	Remove check	2025-10-09 17:10:09 -04:00
Andrew Ho	8ae9b217f9	Use CUDA_DIR instead of CUDA_HOME	2025-10-09 13:51:57 -07:00
Victor A. P. Magri	3f0a5ad501	Add check for HYPRE_USING_UMPIRE	2025-10-09 12:49:42 -04:00
Ketan Mittal	dfb697b9ae	Merge branch 'master' into stale-action-update	2025-10-07 17:15:22 -07:00
Will PaznerandAndrew Ho	cf01c6cc8b	Add -isystem to NVCC_FLAGS This will silence compiler warnings that are internal to CUDA headers Co-authored-by: Andrew Ho <ho37@llnl.gov>	2025-10-07 16:21:27 -07:00
Will Pazner	e4dc9a3941	Small change to indentation style	2025-10-07 16:20:06 -07:00
Andrew Ho	6a104e94e1	Merge branch 'master' into cusparseXcsrsort	2025-10-07 15:00:06 -07:00
Jan Nikl	b7a0b2cf9a	Renamed ParallelAssembleInternal(Matrix) and added more documentation.	2025-10-07 12:28:31 -07:00
nbeams	27a60a4ca3	Change name of output variable in batched inverse unit test	2025-10-06 22:20:10 +00:00
nbeams	b1d5104fed	Add unit test for batched matrix inverse	2025-10-06 21:50:54 +00:00
Tzanio Kolev	26eecf7575	Merge branch 'master' into cusparseXcsrsort	2025-10-05 13:32:37 -07:00
Tzanio Kolev	bce6e2ca76	Merge branch 'master' into najlkin/parform-impro	2025-10-05 13:25:11 -07:00
Tzanio Kolev	26f8f0e955	Merge branch 'master' into fix-magma-inverse	2025-10-05 13:20:43 -07:00
Mittal, Ketan	30249a4857	minor	2025-10-01 14:37:32 -07:00
Mittal, Ketan	8a9bdf863a	minor	2025-10-01 14:34:29 -07:00
Mittal, Ketan	b72368c698	add mesh and modify distance miniapp to use the new mesh	2025-10-01 14:30:22 -07:00
adam-sim-dev	e6224e00aa	Merge branch 'master' into cusparseXcsrsort	2025-09-30 10:09:49 +08:00
Tom Stitt	aa5b2d4b95	wrong arg	2025-09-29 09:04:33 -07:00
adam-sim-dev	1eb3d185d5	Merge branch 'master' into cusparseXcsrsort	2025-09-25 09:06:24 +08:00
adam-sim-dev	7bc231eb85	Merge branch 'master' into cusparseXcsrsort	2025-09-22 10:31:53 +08:00
Mittal, Ketan	fd01ef188f	Merge branch 'stale-action-update' of https://github.com/mfem/mfem into stale-action-update	2025-09-17 10:14:48 -07:00
Mittal, Ketan	f2a29326c0	Merge branch 'master' of https://github.com/mfem/mfem into stale-action-update	2025-09-17 10:14:35 -07:00
Mittal, Ketan	1f4a40d9a0	stale action for wip PRs	2025-09-17 10:14:26 -07:00
Tzanio Kolev	e671fa5737	Merge branch 'master' into stale-action-update	2025-09-17 03:32:29 -07:00
Tzanio Kolev	40d1550fd6	Merge branch 'master' into najlkin/parform-impro	2025-09-17 03:31:41 -07:00
Tzanio Kolev	ae33f831eb	Merge branch 'master' into cusparseXcsrsort	2025-09-17 03:28:05 -07:00
Jan Nikl	193f8a6801	Partially reverted const modifiers in HypreParMatrix::Copy(Bool)CSR().	2025-09-16 16:19:03 -07:00
Jan Nikl	110720dd04	Fixed documentation of Par(Mixed)BilinearForm::KeepNbrBlock().	2025-09-16 15:39:35 -07:00
Jan Nikl	3fd335c77b	Fixed usage of EliminateVDofsInRHS().	2025-09-16 15:36:04 -07:00
Ketan Mittal	90e1f08bc7	Merge branch 'master' into stale-action-update	2025-09-15 17:14:09 -07:00
Mittal, Ketan	37d004206a	minor	2025-09-12 11:09:22 -07:00
Mittal, Ketan	e1d9b7cfff	stale in-review prs	2025-09-12 11:05:00 -07:00
Will Pazner	8425377bd1	Rename MFEM_REAL_T to MFEM_CUDA_or_HIP_REAL_T	2025-09-09 16:47:35 -07:00
adam-sim-dev	992c5dfc4c	Merge branch 'master' into cusparseXcsrsort	2025-09-09 09:28:10 +08:00
adam-sim-dev	ec94332c2a	Merge branch 'master' into cusparseXcsrsort	2025-09-01 07:35:11 +08:00
Will Pazner	e0c1018c9c	Adjust cuda versions for CUSPARSE-1897 workaround	2025-08-07 16:32:54 -07:00
Will Pazner	7f263221dd	Workaround and unit test for bug CUSPARSE-1897	2025-08-07 15:21:24 -07:00
adam-sim-dev	f82e8c04f9	Merge branch 'master' into cusparseXcsrsort	2025-08-04 11:31:33 +08:00
adam-sim-dev	c731dee195	Merge branch 'master' into cusparseXcsrsort	2025-07-28 10:23:43 +08:00
adam-sim-dev	64c78ff17e	Merge branch 'master' into cusparseXcsrsort	2025-07-23 10:01:12 +08:00
adam-sim-dev	07c0f31e37	Merge branch 'master' into cusparseXcsrsort	2025-07-15 09:21:19 +08:00
adam-sim-dev	5cf82dc903	Merge branch 'master' into cusparseXcsrsort	2025-06-27 09:30:04 +08:00
adam-sim-dev	f19ec00b0a	Merge branch 'master' into cusparseXcsrsort	2025-06-24 12:03:46 +08:00
Andrew Ho	9480e5c5bb	Merge branch 'master' into cusparseXcsrsort	2025-06-16 12:21:18 -07:00
Jan Nikl	058fdaae3f	Implemented gradient of ParBlockNonlinearForm with shared face contributions.	2025-06-11 16:23:42 -07:00
Jan Nikl	9a92e4875b	Implemented Mult of ParBlockNonlinearForm with shared face contributions.	2025-06-11 16:22:16 -07:00
adam-sim-dev	cc21011998	Merge pull request #3 from mfem/cusparseXcsrsort-update cusparseXcsrsort-update by pazner	2025-06-05 14:04:25 +08:00
Will Pazner	09aa9374a9	Use MFEM_REAL_T macro for CUDA/HIP datatype	2025-06-04 20:34:10 -07:00
Will Pazner	78e8cdd7e8	Update SparseMatrix::SortColumnIndices (CUDA or HIP) Use non-deprecated APIs, and unify the CUDA and HIP implementations.	2025-06-04 14:09:20 -07:00
Andrew Ho	04a3fd9bb2	Merge branch 'master' into cusparseXcsrsort	2025-06-02 11:38:00 -07:00
Jan Nikl	9dab032bd0	Merge branch 'master' into najlkin/parform-impro	2025-04-24 15:40:44 -07:00
Jan Nikl	cdfe8102ae	Added gradient of ParNonlinearForm with face integrators.	2025-03-18 16:44:51 -07:00
Jan Nikl	5b82bf0328	Revert "WIP: Added support for trace face integrators in ParMixedBilinearForm." This reverts commit `323ee572b6`.	2025-03-06 06:00:07 -08:00
adam-sim-dev	90f33b1a86	Merge branch 'master' into cusparseXcsrsort	2025-02-19 08:36:21 +08:00
adam-sim-dev	136ae2d98f	Merge branch 'master' into cusparseXcsrsort	2025-02-13 09:53:40 +08:00
adam-sim-dev	449a57d9ad	Merge branch 'master' into cusparseXcsrsort	2025-02-05 09:41:21 +08:00
Tzanio Kolev	41b65d6333	Merge branch 'master' into najlkin/parform-impro	2025-02-04 14:55:19 -08:00
adam-sim-dev	7ac5b61955	Merge branch 'master' into cusparseXcsrsort	2025-01-31 18:28:18 +08:00
adam-sim-dev	c98f117689	Merge branch 'master' into cusparseXcsrsort	2025-01-27 21:48:43 +08:00
Tzanio Kolev	337a0298bf	Merge branch 'master' into cusparseXcsrsort	2025-01-19 13:08:07 -08:00
adam-sim-dev	7930b9b3ca	Merge branch 'master' into cusparseXcsrsort	2025-01-13 12:09:24 +08:00
Jan Nikl	2534d2207d	Fixed name of ParallelEliminateTrialEssentialBC().	2025-01-09 10:49:10 -08:00
Jan Nikl	260b817b3c	Merge branch 'master' into najlkin/parform-impro	2025-01-09 10:01:52 -08:00
Jan Nikl	613d5dd826	Fixed constness in some HyperParMatrix constructors.	2025-01-08 17:56:42 -08:00
Jan Nikl	323ee572b6	WIP: Added support for trace face integrators in ParMixedBilinearForm.	2025-01-08 17:56:05 -08:00
Jan Nikl	56ff5ac5bb	Added support for interior face integrators to ParMixedBilinearForm.	2025-01-08 17:55:42 -08:00
Jan Nikl	e1a06bd6c8	Added methods to Par(Mixed)BilinearForm for elimination of essential BCs.	2025-01-08 17:55:08 -08:00
Jan Nikl	9acae54669	Extended ParMixedBilinearForm methods for parallel assembly.	2025-01-08 17:54:28 -08:00
Jan Nikl	7d92e22a45	Added ParallelAssembleInternal() method to ParBilinearForm.	2025-01-08 17:53:21 -08:00
adam-sim-dev	c0d8ba243d	Merge branch 'master' into cusparseXcsrsort	2024-12-30 08:30:06 +08:00
adam-sim-dev	c946067b9b	Merge branch 'master' into cusparseXcsrsort	2024-12-19 04:49:25 +08:00
adam-sim-dev	c35e5ad7fe	Merge branch 'master' into cusparseXcsrsort	2024-12-08 18:57:22 +08:00
adam-sim-dev	059f0acee6	Merge branch 'master' into cusparseXcsrsort	2024-12-06 02:39:06 +08:00
adam-sim-dev	e18d795334	Merge branch 'master' into cusparseXcsrsort	2024-11-27 09:13:48 +08:00
adam-sim-dev	243b234033	Merge branch 'master' into cusparseXcsrsort	2024-11-21 11:22:17 +08:00
adam-sim-dev	5afe2a09a3	Merge branch 'master' into cusparseXcsrsort	2024-11-18 11:06:52 +08:00
adam-sim-dev	64f8ab42c1	Merge branch 'master' into cusparseXcsrsort	2024-11-13 10:55:58 +08:00
adam-sim-dev	5fb9fc8ec5	Merge branch 'master' into cusparseXcsrsort	2024-11-13 08:41:35 +08:00
adam-sim-dev	587d437f32	Merge branch 'master' into cusparseXcsrsort	2024-11-05 08:43:47 +08:00
adam-sim-dev	c0ae0f0a4b	Merge branch 'master' into cusparseXcsrsort	2024-11-04 08:28:10 +08:00
adam-sim-dev	179e6ec141	Merge branch 'master' into cusparseXcsrsort	2024-10-29 09:42:22 +08:00
adam-sim-dev	2434d4ac71	Merge branch 'master' into cusparseXcsrsort	2024-10-24 16:48:24 +08:00
adam-sim-devandAndrew Ho	a9e5c58897	Correct descrA to matA_descr Co-authored-by: Andrew Ho <ho37@llnl.gov>	2024-10-16 08:21:37 +08:00
adam-sim-dev	35d9fa1f6c	Merge branch 'master' into cusparseXcsrsort	2024-10-15 12:15:04 +08:00
adam-sim-dev	fb9117e9fb	Merge branch 'master' into cusparseXcsrsort	2024-10-10 22:51:24 +08:00
adam-sim-dev	d4d149a5ff	Use cusparseXcsrsort instead of the deprecated cusparseXcsru2csr	2024-09-29 14:27:26 +08:00