Merge branch 'master' into integrate-nontensor-qf

2025-10-08 11:03:36 -07:00 · 2025-10-06 11:44:39 -07:00 · 2025-09-17 11:34:09 -07:00 · 2025-09-05 17:00:24 -07:00
28 changed files with 582 additions and 1983 deletions
@@ -29,47 +29,3 @@ 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,6 +197,16 @@ 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) -isystem "$(CUDA_DIR)/include"
+ -arch=$(CUDA_ARCH)
 # Prefixes for passing flags to the host compiler and linker when using
 # CUDA_CXX=nvcc
 CUDA_XCOMPILER = -Xcompiler=
@@ -1,593 +0,0 @@
-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, bool finalize) const
+Operator &NonlinearForm::GetGradient(const Vector &x) const
 {
   if (ext)
   {
@@ -644,8 +644,6 @@ Operator &NonlinearForm::GetGradient(const Vector &x, bool finalize) const
      }
   }

-   if (!finalize) { return *Grad; }
-
   if (!Grad->Finalized())
   {
      Grad->Finalize(skip_zeros);
@@ -1205,14 +1203,7 @@ const BlockVector &BlockNonlinearForm::Prolongate(const BlockVector &bx) const
      aux1.Update(block_offsets);
      for (int s = 0; s < fes.Size(); s++)
      {
-         if (P[s])
-         {
-            P[s]->Mult(bx.GetBlock(s), aux1.GetBlock(s));
-         }
-         else
-         {
-            aux1.GetBlock(s) = bx.GetBlock(s);
-         }
+         P[s]->Mult(bx.GetBlock(s), aux1.GetBlock(s));
      }
      return aux1;
   }
@@ -1241,16 +1232,11 @@ 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,
-                                                bool finalize) const
+void BlockNonlinearForm::ComputeGradientBlocked(const BlockVector &bx) const
 {
   const int skip_zeros = 0;
   Array<Array<int> *> vdofs(fes.Size());
@@ -1504,7 +1490,7 @@ void BlockNonlinearForm::ComputeGradientBlocked(const BlockVector &bx,
      }
   }

-   if (finalize && !Grads(0,0)->Finalized())
+   if (!Grads(0,0)->Finalized())
   {
      for (int i=0; i<fes.Size(); ++i)
      {
@@ -1543,23 +1529,7 @@ Operator &BlockNonlinearForm::GetGradient(const Vector &x) const
         for (int s2 = 0; s2 < fes.Size(); ++s2)
         {
            delete cGrads(s1, 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;
-            }
+            cGrads(s1, s2) = RAP(*cP[s1], *Grads(s1, s2), *cP[s2]);
            mGrads(s1, s2) = cGrads(s1, s2);
         }
      }
@@ -217,12 +217,7 @@ 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 { 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;
+   Operator &GetGradient(const Vector &x) const override;

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

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

 public:
   /// Construct an empty BlockNonlinearForm. Initialize with SetSpaces().
@@ -151,15 +151,6 @@ 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();
@@ -342,15 +333,6 @@ 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
@@ -362,26 +344,6 @@ 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
 {
@@ -523,7 +485,7 @@ void ParBilinearForm::FormLinearSystem(
      HypreParVector true_X(pfes), true_B(pfes);
      P.MultTranspose(b, true_B);
      R.Mult(x, true_X);
-      ParallelEliminateTDofsInRHS(ess_tdof_list, true_X, true_B);
+      p_mat.EliminateBC(p_mat_e, ess_tdof_list, true_X, true_B);
      R.MultTranspose(true_B, b);
      hybridization->ReduceRHS(true_B, B);
      X.SetSize(B.Size());
@@ -536,11 +498,17 @@ void ParBilinearForm::FormLinearSystem(
      B.SetSize(X.Size());
      P.MultTranspose(b, B);
      R.Mult(x, X);
-      ParallelEliminateTDofsInRHS(ess_tdof_list, X, B);
+      p_mat.EliminateBC(p_mat_e, 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)
 {
@@ -585,7 +553,7 @@ void ParBilinearForm::FormSystemMatrix(const Array<int> &ess_tdof_list,
         mat = NULL;
         delete mat_e;
         mat_e = NULL;
-         ParallelEliminateTDofs(ess_tdof_list);
+         p_mat_e.EliminateRowsCols(p_mat, ess_tdof_list);
      }
      if (hybridization)
      {
@@ -647,180 +615,36 @@ void ParBilinearForm::Update(FiniteElementSpace *nfes)
   p_mat_e.Clear();
 }

-void ParMixedBilinearForm::pAllocMat()
-{
-   const int trial_nbr_size = trial_pfes->GetFaceNbrVSize();
-   const int test_nbr_size = test_pfes->GetFaceNbrVSize();

-   if (keep_nbr_block)
-   {
-      mat = new SparseMatrix(height + test_nbr_size, width + trial_nbr_size);
-   }
-   else
-   {
-      mat = new SparseMatrix(height, width + trial_nbr_size);
-   }
+HypreParMatrix *ParMixedBilinearForm::ParallelAssemble()
+{
+   // 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);
+
+   HypreParMatrix *rap = RAP(test_pfes->Dof_TrueDof_Matrix(), A,
+                             trial_pfes->Dof_TrueDof_Matrix());
+
+   delete A;
+
+   return rap;
 }

-void ParMixedBilinearForm::AssembleSharedFaces(int skip_zeros)
+void ParMixedBilinearForm::ParallelAssemble(OperatorHandle &A)
 {
-   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);
-   }
+   // 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);

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

@@ -846,44 +670,6 @@ 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,
@@ -904,8 +690,10 @@ void ParMixedBilinearForm::FormRectangularSystemMatrix(
      mat = NULL;
      delete mat_e;
      mat_e = NULL;
-      ParallelEliminateTrialTDofs(trial_tdof_list);
-      ParallelEliminateTestTDofs(test_tdof_list);
+      HypreParMatrix *temp =
+         p_mat.As<HypreParMatrix>()->EliminateCols(trial_tdof_list);
+      p_mat.As<HypreParMatrix>()->EliminateRows(test_tdof_list);
+      p_mat_e.Reset(temp, true);
   }

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

-   ParallelEliminateTrialTDofsInRHS(trial_tdof_list, X, B);
+   p_mat_e.As<HypreParMatrix>()->Mult(-1.0, X, 1.0, 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,14 +101,6 @@ 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); }
@@ -154,13 +146,6 @@ 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
@@ -172,12 +157,6 @@ 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
@@ -190,28 +169,6 @@ 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;
@@ -281,6 +238,8 @@ public:

   void Update(FiniteElementSpace *nfes = NULL) override;

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

@@ -298,13 +257,6 @@ 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 &);
@@ -324,7 +276,6 @@ public:
   {
      trial_pfes = trial_fes;
      test_pfes  = test_fes;
-      keep_nbr_block = false;
   }

   /** @brief Create a ParMixedBilinearForm on the given FiniteElementSpace%s
@@ -344,89 +295,15 @@ 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.
-   /** 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);
+   HypreParMatrix *ParallelAssemble();

   /** @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) { 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);
+   void ParallelAssemble(OperatorHandle &A);

   using MixedBilinearForm::FormRectangularSystemMatrix;
   using MixedBilinearForm::FormRectangularLinearSystem;
@@ -105,59 +105,6 @@ 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); }
@@ -165,61 +112,19 @@ Operator &ParNonlinearForm::GetGradient(const Vector &x) const
   ParFiniteElementSpace *pfes = ParFESpace();

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

-   if (fnfi.Size())
+   NonlinearForm::GetGradient(x); // (re)assemble Grad, no b.c.
+
+   OperatorHandle dA(pGrad.Type()), Ph(pGrad.Type());
+
+   if (fnfi.Size() == 0)
   {
-      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);
+      dA.MakeSquareBlockDiag(pfes->GetComm(), pfes->GlobalVSize(),
+                             pfes->GetDofOffsets(), Grad);
   }
   else
   {
-      NonlinearForm::GetGradient(x); // (re)assemble Grad, no b.c.
-
-      dA.MakeSquareBlockDiag(pfes->GetComm(), pfes->GlobalVSize(),
-                             pfes->GetDofOffsets(), Grad);
+      MFEM_ABORT("TODO: assemble contributions from shared face terms");
   }

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

   if (fnfi.Size() > 0)
   {
-      // 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];
-      }
+      MFEM_ABORT("TODO: assemble contributions from shared face terms");
   }

   for (int s=0; s<fes.Size(); ++s)
@@ -486,106 +328,6 @@ 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)
@@ -605,155 +347,49 @@ BlockOperator & ParBlockNonlinearForm::GetGradient(const Vector &x) const
      }
   }

-   // 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));
-   }
+   GetLocalGradient(x); // gradients are stored in 'Grads'

   if (fnfi.Size() > 0)
   {
-      const int skip_zeros = 0;
-
-      for (int s=0; s<fes.Size(); ++s)
-      {
-         const_cast<ParFiniteElementSpace*>(pfes[s])->ExchangeFaceNbrData();
-      }
-
-      for (int s1=0; s1<fes.Size(); ++s1)
-      {
-         for (int s2=0; s2<fes.Size(); ++s2)
-         {
-            if (Grads(s1,s2) == NULL)
-            {
-               int nbr_size = pfes[s2]->GetFaceNbrVSize();
-               Grads(s1,s2) = new SparseMatrix(pfes[s1]->GetVSize(),
-                                               pfes[s2]->GetVSize() + nbr_size);
-            }
-         }
-      }
-
-      // (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)
-         {
-            Grads(s1,s2)->Finalize(skip_zeros);
-         }
-
-      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());
-         }
-      }
+      MFEM_ABORT("TODO: assemble contributions from shared face terms");
   }
-   else
+
+   for (int s1=0; s1<fes.Size(); ++s1)
   {
-      // (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)
      {
-         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)
         {
-            OperatorHandle dA(phBlockGrad(s1,s2)->Type()),
-                           Ph(phBlockGrad(s1,s2)->Type()),
-                           Rh(phBlockGrad(s1,s2)->Type());
+            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);

-            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());
+            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,8 +29,6 @@ protected:
   mutable ParGridFunction X, Y;
   mutable OperatorHandle pGrad;

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

@@ -83,8 +81,6 @@ 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;
@@ -315,13 +315,36 @@ void FaceQuadratureSpace::Save(std::ostream &os) const

 const Vector &FaceQuadratureSpace::GetGeometricFactorWeights() const
 {
-   auto flags = FaceGeometricFactors::DETERMINANTS;
-   // TODO: assumes only one integration rule. This should be fixed once
-   // Mesh::GetFaceGeometricFactors acceps a QuadratureSpace instead of
-   // IntegrationRule.
-   const IntegrationRule &ir = GetIntRule(0);
-   auto *geom = mesh.GetFaceGeometricFactors(ir, flags, face_type);
-   return geom->detJ;
+   if (mesh.MeshGenerator() == (1 << 1)) // only tensor product elements
+   {
+      auto flags = FaceGeometricFactors::DETERMINANTS;
+      const IntegrationRule &ir = GetIntRule(0);
+      auto *geom = mesh.GetFaceGeometricFactors(ir, flags, face_type);
+      return geom->detJ;
+   }
+   else
+   {
+      if (non_tensor_weights.Size() == 0)
+      {
+         non_tensor_weights.SetSize(size);
+         int idx = 0;
+         for (int f = 0; f < face_indices.Size(); ++f)
+         {
+            const int f_idx = face_indices[f];
+            const Geometry::Type geom = mesh.GetFaceGeometry(f_idx);
+            ElementTransformation &T = *mesh.GetFaceTransformation(f_idx);
+            const IntegrationRule &ir = *int_rule[geom];
+            for (int q = 0; q < ir.Size(); ++q)
+            {
+               const IntegrationPoint &ip = ir.IntPoint(q);
+               T.SetIntPoint(&ip);
+               non_tensor_weights[idx] = T.Weight();
+               idx += 1;
+            }
+         }
+      }
+      return non_tensor_weights;
+   }
 }

 } // namespace mfem
@@ -220,6 +220,9 @@ class FaceQuadratureSpace : public QuadratureSpaceBase
   /// Inverse of the map @a face_indices.
   const std::unordered_map<int,int> &face_indices_inv;

+   /// Integration weights (used only in the non-tensor case).
+   mutable Vector non_tensor_weights;
+
   const Vector &GetGeometricFactorWeights() const override;
   void ConstructOffsets();
   void Construct();
@@ -0,0 +1,214 @@
+// 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,13 +16,9 @@
 #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>
@@ -139,8 +135,6 @@ 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(); }

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

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

   /** @brief Find the minimal element in the array, using the comparison
       operator `<` for class T. */
-   inline T Min() const;
+   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); }
@@ -303,22 +297,22 @@ public:
   }

   /// Return 1 if the array is sorted from lowest to highest.  Otherwise return 0.
-   inline int IsSorted() const;
+   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.
-   inline bool IsConstant() const;
+   bool IsConstant() const;

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

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

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

   /// Set all entries of the array to the provided constant.
   inline void operator=(const T &a);
@@ -803,14 +797,8 @@ template <typename T> template <typename CT>
 inline Array<T> &Array<T>::operator=(const Array<CT> &src)
 {
   SetSize(src.Size());
-
-   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];
-   });
+   for (int i = 0; i < size; i++) { (*this)[i] = T(src[i]); }
+   return *this;
 }

 template <class T>
@@ -1026,24 +1014,19 @@ template <class T>
 inline void Array<T>::GetSubArray(int offset, int sa_size, Array<T> &sa) const
 {
   sa.SetSize(sa_size);
-   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)
+   for (int i = 0; i < sa_size; i++)
   {
-      y[i] = x[offset + i];
-   });
+      sa[i] = (*this)[offset+i];
+   }
 }

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

 template <class T>
@@ -1052,153 +1035,6 @@ 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
@@ -1238,40 +1074,6 @@ 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,6 +12,7 @@
 #ifndef MFEM_REDUCERS_HPP
 #define MFEM_REDUCERS_HPP

+#include "array.hpp"
 #include "forall.hpp"

 #include <cmath>
@@ -513,33 +514,6 @@ 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;
-   }
-};
-
 }

 /**
@@ -555,7 +529,8 @@ public:
 @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)
+void reduce(int N, T &res, B &&body, const R &reducer, bool use_dev,
+            Array<T> &workspace)
 {
   if (N == 0)
   {
@@ -592,7 +567,13 @@ 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 work = internal::ReductionWorkspace<T>::Get(nblocks);
+      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();
      red.work = work;
      forall_2D(nblocks, block_size, 1, std::move(red));
      // wait for results
@@ -28,37 +28,8 @@

 namespace mfem
 {
-
-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<>{})
+/// Equivalent to InclusiveScan(use_dev, d_in, d_out, num_items, workspace,
+/// std::plus<>{})
 template <class InputIt, class OutputIt>
 void InclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items)
 {
@@ -66,12 +37,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))
   {
-      using internal::ScanWorkspace;
-      size_t bytes = ScanWorkspace::Size();
-      if (bytes > 0)
+      static Array<std::byte> workspace;
+      size_t bytes = workspace.Size();
+      if (bytes)
      {
         auto err = MFEM_CUB_NAMESPACE::DeviceScan::InclusiveSum(
-                       ScanWorkspace::Get(bytes), bytes, d_in, d_out, num_items);
+                       workspace.Write(), bytes, d_in, d_out, num_items);
 #if defined(MFEM_USE_CUDA)
         if (err == cudaSuccess)
         {
@@ -86,12 +57,11 @@ 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));
-      // resize buffer (in ScanWorkspace::Get) and try again
+      workspace.SetSize(bytes);
      MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::InclusiveSum(
-                        ScanWorkspace::Get(bytes), bytes, d_in, d_out, num_items));
+                        workspace.Write(), bytes, d_in, d_out, num_items));
      return;
   }
 #endif
@@ -131,13 +101,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))
   {
-      using internal::ScanWorkspace;
-      size_t bytes = ScanWorkspace::Size();
-      if (bytes > 0)
+      static Array<std::byte> workspace;
+      size_t bytes = workspace.Size();
+      if (bytes)
      {
         auto err = MFEM_CUB_NAMESPACE::DeviceScan::InclusiveScan(
-                       ScanWorkspace::Get(bytes), bytes, d_in, d_out, scan_op,
-                       num_items);
+                       workspace.Write(), bytes, d_in, d_out, scan_op, num_items);
 #if defined(MFEM_USE_CUDA)
         if (err == cudaSuccess)
         {
@@ -154,9 +123,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(
-                        ScanWorkspace::Get(bytes), bytes, d_in, d_out, scan_op,
-                        num_items));
+                        workspace.Write(), bytes, d_in, d_out, scan_op, num_items));
      return;
   }
 #endif
@@ -195,13 +164,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))
   {
-      using internal::ScanWorkspace;
-      size_t bytes = ScanWorkspace::Size();
+      static Array<std::byte> workspace;
+      size_t bytes = workspace.Size();
      if (bytes)
      {
         auto err = MFEM_CUB_NAMESPACE::DeviceScan::ExclusiveScan(
-                       ScanWorkspace::Get(bytes), bytes, d_in, d_out, scan_op,
-                       init_value, num_items);
+                       workspace.Write(), bytes, d_in, d_out, scan_op, init_value,
+                       num_items);
 #if defined(MFEM_USE_CUDA)
         if (err == cudaSuccess)
         {
@@ -218,9 +187,10 @@ 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(
-                        ScanWorkspace::Get(bytes), bytes, d_in, d_out, scan_op,
-                        init_value, num_items));
+                        workspace.Write(), bytes, d_in, d_out, scan_op, init_value,
+                        num_items));
      return;
   }
 #endif
@@ -243,7 +213,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,
-/// std::plus<>{})
+/// workspace, 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_LU_ptrs, n, d_P_ptrs, info_array.Write(), n_mat,
+               n, n, d_A_ptrs, n, d_P_ptrs, info_array.Write(), n_mat,
               Magma::Queue());
   MFEM_VERIFY(status == MAGMA_SUCCESS, "");

@@ -561,8 +561,7 @@ 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(const Memory<SrcT> &src, MemoryClass dst_mc,
-                       Memory<DstT> &dst)
+void CopyConvertMemory(Memory<SrcT> &src, MemoryClass dst_mc, Memory<DstT> &dst)
 {
   auto capacity = src.Capacity();
   dst.New(capacity, GetMemoryType(dst_mc));
@@ -843,8 +842,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, const HYPRE_BigInt *rows,
-                                    const HYPRE_BigInt *cols)
+static bool RowAndColStartsAreEqual(MPI_Comm comm, HYPRE_BigInt *rows,
+                                    HYPRE_BigInt *cols)
 {
   const int part_size = GetPartitioningArraySize(comm);
   bool are_equal = true;
@@ -1132,7 +1131,7 @@ HypreParMatrix::HypreParMatrix(
 HypreParMatrix::HypreParMatrix(MPI_Comm comm,
                               HYPRE_BigInt *row_starts,
                               HYPRE_BigInt *col_starts,
-                               const SparseMatrix *sm_a)
+                               SparseMatrix *sm_a)
 {
   MFEM_ASSERT(sm_a != NULL, "invalid input");
   MFEM_VERIFY(!HYPRE_AssumedPartitionCheck(),
@@ -1146,7 +1145,7 @@ HypreParMatrix::HypreParMatrix(MPI_Comm comm,

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

   // NOTE: this call creates a matrix on host even when device support is
@@ -1308,11 +1307,10 @@ HypreParMatrix::HypreParMatrix(MPI_Comm comm, int id, int np,
 HypreParMatrix::HypreParMatrix(MPI_Comm comm, int nrows,
                               HYPRE_BigInt glob_nrows,
                               HYPRE_BigInt glob_ncols,
-                               const int *I,
-                               const HYPRE_BigInt *J,
-                               const real_t *data,
-                               const HYPRE_BigInt *rows,
-                               const HYPRE_BigInt *cols)
+                               int *I, HYPRE_BigInt *J,
+                               real_t *data,
+                               HYPRE_BigInt *rows,
+                               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,
-                  const SparseMatrix *a); // constructor with 4 arguments, v2
+                  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, const int *I, const HYPRE_BigInt *J,
-                  const real_t *data, const HYPRE_BigInt *rows,
-                  const HYPRE_BigInt *cols); // constructor with 9 arguments
+                  HYPRE_BigInt glob_ncols, int *I, HYPRE_BigInt *J,
+                  real_t *data, HYPRE_BigInt *rows,
+                  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,12 +46,6 @@
 #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
 {

@@ -63,10 +57,8 @@ 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()
@@ -472,67 +464,109 @@ void SparseMatrix::SortColumnIndices()
   }

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

-      // 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);
+      cusparseMatDescr_t matA_descr;
+      cusparseCreateMatDescr( &matA_descr );
+      cusparseSetMatIndexBase( matA_descr, CUSPARSE_INDEX_BASE_ZERO );
+      cusparseSetMatType( matA_descr, CUSPARSE_MATRIX_TYPE_GENERAL );

-      // 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);
+      cusparseCreateCsru2csrInfo( &sortInfoA );

-      // 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; });
+#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

-      // 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);
+      CuMemAlloc( &pBuffer, pBufferSizeInBytes );

-      // 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));
+#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 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.
+      // 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.
      MFEM_STREAM_SYNC;

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

-      MFEM_Cu_or_Hip(MemFree)(d_a_unsorted);
-      MFEM_Cu_or_Hip(MemFree)(pBuffer);
+      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
   }
   else
 #endif // MFEM_USE_CUDA_OR_HIP
@@ -787,15 +821,27 @@ 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),
-            MFEM_CUDA_or_HIP_REAL_T);
+#ifdef MFEM_USE_SINGLE
+            MFEM_CUDA_or_HIP(_R_32F));
+#else
+            MFEM_CUDA_or_HIP(_R_64F));
+#endif

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

@@ -841,7 +891,11 @@ void SparseMatrix::AddMult(const Vector &x, Vector &y, const real_t a) const
         vecX_descr,
         &beta,
         vecY_descr,
-         MFEM_CUDA_or_HIP_REAL_T,
+#ifdef MFEM_USE_SINGLE
+         MFEM_CUDA_or_HIP(_R_32F),
+#else
+         MFEM_CUDA_or_HIP(_R_64F),
+#endif
         MFEM_GPUSPARSE_ALG,
         dBuffer);
 #else
@@ -4318,14 +4372,6 @@ 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)
@@ -4333,14 +4379,12 @@ 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,17 +98,9 @@ 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;
-#endif
+   mutable bool initBuffers = false;

 #if defined(MFEM_USE_CUDA)
   cusparseStatus_t status;
@@ -92,6 +92,18 @@ 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();
@@ -979,7 +991,7 @@ real_t Vector::Norml2() const
         }
      }
   },
-   L2Reducer{}, UseDevice());
+   L2Reducer{}, UseDevice(), Lpvector_workspace());
   // final answer
   return res.second * sqrt(res.first);
 }
@@ -994,7 +1006,7 @@ real_t Vector::Normlinf() const
   {
      r = fmax(r, fabs(m_data[i]));
   },
-   MaxReducer<real_t> {}, UseDevice());
+   MaxReducer<real_t> {}, UseDevice(), vector_workspace());
   return res;
 }

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

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

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

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

@@ -1236,7 +1248,7 @@ real_t Vector::Sum() const
   {
      r += m_data[i];
   },
-   SumReducer<real_t> {}, UseDevice());
+   SumReducer<real_t> {}, UseDevice(), vector_workspace());
   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->ParallelEliminateTDofsInRHS(ess_tdof_list, u, z);
+   K->EliminateVDofsInRHS(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->ParallelEliminateTDofsInRHS(ess_tdof_list, u, z);
+   K->EliminateVDofsInRHS(ess_tdof_list, u, z);

   T_solver.Mult(z, du_dt);
   du_dt.SetSubVector(ess_tdof_list, 0.0);
@@ -87,8 +87,6 @@
 //
 //   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>
@@ -235,8 +233,7 @@ 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.\n\t"
-                  "5: Point source in MFEM mesh.");
+                  "4: Combo of a doughnut and swiss cheese shapes in 3D.");
   args.AddOption(&rs_levels, "-rs", "--refine-serial",
                  "Number of times to refine the mesh uniformly in serial.");
   args.AddOption(&order, "-o", "--order",
@@ -302,11 +299,6 @@ 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);
@@ -314,7 +306,7 @@ int main(int argc, char *argv[])
   if (solver_type == 0)
   {
      auto ds = new HeatDistanceSolver(t_param * dx * dx);
-      if (problem == 0 || problem == 5)
+      if (problem == 0)
      {
         ds->transform = false;
      }
@@ -342,7 +334,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 && problem != 5)
+   if (problem != 0)
   {
      real_t filter_weight = dx;
      // The normalization-based solver needs a more diffused input.
@@ -160,8 +160,10 @@ TEST_CASE("Quadrature Function Integration", "[QuadratureFunction][GPU]")
   auto fname = GENERATE(
                   "../../data/star.mesh",
                   "../../data/star-q3.mesh",
+                   "../../data/beam-tri.mesh",
                   "../../data/fichera.mesh",
-                   "../../data/fichera-q3.mesh"
+                   "../../data/fichera-q3.mesh",
+                   "../../data/beam-tet.mesh"
                );
   const int order = GENERATE(1, 2, 3);

@@ -22,6 +22,7 @@ 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)
@@ -35,7 +36,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);
+      SumReducer<int> {}, use_dev, workspace);
      // correct for even-length summations
      int expected = (AsConst(a)[0] + AsConst(a)[a.Size() - 1]) * a.Size() / 2;
      CAPTURE(use_dev);
@@ -45,6 +46,7 @@ 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)
@@ -62,7 +64,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);
+         MultReducer<long long> {}, use_dev, workspace);
         long long expected = 0;
         CAPTURE(use_dev);
         REQUIRE(res == expected);
@@ -74,7 +76,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);
+         MultReducer<long long> {}, use_dev, workspace);
         long long expected = 21936950640377856;
         CAPTURE(use_dev);
         REQUIRE(res == expected);
@@ -84,6 +86,7 @@ TEST_CASE("Reduce Mult", "[Reduction],[GPU]")

 TEST_CASE("Reduce BAnd", "[Reduction],[GPU]")
 {
+   Array<unsigned> workspace;
   Array<unsigned> a(10);
   SECTION("{ Bit unset }")
   {
@@ -105,7 +108,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);
+         BAndReducer<unsigned> {}, use_dev, workspace);
         CAPTURE(use_dev);
         REQUIRE(res == ((~1u) & ~(1u << unset_bit)));
         REQUIRE((res & (1u << unset_bit)) == 0);
@@ -129,7 +132,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);
+         BAndReducer<unsigned> {}, use_dev, workspace);
         CAPTURE(use_dev);
         REQUIRE(res == (1u << set_bit));
      }
@@ -138,6 +141,7 @@ 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)
@@ -153,7 +157,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);
+      BOrReducer<unsigned> {}, use_dev, workspace);
      CAPTURE(use_dev);
      REQUIRE(res == 0x3ffu);
   }
@@ -161,6 +165,7 @@ 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)
@@ -185,7 +190,7 @@ TEST_CASE("Reduce Min", "[Reduction],[GPU]")
            r = dptr[i];
         }
      },
-      MinReducer<int> {}, use_dev);
+      MinReducer<int> {}, use_dev, workspace);
      CAPTURE(use_dev);
      REQUIRE(res == -10);
   }
@@ -193,6 +198,7 @@ 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)
@@ -217,7 +223,7 @@ TEST_CASE("Reduce Max", "[Reduction],[GPU]")
            r = dptr[i];
         }
      },
-      MaxReducer<int> {}, use_dev);
+      MaxReducer<int> {}, use_dev, workspace);
      CAPTURE(use_dev);
      REQUIRE(res == 999 - 10);
   }
@@ -225,6 +231,7 @@ 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)
@@ -255,7 +262,7 @@ TEST_CASE("Reduce MinMax", "[Reduction],[GPU]")
            r.second = dptr[i];
         }
      },
-      MinMaxReducer<int> {}, use_dev);
+      MinMaxReducer<int> {}, use_dev, workspace);
      CAPTURE(use_dev);
      REQUIRE(res.first == -10);
      REQUIRE(res.second == a.Size() - 11);
@@ -264,6 +271,7 @@ 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)
@@ -289,7 +297,7 @@ TEST_CASE("Reduce ArgMin", "[Reduction],[GPU]")
            r.second = i;
         }
      },
-      ArgMinReducer<double, int> {}, use_dev);
+      ArgMinReducer<double, int> {}, use_dev, workspace);
      CAPTURE(use_dev);
      REQUIRE(res.first == -10);
      REQUIRE(res.second >= 0);
@@ -300,6 +308,7 @@ 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();
@@ -328,7 +337,7 @@ TEST_CASE("Reduce ArgMax", "[Reduction],[GPU]")
            r.second = i;
         }
      },
-      ArgMaxReducer<double, int> {}, use_dev);
+      ArgMaxReducer<double, int> {}, use_dev, workspace);
      CAPTURE(use_dev);
      REQUIRE(res.first == a.Size() - 11);
      REQUIRE(res.second >= 0);
@@ -339,6 +348,7 @@ 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)
@@ -373,7 +383,7 @@ TEST_CASE("Reduce ArgMinMax", "[Reduction],[GPU]")
            r.max_loc = i;
         }
      },
-      ArgMinMaxReducer<double, int> {}, use_dev);
+      ArgMinMaxReducer<double, int> {}, use_dev, workspace);
      CAPTURE(use_dev);
      REQUIRE(res.min_val == -10);
      REQUIRE(res.min_loc >= 0);
@@ -373,7 +373,6 @@ 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;
@@ -405,7 +404,6 @@ 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
@@ -465,33 +463,6 @@ 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,54 +240,4 @@ 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
Andrew Ho	633db43d52	Merge branch 'master' into integrate-nontensor-qf	2025-10-08 11:03:36 -07:00
Andrew Ho	4a16deea26	Merge branch 'master' into integrate-nontensor-qf	2025-10-06 11:44:39 -07:00
Andrew Ho	9a266d8248	Merge branch 'master' into integrate-nontensor-qf	2025-09-17 11:34:09 -07:00
Will Pazner	4f08e58a95	Support non-tensor-product meshes in QuadratureFunction::Integrate on faces	2025-09-05 17:00:24 -07:00