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

Fix comments in scan.hpp
Add more device support to Array<T>
2025-10-21 09:57:37 -07:00 · 2025-10-21 09:57:37 -07:00 · 2025-10-21 09:57:37 -07:00 · 2025-10-21 09:57:37 -07:00 · 2025-10-20 15:38:54 -07:00 · 2025-10-18 10:48:44 -07:00
25 changed files with 1211 additions and 462 deletions
@@ -12,7 +12,7 @@ coverage:
        threshold: 0%
        base: auto
        branches:
-          - main
+          - master
        if_ci_failed: error
        informational: true
        only_pulls: true
@@ -22,7 +22,7 @@ coverage:
        threshold: 1%  # allows variations around the target
        base: auto
        branches:
-          - main
+          - master
        if_ci_failed: error
        only_pulls: true

@@ -40,8 +40,6 @@ Discretization improvements
  provided that neighboring hexahedra are not refined in conflicting directions.
  A new ParMesh method is added to check for such conflicts, before refinement.

- Renamed the default GitHub branch from "master" to "main".
-
 Meshing improvements
 --------------------

@@ -652,7 +652,7 @@ if (MFEM_USE_CUDA)
 endif()

 add_subdirectory(config)
-set(MAIN_HEADERS
+set(MASTER_HEADERS
  ${PROJECT_SOURCE_DIR}/mfem.hpp
  ${PROJECT_SOURCE_DIR}/mfem-performance.hpp)

@@ -689,7 +689,7 @@ set(MFEM_SOURCE_DIR ${CMAKE_CURRENT_SOURCE_DIR})
 set(MFEM_INSTALL_DIR ${CMAKE_INSTALL_PREFIX})

 # Declaring the library
-mfem_add_library(mfem ${SOURCES} ${HEADERS} ${MAIN_HEADERS})
+mfem_add_library(mfem ${SOURCES} ${HEADERS} ${MASTER_HEADERS})
 # message(STATUS "TPL_LIBRARIES = ${TPL_LIBRARIES}")
 target_link_libraries(mfem PUBLIC ${TPL_LIBRARIES} ${TPL_TARGETS})
 if (TPL_TARGETS)
@@ -888,12 +888,12 @@ install(TARGETS ${PROJECT_NAME}
  LIBRARY DESTINATION ${INSTALL_LIB_DIR}
  ARCHIVE DESTINATION ${INSTALL_LIB_DIR})

-# Install the main headers
+# Install the master headers
 foreach(Header mfem.hpp mfem-performance.hpp)
  install(FILES ${PROJECT_BINARY_DIR}/InstallHeaders/${Header}
    DESTINATION ${INSTALL_INCLUDE_DIR})
 endforeach()
-install(FILES ${MAIN_HEADERS} DESTINATION ${INSTALL_INCLUDE_DIR}/mfem)
+install(FILES ${MASTER_HEADERS} DESTINATION ${INSTALL_INCLUDE_DIR}/mfem)

 # Install the headers (except common miniapp which is installed from its subdir)
 install(DIRECTORY ${MFEM_SOURCE_DIRS}
@@ -3,10 +3,10 @@
 </p>

 <p align="center">
-<a href="https://github.com/mfem/mfem/blob/main/LICENSE"><img alt="License" src="https://img.shields.io/badge/License-BSD-brightgreen.svg"></a>
-<a href="https://github.com/mfem/mfem/actions?query=workflow%3Arepo-check+branch%3Amain"><img alt="Repo check" src="https://github.com/mfem/mfem/actions/workflows/repo-check.yml/badge.svg?branch=main"></a>
-<a href="https://github.com/mfem/mfem/actions?query=workflow%3Abuild-analysis+branch%3Amain"><img alt="Build Analysis" src="https://github.com/mfem/mfem/actions/workflows/mfem-analysis.yml/badge.svg?branch=main"></a>
-<a href="https://github.com/mfem/mfem/actions?query=workflow%3Abuilds-and-tests+branch%3Amain"><img alt="Builds and Tests" src="https://github.com/mfem/mfem/actions/workflows/builds-and-tests.yml/badge.svg?branch=main"></a>
+<a href="https://github.com/mfem/mfem/blob/master/LICENSE"><img alt="License" src="https://img.shields.io/badge/License-BSD-brightgreen.svg"></a>
+<a href="https://github.com/mfem/mfem/actions?query=workflow%3Arepo-check+branch%3Amaster"><img alt="Repo check" src="https://github.com/mfem/mfem/actions/workflows/repo-check.yml/badge.svg?branch=master"></a>
+<a href="https://github.com/mfem/mfem/actions?query=workflow%3Abuild-analysis+branch%3Amaster"><img alt="Build Analysis" src="https://github.com/mfem/mfem/actions/workflows/mfem-analysis.yml/badge.svg?branch=master"></a>
+<a href="https://github.com/mfem/mfem/actions?query=workflow%3Abuilds-and-tests+branch%3Amaster"><img alt="Builds and Tests" src="https://github.com/mfem/mfem/actions/workflows/builds-and-tests.yml/badge.svg?branch=master"></a>
 <a href="https://ci.appveyor.com/project/mfem/mfem"><img alt="Build Status" src="https://ci.appveyor.com/api/projects/status/19non9sqm6msi2wy?svg=true"></a>
 <a href="https://docs.mfem.org/html/index.html"><img alt="Doxygen" src="https://img.shields.io/badge/code-documented-brightgreen.svg"></a>
 </p>
@@ -27,7 +27,7 @@ in the MFEM community, you agree to abide by its rules.
 If you plan on contributing to MFEM, consider reviewing the
 [issue tracker](https://github.com/mfem/mfem/issues) first to check if a thread
 already exists for your desired feature or the bug you ran into. Use a pull
-request (PR) toward the `mfem:main` branch to propose your contribution. If
+request (PR) toward the `mfem:master` branch to propose your contribution. If
 you are planning significant code changes or have questions, you may want to
 open an [issue](https://github.com/mfem/mfem/issues) before issuing a PR. In
 addition to technical contributions, we are also interested in your results and
@@ -47,7 +47,7 @@ back to them before issuing pull requests:
  - [Pull Requests](#pull-requests)
  - [MFEM PR Rules](#mfem-pr-rules)
  - [Pull Request Checklist](#pull-request-checklist)
-  - [Main/Next Workflow](#mainnext-workflow)
+  - [Master/Next Workflow](#masternext-workflow)
  - [Releases](#releases)
  - [Release Checklist](#release-checklist)
 - [LLNL Workflow](#llnl-workflow)
@@ -66,12 +66,12 @@ Origin](#developers-certificate-of-origin-11) at the end of this file.*
 ## Quick Summary

 - We encourage you to [join the MFEM organization](#mfem-organization) and create
-  development branches off `mfem:main`.
+  development branches off `mfem:master`.
 - Please follow the [developer guidelines](#developer-guidelines), in particular
  with regards to documentation and code styling.
 - Please do not commit large/binary files to the central repository (use a fork
  instead).
- Pull requests should be issued toward `mfem:main`. Make sure
+- Pull requests should be issued toward `mfem:master`. Make sure
  to check the items off the [Pull Request Checklist](#pull-request-checklist) and
  follow the [MFEM PR Rules](#mfem-pr-rules).
 - When your contribution is fully working and ready to be reviewed, add
@@ -81,8 +81,8 @@ Origin](#developers-certificate-of-origin-11) at the end of this file.*
 - The reviewers have 3 weeks to evaluate the PR and work with the author to
  fix issues and implement improvements.
 - During review there should be no force pushes/rewriting history in the branch.
- After approval, MFEM developers merge the PR manually in the [mfem:next branch](#mainnext-workflow).
- After a week of testing in `mfem:next`, the original PR is merged in `mfem:main`.
+- After approval, MFEM developers merge the PR manually in the [mfem:next branch](#masternext-workflow).
+- After a week of testing in `mfem:next`, the original PR is merged in `mfem:master`.
 - We use [milestones](https://github.com/mfem/mfem/milestones) to coordinate the
  work on different PRs toward a release.
 - Don't hesitate to [contact us](#contact-information) if you have any questions.
@@ -289,7 +289,7 @@ Before you can start, you need a GitHub account, here are a few suggestions:
  author, is willing to work on it and be its champion.

 - The author creates a branch for the new feature (with suffix `-dev`), off
-  the `main` branch, or another existing feature branch, for example:
+  the `master` branch, or another existing feature branch, for example:

  ```
  # Clone assuming you have setup your ssh keys on GitHub:
@@ -298,8 +298,8 @@ Before you can start, you need a GitHub account, here are a few suggestions:
  # Alternatively, clone using the "https" protocol:
  git clone https://github.com/mfem/mfem.git

-  # Create a new feature branch starting from "main":
-  git checkout main
+  # Create a new feature branch starting from "master":
+  git checkout master
  git pull
  git checkout -b feature-dev

@@ -375,7 +375,7 @@ Before you can start, you need a GitHub account, here are a few suggestions:
 ### Pull Requests

 - When your branch is ready for other developers to review / comment on
-  the code, create a pull request towards `mfem:main`.
+  the code, create a pull request towards `mfem:master`.

 - Pull request typically have titles like:

@@ -411,7 +411,7 @@ Before you can start, you need a GitHub account, here are a few suggestions:
  assigned, the PR is considered under review. To help with the review process
  there should be no force pushes/rewriting history in the branch.

- After approval, the PR is [tested](#mainnext-workflow) for a week with
+- After approval, the PR is [tested](#masternext-workflow) for a week with
  other approved PRs in the `mfem:next` branch.

 - Consider manually running the tests in `tests/scripts` before merging in
@@ -447,7 +447,7 @@ The Pull Request (PR) approval process in MFEM is similar to the approval of pap

 3. A PR can be (manually) merged in the *next* branch only if 2 of the assigned reviewers have approved it and it has passed internal testing. This merge can be performed by any of the assigned reviewers or by any of the editors.

-4. A PR can be merged in the *main* branch only if it has been tested successfully for a week in *next* and an editor has (optionally) taken a final look. This merge can be performed only by one of the editors.
+4. A PR can be merged in the *master* branch only if it has been tested successfully for a week in *next* and an editor has (optionally) taken a final look. This merge can be performed only by one of the editors.

 #### Responsibilities of Editors

@@ -470,7 +470,7 @@ The current list of MFEM editors is:

 5. To remind the reviewers about timely completion of their review.

-6. To take a final look and complete the PR merge in *main*. The final look step is optional and shouldn't take more than 3 days.
+6. To take a final look and complete the PR merge in *master*. The final look step is optional and shouldn't take more than 3 days.

 7. The assignment of bugfixes should be expedited proportional to their importance, e.g. in some cases the editor can assign much shorter review window.

@@ -492,7 +492,7 @@ Everyone on the MFEM team can be asked to serve as a reviewer on a PR in their a

 5. To test the PR thoroughly before merging in *next*. The PR author is also encouraged to perform testing and inform the reviewers about the results.

-6. To monitor the PR impact on the testing in the *next* branch and alert the editors that the PR is ready for merging in *main*.
+6. To monitor the PR impact on the testing in the *next* branch and alert the editors that the PR is ready for merging in *master*.

 7. The review of bugfixes should be expedited proportional to their importance. The review window can be much less than three weeks in such cases.

@@ -582,11 +582,11 @@ Before a PR can be merged, it should satisfy the following:
   - [ ] Update internal tests to include the new features.


-### Main/Next Workflow
+### Master/Next Workflow

-MFEM uses a `main`/`next`-branch workflow as described below:
+MFEM uses a `master`/`next`-branch workflow as described below:

- The `main` branch should always be of release quality and changes should not
+- The `master` branch should always be of release quality and changes should not
  be merged until they have been fully tested. This branch is protected, and
  changes can only be made through pull requests.

@@ -613,20 +613,20 @@ MFEM uses a `main`/`next`-branch workflow as described below:

 - After a week of testing in `next` (excluding bugfixes), both on GitHub, as
  well as [internally](#tests-at-llnl) at LLNL, the original PR is merged into
-  `main` (provided there are no issues).
+  `master` (provided there are no issues).

 - After the merge, the feature branch is deleted (unless it is a long-term
  project with periodic PRs).

 - The `next` branch is used just for integrated testing of all PRs approved for
-  merging into `main` to verify that each works individually and that all of
+  merging into `master` to verify that each works individually and that all of
  them work as a group. This branch can be discarded at any time, though we
  typically do that only at the end of a [release cycle](#releases).


 ### Releases

- Releases are just tags in the `main` branch, e.g. https://github.com/mfem/mfem/releases/tag/v3.3.2,
+- Releases are just tags in the `master` branch, e.g. https://github.com/mfem/mfem/releases/tag/v3.3.2,
  and have a version that ends in an even "patch" number, e.g. `v3.2.2` or
  `v3.4` (by convention `v3.4` is the same as `v3.4.0`.)  Between releases, the
  version ends in an odd "patch" number, e.g. `v3.3.3`.
@@ -679,19 +679,19 @@ MFEM uses a `main`/`next`-branch workflow as described below:

 ### Mirroring on Bitbucket

- The GitHub `main` and `next` branches are mirrored to the LLNL institutional
-  Bitbucket repository as `gh-main` and `gh-next`.
+- The GitHub `master` and `next` branches are mirrored to the LLNL institutional
+  Bitbucket repository as `gh-master` and `gh-next`.

- `gh-main` is merged into LLNL's internal `main` through pull requests; write
-  permissions to `main` are restricted to ensure this is the only way in which it
+- `gh-master` is merged into LLNL's internal `master` through pull requests; write
+  permissions to `master` are restricted to ensure this is the only way in which it
  gets updated.

 - We never push directly from LLNL to GitHub.

 - Versions of the code on LLNL's internal server, from most to least stable:
  - MFEM official release on mfem.org -- Most stable, tested in many apps.
-  - `mfem:main` -- Recent development version, guaranteed to work.
-  - `mfem:gh-main` -- Stable development version, passed testing, you can use
+  - `mfem:master` -- Recent development version, guaranteed to work.
+  - `mfem:gh-master` -- Stable development version, passed testing, you can use
     it to build your code between releases.
  - `mfem:gh-next` -- Bleeding-edge development version, may be broken, use at
     your own risk.
@@ -721,14 +721,14 @@ directory.


 ### Linux and Mac smoke tests
-We use GitHub Actions to drive the default tests on the `main` and `next`
+We use GitHub Actions to drive the default tests on the `master` and `next`
 branches. See the `.github/workflows` files and the logs at
 [https://github.com/mfem/mfem/actions](https://github.com/mfem/mfem/actions).

 Testing using GitHub Actions should be kept lightweight, as there is a time
 constraint on jobs. Two virtual machines are configured - Mac (OS X) and Linux.

- Tests on the `main` branch are triggered whenever a PR is issued on this branch.
+- Tests on the `master` branch are triggered whenever a PR is issued on this branch.
 - Tests on the `next` branch are currently scheduled to run each night.


@@ -744,7 +744,7 @@ and debug build is performed with a simple run of `ex1` to verify the executable

 ### Tests at LLNL

- We mirror the `main` and `next` branches internally (to `gh-main` and
+- We mirror the `master` and `next` branches internally (to `gh-master` and
  `gh-next`) and run longer nightly tests via cron. On the weekends, a more
  extensive test is run which extracts and executes all the different sample
  runs from each example and most miniapps.
@@ -33,8 +33,8 @@ RUN mkdir -p /opt/mfem-env \
 RUN cd /opt/mfem-env && \
    . /opt/spack/share/spack/setup-env.sh && \
    spack env activate . && \
-    spack develop --path /code mfem@main+examples+miniapps && \
-    spack add mfem@main+examples+miniapps # && \
+    spack develop --path /code mfem@master+examples+miniapps && \
+    spack add mfem@master+examples+miniapps && \
    spack install

 # ensure mfem always on various paths
@@ -1,11 +1,11 @@
 spack:
-  specs: [mfem@main+examples+miniapps]
-  view:
+  specs: [mfem@master+examples+miniapps]
+  view: 
    mfem:
-      root: /opt/mfem-view
+      root: /opt/mfem-view 
      link_type: copy
  concretization: together
  develop:
    mfem:
      path: /code
-      spec: mfem@main+examples+miniapps
+      spec: mfem@master+examples+miniapps
@@ -115,7 +115,7 @@ fi

 # branch-history
 if [[ "${option}" == "--history" || "${option}" == "" ]]; then
-  git fetch origin main:main
+  git fetch origin master:master
  cd tests/scripts
  if ! ./runtest branch-history; then code=1; fi
  cd -
@@ -1,6 +1,6 @@
 # Jupyter Notebooks using xeus-cling

-[![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/mfem/mfem/main?filepath=examples%2Fjupyter%2Fex.ipynb)
+[![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/mfem/mfem/master?filepath=examples%2Fjupyter%2Fex.ipynb)

 [xeus-cling](https://github.com/jupyter-xeus/xeus-cling) is a C++ Jupyter Kernel based on [cling](https://github.com/root-project/cling),
 which can be used to create interactive C++ MFEM and GLVis notebooks.
@@ -436,7 +436,7 @@ void NonlinearForm::Mult(const Vector &x, Vector &y) const
   // In parallel, the result is in 'py' which is an alias for 'aux2'.
 }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

   void Update(FiniteElementSpace *nfes = NULL) override;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

   /** @brief Copy constructor for a ParCSR matrix which creates a deep copy of
       structure and data from @a P. */
@@ -92,18 +92,6 @@ struct LpReducer
   }
 };

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

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

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

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

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

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

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

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

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

   T_solver.Mult(z, du_dt);
   du_dt.SetSubVector(ess_tdof_list, 0.0);
@@ -89,7 +89,7 @@ The currently available test scripts are:
 - gitignore: Checks for compliance with MFEM's .gitignore rules

 - branch-history: Checks if the current branch history (the commits that will be
-  merged in main) contains unusually large files, unusually large number of
+  merged in master) contains unusually large files, unusually large number of
  changes in a commit, unusually large number of commits, etc.

 Adding new test scripts is easy and we welcome contributions from users and
@@ -27,7 +27,7 @@ sub usage {
          -h|--help   prints this usage information and exits

   This script checks if the current branch history, defined as the commits that
-   will be merged in main (those shown in a GitHub PR), contains unusually
+   will be merged in master (those shown in a GitHub PR), contains unusually
   large files, unusually large number of changes in a commit, unusually large
   number of commits, etc.

@@ -72,7 +72,7 @@ my $max_branch_kb = 1000;
 my $status = 0; # Return code

 # Get SHA hash of all commits in this branch
-my @commits = split /\n/, `git log --pretty=format:%H main..$branch`;
+my @commits = split /\n/, `git log --pretty=format:%H master..$branch`;

 # Check if total number of commits in this branch exceeds the maximum allowable
 my $ncommits = scalar @commits;
@@ -22,7 +22,6 @@ using namespace mfem;

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

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

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

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

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

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

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

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

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

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

 TEST_CASE("Reduce ArgMinMax", "[Reduction],[GPU]")
 {
-   Array<MinMaxLocScalar<double, int>> workspace;
   Array<double> a(1000);
   auto hptr = a.HostReadWrite();
   for (int i = 0; i < a.Size(); ++i)
@@ -383,7 +373,7 @@ TEST_CASE("Reduce ArgMinMax", "[Reduction],[GPU]")
            r.max_loc = i;
         }
      },
-      ArgMinMaxReducer<double, int> {}, use_dev, workspace);
+      ArgMinMaxReducer<double, int> {}, use_dev);
      CAPTURE(use_dev);
      REQUIRE(res.min_val == -10);
      REQUIRE(res.min_loc >= 0);
Author	SHA1	Message	Date
Will Pazner	552d6857cb	Add GPU scan support for Array<T>::PartialSum	2025-10-21 09:57:37 -07:00
Will Pazner	a37d46e917	Fix comments in scan.hpp	2025-10-21 09:57:37 -07:00
Will Pazner	4acdb072b6	Add more device support to Array<T>	2025-10-21 09:57:37 -07:00
Will Pazner	9dbb184537	Remove need to explicitly pass workspace array to reducers	2025-10-21 09:57:37 -07:00
Will Pazner	d67762a1c9	Move contents of array.cpp to array.hpp Remove explicit template instantiations	2025-10-20 15:38:54 -07:00
Tzanio Kolev	cd4e583f9f	Merge pull request #4659 from mfem/najlkin/parform-impro Improvements of Par(Mixed)BilinearForm and Par(Block)NonlinearForm	2025-10-18 10:48:44 -07:00
Tzanio Kolev	ee94776558	Merge pull request #5066 from mfem/najlkin/fix-nc-blknlform [BUG] Non-conforming meshes in BlockNonlinearForm	2025-10-18 10:48:25 -07:00
Tzanio Kolev	af478afd00	Merge branch 'master' into najlkin/fix-nc-blknlform	2025-10-15 16:11:05 -07:00
Jan Nikl	e13d1a1d53	Fixed support of non-conforming meshes in BlockNonlinearForm.	2025-10-10 13:32:53 -07:00
Jan Nikl	b7a0b2cf9a	Renamed ParallelAssembleInternal(Matrix) and added more documentation.	2025-10-07 12:28:31 -07:00
Tzanio Kolev	bce6e2ca76	Merge branch 'master' into najlkin/parform-impro	2025-10-05 13:25:11 -07:00
Tzanio Kolev	40d1550fd6	Merge branch 'master' into najlkin/parform-impro	2025-09-17 03:31:41 -07:00
Jan Nikl	193f8a6801	Partially reverted const modifiers in HypreParMatrix::Copy(Bool)CSR().	2025-09-16 16:19:03 -07:00
Jan Nikl	110720dd04	Fixed documentation of Par(Mixed)BilinearForm::KeepNbrBlock().	2025-09-16 15:39:35 -07:00
Jan Nikl	3fd335c77b	Fixed usage of EliminateVDofsInRHS().	2025-09-16 15:36:04 -07:00
Jan Nikl	058fdaae3f	Implemented gradient of ParBlockNonlinearForm with shared face contributions.	2025-06-11 16:23:42 -07:00
Jan Nikl	9a92e4875b	Implemented Mult of ParBlockNonlinearForm with shared face contributions.	2025-06-11 16:22:16 -07:00
Jan Nikl	9dab032bd0	Merge branch 'master' into najlkin/parform-impro	2025-04-24 15:40:44 -07:00
Jan Nikl	cdfe8102ae	Added gradient of ParNonlinearForm with face integrators.	2025-03-18 16:44:51 -07:00
Jan Nikl	5b82bf0328	Revert "WIP: Added support for trace face integrators in ParMixedBilinearForm." This reverts commit `323ee572b6`.	2025-03-06 06:00:07 -08:00
Tzanio Kolev	41b65d6333	Merge branch 'master' into najlkin/parform-impro	2025-02-04 14:55:19 -08:00
Jan Nikl	2534d2207d	Fixed name of ParallelEliminateTrialEssentialBC().	2025-01-09 10:49:10 -08:00
Jan Nikl	260b817b3c	Merge branch 'master' into najlkin/parform-impro	2025-01-09 10:01:52 -08:00
Jan Nikl	613d5dd826	Fixed constness in some HyperParMatrix constructors.	2025-01-08 17:56:42 -08:00
Jan Nikl	323ee572b6	WIP: Added support for trace face integrators in ParMixedBilinearForm.	2025-01-08 17:56:05 -08:00
Jan Nikl	56ff5ac5bb	Added support for interior face integrators to ParMixedBilinearForm.	2025-01-08 17:55:42 -08:00
Jan Nikl	e1a06bd6c8	Added methods to Par(Mixed)BilinearForm for elimination of essential BCs.	2025-01-08 17:55:08 -08:00
Jan Nikl	9acae54669	Extended ParMixedBilinearForm methods for parallel assembly.	2025-01-08 17:54:28 -08:00
Jan Nikl	7d92e22a45	Added ParallelAssembleInternal() method to ParBilinearForm.	2025-01-08 17:53:21 -08:00