Merge branch 'master' of https://github.com/mfem/mfem into mesh-dist

initial commit
Merge pull request #5045 from mfem/stefanozampini/petsc-requires-hypre
2026-06-09 10:00:05 -07:00 · 2026-06-09 09:59:57 -07:00 · 2026-06-06 16:52:11 -07:00 · 2026-06-05 16:17:51 -07:00 · 2026-06-05 16:13:18 -07:00 · 2026-06-05 16:10:03 -07:00
35 changed files with 1466 additions and 461 deletions
@@ -0,0 +1,42 @@
+# MFEM Pull Request Review Agent Guide
+
+## Purpose and scope
+Review MFEM PRs for correctness, maintainability, performance, portability, test coverage, and MFEM consistency. Use the diff and PR context; reference source files, tests, and CI results when available. Follow `CONTRIBUTING.md`, especially Developer Guidelines, PR rules, checklist, and testing.
+
+## Critical review pillars
+- Correctness and numerical behavior
+- API and user-facing impact
+- Performance implications
+- Maintainability and portability
+
+## Review workflow
+1. Read the PR description, linked issues, and intended behavior.
+2. Inspect the diff before commenting.
+3. Identify affected MFEM components, examples, tests, build or docs changes, and downstream APIs.
+4. Analyze the code against the critical review pillars.
+5. Compare the change against nearby code and MFEM patterns; flag unmotivated deviations.
+6. Check whether tests and documentation were updated appropriately.
+7. Review CI results and suggest actions.
+8. Produce a structured review with prioritized findings.
+9. Always limit conclusions to available evidence.
+
+## MFEM-specific review checklist
+- Component-aware scope: identify the touched subsystem (FEM, solvers, preconditioners, linear algebra, mesh, examples, miniapps, build, or docs) and assess its impact against the review pillars.
+- Numerical and algorithmic behavior: assess issues in convergence, stability, tolerances, precision, iteration limits, and failure handling. If clear opportunities exist to improve the algorithmic approach, call them out with expected impact.
+- API and user-facing impact: assess backward compatibility, user-visible behavior and default changes, migration impact, deprecations, and whether documentation clearly explains user-facing API changes.
+- Data structure and memory semantics: assess ownership, lifetime, aliasing, container behavior, and device-host synchronization.
+- Parallel and serial behavior: assess whether the change preserves equivalent semantics in serial and parallel modes where applicable; if logic is currently mode-specific, check whether extension to the other mode is straightforward (clear abstractions, no hard-wired assumptions), document constraints, and call out expected behavior differences explicitly.
+- Backend and portability impact: assess likely cross-backend risks in CPU, CUDA, HIP, OCCA, RAJA, partial assembly, fallback paths, compiler compatibility, and platform assumptions.
+- Build, dependency, and configuration impact: assess CMake or make changes, optional dependency behavior, and feature-flag interactions.
+- Tests and docs alignment: check available regression or unit coverage evidence for changed behavior, and ensure docs are updated for new flags, APIs, options, or behavior changes.
+- MFEM developer-guideline fit: keep code lean, simple, general, logically separated, and portable; suggest C++17 improvements when they clearly improve safety, clarity, or maintainability.
+- New source files, examples, or miniapps: if a PR adds source/header files, verify they are properly wired into the relevant `makefile` and `CMakeLists.txt`, referenced in docs where applicable (including `doc/CodeDocumentation.dox`), and added to top-level `.gitignore` only when generated artifacts require it.
+- Changelog: verify `CHANGELOG` is updated if the PR introduces significant new features or user-facing changes.
+- MFEM conventions: use `real_t`; use `mfem::out`/`mfem::err` instead of `std::cout`/`std::cerr` in library code; flag large/binary files; if AI assistance is apparent but undisclosed, suggest following `CONTRIBUTING.md`.
+- Edge cases: if the PR touches complex or error-prone areas, suggest additional tests for edge cases, failure modes, and parallel behavior.
+
+## Commenting guidelines
+- Keep comments concise, actionable, and grounded in the diff.
+- Focus on correctness, behavior changes, and user impact over style nits.
+- Be professional, concise, collaborative, technically precise, and avoid unsupported assumptions.
+
@@ -15,6 +15,10 @@ Version 4.9.1 (development)

 Discretization improvements
 ---------------------------
+- Added NVIDIA cuDSS library interface. Implementation examples have been 
+  added to ex1 and ex1p. See https://developer.nvidia.com/cudss for more 
+  details. Supported versions >= 0.6.0.
+
 - Extend FindPointsGSLIB to support surface meshes.

 - Replaced legacy simplex quadrature rules with symmetric positive-weight
@@ -39,11 +43,24 @@ Meshing improvements
 - Improved support for 1D NURBS meshes with variable order, including using
  the patches construct for 1D NURBS meshes.

+- Added the option to include material interfaces (faces separating elements
+  with different element attributes) as additional boundary elements, for
+  parallel visualization, e.g. with GLVis. This is supported by both the Print
+  and PrintAsOne methods of ParMesh. See ParMesh::SetPrintInterfaces().
+
 New and updated examples and miniapps
 -------------------------------------
 - Electromagnetics/lorentz miniapp has been updated to leverage the ParticleSet
  capability.

+Miscellaneous
+-------------
+- Fixed signed DOF handling in parallel grid-function reading (read constructor)
+  and saving via ParGridFunction::SaveAsOne(). Simplified the process of
+  applying the DOF signs by using the new method ApplyDofSigns() in class
+  ParFiniteElementSpace -- the method will return immediately if no sign flips
+  are needed.
+

 Version 4.9, released on Dec 11, 2025
 =====================================
@@ -433,6 +433,15 @@ if (MFEM_USE_STRUMPACK)
  endif()
 endif()

+# cuDSS can only be enabled in CUDA
+if (MFEM_USE_CUDSS)
+  if (MFEM_USE_CUDA)
+    find_package(CUDSS REQUIRED)
+  else()
+    message(FATAL_ERROR " *** cuDSS requires that CUDA be enabled.")
+  endif()
+endif()
+
 # GnuTLS
 if (MFEM_USE_GNUTLS)
  find_package(_GnuTLS REQUIRED)
@@ -631,7 +640,7 @@ find_package(Threads REQUIRED)
 set(MFEM_TPLS OPENMP HYPRE LAPACK BLAS SuperLUDist STRUMPACK METIS SuiteSparse
    SUNDIALS PETSC SLEPC MUMPS AXOM FMS CONDUIT Ginkgo GNUTLS GSLIB HDF5
    NETCDF MPFR PUMI HIOP POSIXCLOCKS MFEMBacktrace ZLIB OCCA CEED RAJA UMPIRE
-    ADIOS2 MKL_CPARDISO MKL_PARDISO AMGX MAGMA CUSPARSE CUBLAS CALIPER CODIPACK
+    ADIOS2 MKL_CPARDISO MKL_PARDISO AMGX MAGMA CUSPARSE CUBLAS CUDSS CALIPER CODIPACK
    BENCHMARK PARELAG TRIBOL MPI_CXX HIP HIPBLAS HIPSPARSE MOONOLITH BLITZ
    ALGOIM ENZYME CUDA::cudart)

@@ -3,12 +3,13 @@
 </p>

 <p align="center">
-<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://github.com/mfem/mfem/blob/master/LICENSE"><img alt="License" src="https://img.shields.io/badge/License-BSD-blue.svg"></a>
+<a href="https://github.com/mfem/mfem/releases/latest"><img alt="GitHub release" src="https://img.shields.io/github/v/release/mfem/mfem"></a>
+<a href="https://github.com/mfem/mfem/actions/workflows/repo-check.yml?query=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/workflows/mfem-analysis.yml?query=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/workflows/builds-and-tests.yml?query=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>
+<a href="https://docs.mfem.org/html/index.html"><img alt="Documentation" src="https://img.shields.io/badge/code-documented-brightgreen.svg"></a>
 </p>


@@ -84,7 +85,7 @@ Origin](#developers-certificate-of-origin-11) at the end of this file.*
  follow the [MFEM PR Rules](#mfem-pr-rules).
 - When your contribution is fully working and ready to be reviewed, add
  the `ready-for-review` label.
- PRs are treated similarly to journal submission with an "editor" assigning two
+- PRs are treated similarly to journal submission, with an "editor" assigning two
  reviewers to evaluate the changes.
 - The reviewers have 3 weeks to evaluate the PR and work with the author to
  fix issues and implement improvements.
@@ -125,7 +126,7 @@ The MFEM source code has the following structure:
  │   ├── petsc
  │   ├── pumi
  │   ├── sundials
-  |   └── superlu
+  │   └── superlu
  ├── fem
  │   ├── ceed
  │   ├── dfem
@@ -137,10 +138,6 @@ The MFEM source code has the following structure:
  │   ├── moonolith
  │   ├── qinterp
  │   └── tmop
-  │   |   ├── assemble
-  │   |   ├── metrics
-  │   |   ├── mult
-  │   |   └── tools
  ├── general
  ├── linalg
  │   ├── batched
@@ -153,11 +150,10 @@ The MFEM source code has the following structure:
  │   ├── common
  │   ├── contact
  │   ├── dfem
+  │   ├── diag-smoothers
  │   ├── dpg
  │   ├── electromagnetics
  │   ├── fluids
-  │   │   ├── navier
-  │   │   └── schrodinger-flow
  │   ├── gslib
  │   ├── hdiv-linear-solver
  │   ├── hooke
@@ -167,6 +163,7 @@ The MFEM source code has the following structure:
  │   ├── nurbs
  │   ├── parelag
  │   ├── performance
+  │   ├── plasma
  │   ├── shifted
  │   ├── solvers
  │   ├── spde
@@ -197,15 +194,15 @@ respectively.

 - The main finite element classes are:
  + [`FiniteElement`](https://docs.mfem.org/html/classmfem_1_1FiniteElement.html)
-  + [`FiniteElementCollection`](https://docs.mfem.org/html/classmfem_1_1FiniteElement.html)
+  + [`FiniteElementCollection`](https://docs.mfem.org/html/classmfem_1_1FiniteElementCollection.html)
  + [`FiniteElementSpace`](https://docs.mfem.org/html/classmfem_1_1FiniteElementSpace.html)
  + [`GridFunction`](https://docs.mfem.org/html/classmfem_1_1GridFunction.html)
  + [`BilinearFormIntegrator`](https://docs.mfem.org/html/classmfem_1_1BilinearFormIntegrator.html) and [`LinearFormIntegrator`](https://docs.mfem.org/html/classmfem_1_1LinearFormIntegrator.html)
-  + [`LinearForm`](https://docs.mfem.org/html/classmfem_1_1LinearFormIntegrator.html), [`BilinearForm`](https://docs.mfem.org/html/classmfem_1_1BilinearForm.html) and [`MixedBilinearForm`](https://docs.mfem.org/html/classmfem_1_1MixedBilinearForm.html)
+  + [`LinearForm`](https://docs.mfem.org/html/classmfem_1_1LinearForm.html), [`BilinearForm`](https://docs.mfem.org/html/classmfem_1_1BilinearForm.html) and [`MixedBilinearForm`](https://docs.mfem.org/html/classmfem_1_1MixedBilinearForm.html)

 - The main linear algebra classes and sources are
  + [`Operator`](https://docs.mfem.org/html/classmfem_1_1Operator.html) and [`BilinearForm`](https://docs.mfem.org/html/classmfem_1_1BilinearForm.html)
-  + [`Vector`](https://docs.mfem.org/html/classmfem_1_1BilinearForm.html) and [`LinearForm`](https://docs.mfem.org/html/classmfem_1_1LinearForm.html)
+  + [`Vector`](https://docs.mfem.org/html/classmfem_1_1Vector.html) and [`LinearForm`](https://docs.mfem.org/html/classmfem_1_1LinearForm.html)
  + [`DenseMatrix`](https://docs.mfem.org/html/classmfem_1_1DenseMatrix.html) and [`SparseMatrix`](https://docs.mfem.org/html/classmfem_1_1SparseMatrix.html)
  + Sparse [smoothers](https://docs.mfem.org/html/sparsesmoothers_8hpp.html) and linear [solvers](https://docs.mfem.org/html/solvers_8hpp.html)

@@ -217,8 +214,8 @@ shared geometric entities between different tasks. The parallel source files
 have a `p` prefix, e.g. `pmesh.cpp` vs. the serial `mesh.cpp`.

 - The main parallel classes are
-  + [`ParMesh`](https://docs.mfem.org/html/solvers_8hpp.html)
-  + [`ParNCMesh`](https://docs.mfem.org/html/classmfem_1_1ParMesh.html)
+  + [`ParMesh`](https://docs.mfem.org/html/classmfem_1_1ParMesh.html)
+  + [`ParNCMesh`](https://docs.mfem.org/html/classmfem_1_1ParNCMesh.html)
  + [`ParFiniteElementSpace`](https://docs.mfem.org/html/classmfem_1_1ParFiniteElementSpace.html)
  + [`ParGridFunction`](https://docs.mfem.org/html/classmfem_1_1ParGridFunction.html)
  + [`ParBilinearForm`](https://docs.mfem.org/html/classmfem_1_1ParBilinearForm.html) and [`ParLinearForm`](https://docs.mfem.org/html/classmfem_1_1ParLinearForm.html)
@@ -228,14 +225,14 @@ have a `p` prefix, e.g. `pmesh.cpp` vs. the serial `mesh.cpp`.
 #### GPU and general device support

 GPU and multi-core CPU support is based on device kernels supporting different
-backends (CUDA, OCCA, RAJA, OpenMP, etc.) and an internal lightweight
+backends (CUDA, HIP, OCCA, RAJA, OpenMP, etc.) and an internal lightweight
 device/host memory manager.

 - The main device-relevant classes and sources are:
  + [`Device`](https://docs.mfem.org/html/device_8hpp.html)
  + [`MemoryManager`](https://docs.mfem.org/html/mem_manager_8hpp.html)
  + the [`mfem::forall`](https://docs.mfem.org/html/forall_8hpp.html) function
-  + the [`cuda.hpp`](https://docs.mfem.org/html/cuda_8hpp.html) and [`occa.hpp`](https://docs.mfem.org/html/occa_8hpp.html) files
+  + the [`cuda.hpp`](https://docs.mfem.org/html/cuda_8hpp.html), [`hip.hpp`](https://docs.mfem.org/html/hip_8hpp.html) and [`occa.hpp`](https://docs.mfem.org/html/occa_8hpp.html) files

 #### Utilities, building and documentation
 - The `general/` directory contains C++ classes that serve as utilities for
@@ -249,8 +246,8 @@ device/host memory manager.
 - `examples` and `miniapps` respectively gather simple and more fully-featured
  demonstrations of the usage on MFEM. They both rely on `data/` for the
  collection of meshes.
- The `tests/` directory contains a unit test suite and will later contain more
-  tests that run example codes.
+- The `tests/` directory contains a unit test suite, additional tests, and
+  benchmarks.

 See also the [code overview](https://mfem.org/code-overview/) section on the MFEM
 website.
@@ -284,8 +281,8 @@ Before you can start, you need a GitHub account, here are a few suggestions:
  the top of https://github.com/mfem.
 - Consider making your membership public by going to https://github.com/orgs/mfem/people
  and clicking on the organization visibility drop box next to your name.
- Project discussions and announcements will be posted at
-  https://github.com/orgs/mfem/teams/everyone.
+- Project discussions and announcements will be posted at https://github.com/orgs/mfem/discussions,
+  tagging the `@mfem/everyone` team when appropriate.

 #### Structure
 - The MFEM source code is in the [mfem](https://github.com/mfem/mfem)
@@ -345,11 +342,12 @@ Before you can start, you need a GitHub account, here are a few suggestions:
  - Well-designed simple code is frequently more general and powerful.
  - Lean code base is easier to understand by new collaborators.
  - New features should be added only if they are necessary or generally useful.
-  - Introduction of language constructions not currently used in MFEM should be
+  - Introduction of language constructs not currently used in MFEM should be
    justified and generally avoided (to maintain portability to various systems
    and compilers, including early access hardware).
-  - We prefer basic C++ and the C++03 standard, to keep the code readable by
-    a large audience and to make sure it compiles anywhere.
+  - We prefer basic C++. Use C++17 features judiciously, prioritizing readability,
+    consistency with existing MFEM code, and portability to different systems,
+    compilers and device backends.

 - *Keep the code general and reasonably efficient*
  - The main goal is fast prototyping for research and application development.
@@ -392,7 +390,7 @@ Before you can start, you need a GitHub account, here are a few suggestions:
 - When your branch is ready for other developers to review / comment on
  the code, create a pull request towards `mfem:master`.

- Pull request typically have titles like:
+- Pull requests typically have titles like:

     `Description [new-feature-dev]`

@@ -413,12 +411,12 @@ Before you can start, you need a GitHub account, here are a few suggestions:
 - Add a description, appropriate labels and assign yourself to the PR. The MFEM
  team will add reviewers as appropriate.

- List outstanding TODO items in the description, see PR #222 for an example.
+- List outstanding TODO items in the description.

 - When your contribution is fully working and ready to be reviewed, add
-  the `ready-for-review` label.
+  or request the `ready-for-review` label.

- PRs are treated similarly to journal submission with an "editor" assigning
+- PRs are treated similarly to journal submission, with an "editor" assigning
  two reviewers to evaluate the changes. The reviewers have 3 weeks to evaluate
  the PR and work with the author to implement improvements and fix issues.

@@ -444,7 +442,7 @@ Before you can start, you need a GitHub account, here are a few suggestions:
  checks in GitHub Actions enforce MFEM-specific rules which are explained in
  the error messages and the `tests/scripts` directory.

- Also note that the tests `branch-history` and `repos-checks` found in GitHub
+- Also note that the tests `branch-history` and `repo-check` found in GitHub
  Actions can be triggered automatically before each push using git hooks. See
  the [git hooks README](config/githooks/README.md) for a detailed explanation.

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

 3. To ensure the quality of the PR by making sure that the code adheres to the [Developer Guidelines](#developer-guidelines), e.g. all methods, data members, and functions have documentation, including data ownership and lifetime, new examples/miniapps have a corresponding PR in mfem/web, major features have `CHANGELOG` entries, etc.

-3. To seek help from the editors in case of difficulties.
+4. To seek help from the editors in case of difficulties.

-4. To complete the review in a timely manner: 3 weeks from assignment.
+5. To complete the review in a timely manner: 3 weeks from assignment.

-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 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 *master*.
+7. 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.
+8. The review of bugfixes should be expedited proportional to their importance. The review window can be much less than three weeks in such cases.

 #### Responsibilities of Authors

@@ -535,30 +533,30 @@ Before a PR can be merged, it should satisfy the following:
 - [ ] Code builds.
 - [ ] Code passes `make style`.
 - [ ] Update `CHANGELOG`:
-    - [ ] Is this a new feature users need to be aware of? New or updated example or miniapp?
-    - [ ] Does it make sense to create a new section in the `CHANGELOG` to group with other related features?
+   - [ ] Is this a new feature users need to be aware of? New or updated example or miniapp?
+   - [ ] Does it make sense to create a new section in the `CHANGELOG` to group with other related features?
 - [ ] Update `INSTALL`:
-    - [ ] Had a new optional library been added? If so, what range of versions of this library are required? (*Make sure the external library is compatible with our BSD license, e.g. it is not licensed under GPL!*)
-    - [ ] Have the version ranges for any required or optional libraries changed?
-    - [ ] Does `make` or `cmake` have a new target?
-    - [ ] Did the requirements or the installation process change? *(rare)*
+   - [ ] Has a new optional library been added? If so, what range of versions of this library are required? (*Make sure the external library is compatible with our BSD license, e.g. it is not licensed under GPL!*)
+   - [ ] Have the version ranges for any required or optional libraries changed?
+   - [ ] Does `make` or `cmake` have a new target?
+   - [ ] Did the requirements or the installation process change? *(rare)*
 - [ ] Update continuous integration server configurations if necessary (e.g. with new version requirements for each of MFEM's dependencies)
-    - [ ] `.github`
-    - [ ] `.appveyor.yml`
+   - [ ] `.github`
+   - [ ] `.appveyor.yml`
 - [ ] Update `.gitignore`:
-    - [ ] Check if `make distclean; git status` shows any files that were generated from the source by the project (not an IDE) but we don't want to track in the repository.
-    - [ ] Add new patterns (just for the new files above) and re-run the above test.
+   - [ ] Check if `make distclean; git status` shows any files that were generated from the source by the project (not an IDE) but we don't want to track in the repository.
+   - [ ] Add new patterns (just for the new files above) and re-run the above test.
 - [ ] New examples:
-    - [ ] All sample runs at the top of the example source file work.
-    - [ ] Update `examples/makefile`:
+   - [ ] All sample runs at the top of the example source file work.
+   - [ ] Update `examples/makefile`:
      - [ ] Add the example code to the appropriate `SEQ_EXAMPLES` and `PAR_EXAMPLES` variables.
      - [ ] Add any files generated by it to the `clean` target.
      - [ ] Add the example binary and any files generated by it to the top-level `.gitignore` file.
-    - [ ] Update `examples/CMakeLists.txt`:
+   - [ ] Update `examples/CMakeLists.txt`:
      - [ ] Add the example code to the `ALL_EXE_SRCS` variable.
      - [ ] Make sure `THIS_TEST_OPTIONS` is set correctly for the new example.
   - [ ] List the new example in `doc/CodeDocumentation.dox`.
-   - [ ] If new examples directory (e.g.`examples/pumi`), list it in `doc/CodeDocumentation.conf.in`
+   - [ ] If new examples directory (e.g. `examples/pumi`), list it in `doc/CodeDocumentation.conf.in`
   - [ ] Companion pull request for documentation in [mfem/web](https://github.com/mfem/web) repo:
      - [ ] Update or add example-specific documentation, see e.g. the `src/examples.md`.
      - [ ] Add the description, labels and screenshots in `src/examples.md` and `src/img`.
@@ -575,13 +573,13 @@ Before a PR can be merged, it should satisfy the following:
      - [ ] Add/update the `CMakeLists.txt` file in the new miniapp directory.
      - [ ] Consider adding a new test for the new miniapp.
   - [ ] List the new miniapp in `doc/CodeDocumentation.dox`
-   - [ ] If new miniapps directory (e.g.`miniapps/nurbs`), add it to `MINIAPP_SUBDIRS` in the `makefile`.
-   - [ ] If new miniapps directory (e.g.`miniapps/nurbs`), list it in `doc/CodeDocumentation.conf.in`
+   - [ ] If new miniapps directory (e.g. `miniapps/nurbs`), add it to `MINIAPP_SUBDIRS` in the `makefile`.
+   - [ ] If new miniapps directory (e.g. `miniapps/nurbs`), list it in `doc/CodeDocumentation.conf.in`
   - [ ] Companion pull request for documentation in [mfem/web](https://github.com/mfem/web) repo:
-     - [ ] Update or add miniapp-specific documentation, see e.g. the `src/meshing.md` and `src/electromagnetics.md` files.
-     - [ ] Add the description, labels and screenshots in `src/examples.md` and `src/img`.
-     - [ ] The miniapps go at the end of the page, and are usually listed only under a specific "Application (PDE)" category.
-     - [ ] Add a short description of the miniapp in the "Extensive Examples" section of `features.md`.
+      - [ ] Update or add miniapp-specific documentation, see e.g. the `src/meshing.md` and `src/electromagnetics.md` files.
+      - [ ] Add the description, labels and screenshots in `src/examples.md` and `src/img`.
+      - [ ] The miniapps go at the end of the page, and are usually listed only under a specific "Application (PDE)" category.
+      - [ ] Add a short description of the miniapp in the "Extensive Examples" section of `features.md`.
 - [ ] New capability:
   - [ ] All new public, protected, and private classes, methods, data members, and functions have full Doxygen-style documentation in source comments. Documentation should include descriptions of member data, function arguments and return values, template parameters, and prerequisites for calling new functions.
   - [ ] Pointer arguments and return values must specify whether ownership is being transferred or lent with the call.
@@ -683,7 +681,7 @@ MFEM uses a `master`/`next`-branch workflow as described below:
 - [ ] Update URL shortlinks:
    - [ ] Create a shortlink at [http://bit.ly/](http://bit.ly/) for the release tarball, e.g. https://mfem.github.io/releases/mfem-3.1.tgz.
    - [ ] (LLNL only) Add and commit the new shortlink in the `links` and `links-mfem` files of the internal `mfem/downloads` repo.
-    - [ ] Add the new shortlinks to the MFEM packages in `spack`, `homebrew/science`, `VisIt`, etc.
+    - [ ] Add the new shortlinks to the MFEM package in `spack`.
 - [ ] Update website in `mfem/web` repo:
    - Update version and shortlinks in `src/index.md` and `src/download.md`.
    -  Use [cloc-1.62.pl](http://cloc.sourceforge.net/) and `ls -lh` to estimate the SLOC and the tarball size in `src/download.md`.
@@ -735,22 +733,24 @@ commit or push, see the [README](config/githooks/README.md) in the `config/githo
 directory.


-### Linux and Mac smoke tests
+### GitHub Actions smoke tests
+
 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.
+GitHub Actions testing should be kept lightweight, as there is a time
+constraint on jobs. The current workflows cover Linux, macOS, and Windows
+configurations.

 - 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.

+### Additional Windows smoke test

-### Windows smoke test
-We use Appveyor to test building with the MS Visual C++ compiler in a Windows
-environment, as well as to test the CMake build. See the `.appveyor` file and the
-build logs at
+We also use Appveyor to test building with the MS Visual C++ compiler in a Windows
+environment, as well as to test the CMake build. See the `.appveyor.yml` file
+and the build logs at
 [https://ci.appveyor.com/project/mfem/mfem](https://ci.appveyor.com/project/mfem/mfem).

 CMake is used to generate the MSVC Project files and drive the build.  A release
@@ -38,14 +38,13 @@ the option MFEM_USE_METIS.
 MFEM also includes support for devices such as GPUs, and programming models such
 as CUDA, HIP, OCCA, OpenMP and RAJA.

- Starting with version 4.0, MFEM requires a C++11 compiler. We recommend using
-  a newer compiler, e.g. GCC version 4.9 or higher.
+- Starting with version 4.9, MFEM requires a C++17 compiler.

 - CUDA support requires an NVIDIA GPU and an installation of the CUDA Toolkit
  https://developer.nvidia.com/cuda-toolkit

 - HIP support requires an AMD GPU and an installation of the ROCm software stack
-  https://rocmdocs.amd.com
+  https://rocm.docs.amd.com

 - OCCA support requires the OCCA library
  https://libocca.org
@@ -83,9 +82,9 @@ Serial build:
 Parallel build:
   (download hypre and METIS 4 from above URLs)
   (build METIS 4 in ../metis-4.0 relative to mfem/)
+   (for METIS 5, see https://mfem.org/building/#parallel-build-using-metis-5)
   (build hypre in ../hypre relative to mfem/)
   make parallel -j 4
-(For METIS 5, see https://mfem.org/building/#parallel-build-using-metis-5)

 CUDA build:
   make cuda -j 4
@@ -115,14 +114,14 @@ Serial build:
 Parallel build:
   (download hypre and METIS 4 from above URLs)
   (build METIS 4 in ../metis-4.0 relative to mfem/)
+   (for METIS 5, see https://mfem.org/building/#parallel-build-using-metis-5)
   (build hypre in ../hypre relative to mfem/)
   mkdir <mfem-build-dir> ; cd <mfem-build-dir>
   cmake <mfem-source-dir> -DMFEM_USE_MPI=YES
   make -j 4
-(For METIS 5, see https://mfem.org/building/#parallel-build-using-metis-5)

 Parallel build with fetching of hypre and METIS:
-   mkdir <mfem-buil-dir> ; cd <mfem-build-dir>
+   mkdir <mfem-build-dir> ; cd <mfem-build-dir>
   cmake <mfem-source-dir> -DMFEM_USE_MPI=YES -DMFEM_FETCH_TPLS=YES
   make -j 4

@@ -134,7 +133,8 @@ CUDA build:

 HIP build:
   mkdir <mfem-build-dir> ; cd <mfem-build-dir>
-   cmake <mfem-source-dir> -DMFEM_USE_HIP=YES -DHIP_ARCH=gfx942 -DCMAKE_CXX_COMPILER=amdclang++ -DCMAKE_HIP_COMPILER=amdclang++
+   cmake <mfem-source-dir> -DMFEM_USE_HIP=YES -DHIP_ARCH=gfx942 \
+         -DCMAKE_CXX_COMPILER=amdclang++ -DCMAKE_HIP_COMPILER=amdclang++
   make -j 4

 Example codes (serial/parallel, depending on the build):
@@ -269,6 +269,7 @@ Compilers:
   CXX      - C++ compiler, serial build
   MPICXX   - MPI C++ compiler, parallel build
   CUDA_CXX - The CUDA compiler, 'nvcc' or 'clang++'
+   HIP_CXX  - The HIP compiler, e.g. 'hipcc'

 Compiler options:
   OPTIM_FLAGS - Options for optimized build
@@ -395,6 +396,11 @@ MFEM_USE_STRUMPACK = YES/NO
   classes. When enabled, this option uses the STRUMPACK_* library options, see
   below.

+MFEM_USE_CUDSS = YES/NO
+   Enable MFEM functionality based on the cuDSS library. When using cuDSS, CUDA
+   support must be also enabled in MFEM, i.e. MFEM_USE_CUDA=YES must be set.
+   When enabled, this option uses the CUDSS_* library options, see below.
+
 MFEM_USE_GINKGO = YES/NO
   Enable MFEM functionality based on the Ginkgo library, which provides
   iterative linear solvers and preconditioners with OpenMP, CUDA backends, see
@@ -554,13 +560,13 @@ MFEM_USE_RAJA = YES/NO
 MFEM_USE_OCCA = YES/NO
   Enables support for the OCCA library in MFEM. OCCA is an open-source library
   which aims to make it easy to program different types of devices (e.g. CPU,
-   GPU, FPGA) by providing an unified API for interacting with JIT-compiled
+   GPU, FPGA) by providing a unified API for interacting with JIT-compiled
   backends. In order to use the OCCA CUDA backend, CUDA support must be enabled
   in MFEM as well, i.e. MFEM_USE_CUDA=YES must be set.

 MFEM_USE_GSLIB = YES/NO
   Enables MFEM functionality based on the GSLIB library, and specifically its
-   FindPoints component, which provides a robust algorithms to evaluate finite
+   FindPoints component, which provides robust algorithms to evaluate finite
   element functions in a collection of points in physical space. When enabled,
   the user can use the GSLIB-FindPoints methods as shown in miniapps/gslib.

@@ -719,9 +725,18 @@ The specific libraries and their options are:
  Options: STRUMPACK_OPT, STRUMPACK_LIB.
  Versions: STRUMPACK >= 3.0.0.

- Ginkgo (optional), used when MFEM_USE_GINKGO = YES. Note that Ginkgo needs a
-  C++ compiler that supports the C++-17 standard. For additional requirements
-  and dependencies of specific modules, see the Ginkgo webpage below.
+- CUDSS (optional), used when MFEM_USE_CUDSS = YES. Note that CUDSS requires 
+  CUDA 12.x toolkit and the cuDSS libraries. The supported communication backend
+  is OpenMPI 4.x (default), and OpenMPI 4.x or a later version must be pre-built.
+  The source files in the cuDSS tarball provide guidance for developing custom 
+  MPI implementations.
+  URL: https://developer.nvidia.com/cudss
+       https://docs.nvidia.com/cuda/cudss/advanced_features.html#communication-layer-library-in-cudss
+  Options: CUDSS_OPT, CUDSS_LIB.
+  Versions: cuDSS >= 0.6.0.
+
+- Ginkgo (optional), used when MFEM_USE_GINKGO = YES. Ginkgo may have additional
+  requirements and module-specific dependencies; see the webpage below.
  URL: https://ginkgo-project.github.io
  Options: GINKGO_OPT, GINKGO_LIB, GINKGO_DIR, GINKGO_BUILD_TYPE (Release or
           Debug).
@@ -793,7 +808,7 @@ The specific libraries and their options are:
  Options: CONDUIT_OPT, CONDUIT_LIB.
  Versions: Conduit >= 0.3.1.

- ADIOS2 (optional) used when MFEM_USE_ADIOS2 = YES.
+- ADIOS2 (optional), used when MFEM_USE_ADIOS2 = YES.
  URL: https://adios2.readthedocs.io/
  Versions: ADIOS >= 2.5.0.

@@ -869,7 +884,7 @@ The specific libraries and their options are:
  Options: RAJA_DIR, RAJA_OPT, RAJA_LIB.
  Versions: RAJA >= 2022.10.3.

- Moonolith (optional), use when MFEM_USE_MOONOLITH = YES.
+- Moonolith (optional), used when MFEM_USE_MOONOLITH = YES.
  URL: https://bitbucket.org/zulianp/par_moonolith
  Options: MOONOLITH_DIR
  Versions: MOONOLITH >= 1.1.0.
@@ -957,7 +972,7 @@ CMAKE_BUILD_TYPE which can be set to standard values like "Debug", and "Release"
 To use a specific generator use the "-G <generator>" option of cmake:

   cmake <mfem-source-dir> -G "Xcode"
-   cmake <mfem-source-dir> -G "Visual Studio 12 2013"
+   cmake <mfem-source-dir> -G "Visual Studio 17 2022"
   cmake <mfem-source-dir> -G "MinGW Makefiles"

 With CMake it is possible to build MFEM as a shared library using the standard
@@ -1202,7 +1217,7 @@ larger problems, there are two options:
 Specific options for HIP
 ========================
 MFEM expects the `ROCM_PATH` environment variable to be set to the path of the
-ROCM install, as well as having `$ROCM_PATH/bin` in `PATH`.
+ROCm install, as well as having `$ROCM_PATH/bin` in `PATH`.

 Specific options for RAJA+HIP+MPI
 =================================
@@ -35,6 +35,7 @@ set(MFEM_USE_SUITESPARSE @MFEM_USE_SUITESPARSE@)
 set(MFEM_USE_SUPERLU @MFEM_USE_SUPERLU@)
 set(MFEM_USE_MUMPS @MFEM_USE_MUMPS@)
 set(MFEM_USE_STRUMPACK @MFEM_USE_STRUMPACK@)
+set(MFEM_USE_CUDSS @MFEM_USE_CUDSS@)
 set(MFEM_USE_GINKGO @MFEM_USE_GINKGO@)
 set(MFEM_USE_AMGX @MFEM_USE_AMGX@)
 set(MFEM_USE_MAGMA @MFEM_USE_MAGMA@)
@@ -109,6 +110,10 @@ if (MFEM_USE_RAJA)
   find_dependency(RAJA)
 endif()

+if (MFEM_USE_CUDSS)
+  find_dependency(cudss)
+endif (MFEM_USE_CUDSS)
+
 if (MFEM_USE_UMPIRE)
   find_dependency(umpire)
 endif()
@@ -108,6 +108,15 @@
 // Enable MFEM functionality based on the STRUMPACK library.
 #cmakedefine MFEM_USE_STRUMPACK

+// Enable MFEM functionality based on the cuDSS library.
+#cmakedefine MFEM_USE_CUDSS
+
+// CUDSS communication layer library path
+#cmakedefine MFEM_CUDSS_COMM_LIB "@MFEM_CUDSS_COMM_LIB@"
+
+// CUDSS threading layer library path
+#cmakedefine MFEM_CUDSS_THREADING_LIB "@MFEM_CUDSS_THREADING_LIB@"
+
 // Enable functionality based on the Ginkgo library.
 #cmakedefine MFEM_USE_GINKGO

@@ -0,0 +1,68 @@
+if (NOT cudss_DIR AND CUDSS_DIR)
+  set(cudss_DIR ${CUDSS_DIR}/lib/cmake/cudss)
+endif()
+message(STATUS "Looking for CUDSS ...")
+message(STATUS "   in CUDSS_DIR = ${CUDSS_DIR}")
+message(STATUS "      cudss_DIR = ${cudss_DIR}")
+find_package(cudss)
+set(CUDSS_FOUND ${cudss_FOUND})
+set(CUDSS_LIBRARIES "cudss")
+if (CUDSS_FOUND)
+  message(STATUS
+    "Found CUDSS target: ${CUDSS_LIBRARIES} (version: ${cudss_VERSION})")
+else()
+  set(msg STATUS)
+  if (CUDSS_FIND_REQUIRED)
+    set(msg FATAL_ERROR)
+  endif()
+  message(${msg}
+    "CUDSS not found. Please set CUDSS_DIR to the install prefix.")
+endif()
+
+if(CUDSS_FOUND AND TARGET cudss)
+  get_target_property(CUDSS_LIBRARY_LOCATION cudss IMPORTED_LOCATION)
+  if(NOT CUDSS_LIBRARY_LOCATION)
+    get_target_property(CUDSS_LIBRARY_LOCATION cudss IMPORTED_LOCATION_RELEASE)
+  endif()
+  if(CUDSS_LIBRARY_LOCATION)
+    get_filename_component(CUDSS_LIBRARY_DIR "${CUDSS_LIBRARY_LOCATION}" DIRECTORY)
+  else()
+    message(WARNING "Could not determine the location of the cuDSS library.")
+  endif()
+else()
+    message(WARNING "cuDSS target not available; cannot determine library directory.")
+endif()
+
+# Set the full name of the cuDSS threading library if OpenMP is enabled.
+# The threading layer library (libcudss_mtlayer_gomp.so) is located under the
+# cuDSS library directory by default.
+if (MFEM_USE_OPENMP)
+  find_file(
+    CUDSS_THREADING_LIB
+    NAMES libcudss_mtlayer_gomp.so
+    PATHS ${CUDSS_LIBRARY_DIR}
+    NO_DEFAULT_PATH
+  )
+  if (NOT DEFINED MFEM_CUDSS_THREADING_LIB AND CUDSS_THREADING_LIB)
+    set(MFEM_CUDSS_THREADING_LIB "${CUDSS_THREADING_LIB}")
+  endif()
+  message(STATUS "CUDSS threading layer library: ${MFEM_CUDSS_THREADING_LIB}")
+endif()
+
+# Set the full name of the cuDSS communication library if MFEM use OpenMPI.
+# The communication layer library (libcudss_commlayer_mpi.so) is located under the
+# cuDSS library directory by default.
+# The communication layer library is used pre-built communication layers for OpenMPI 
+# by default.
+if (MFEM_USE_MPI)
+  find_file(
+    CUDSS_COMM_LIB
+    NAMES libcudss_commlayer_openmpi.so
+    PATHS ${CUDSS_LIBRARY_DIR}
+    NO_DEFAULT_PATH
+  )
+  if (NOT DEFINED MFEM_CUDSS_COMM_LIB AND CUDSS_COMM_LIB)
+    set(MFEM_CUDSS_COMM_LIB "${CUDSS_COMM_LIB}")
+  endif()
+  message(STATUS "CUDSS communication layer library: ${MFEM_CUDSS_COMM_LIB}")
+endif()
@@ -157,4 +157,10 @@ constexpr real_t operator""_r(unsigned long long v)
 #endif
 #endif // MFEM_USE_MPI not defined

+#ifndef MFEM_USE_CUDA
+#ifdef MFEM_USE_CUDSS
+#error Building with cuDSS (MFEM_USE_CUDSS=YES) requires CUDA (MFEM_USE_CUDA=YES)
+#endif
+#endif // MFEM_USE_CUDSS not defined
+
 #endif // MFEM_CONFIG_HPP
@@ -108,6 +108,15 @@
 // Enable MFEM functionality based on the STRUMPACK library.
 // #define MFEM_USE_STRUMPACK

+// Enable MFEM functionality based on the cuDSS library.
+// #define MFEM_USE_CUDSS
+
+// CUDSS communication layer library path
+// #define MFEM_CUDSS_COMM_LIB "@MFEM_CUDSS_COMM_LIB@"
+
+// CUDSS threading layer library path
+// #define MFEM_CUDSS_THREADING_LIB "@MFEM_CUDSS_THREADING_LIB@"
+
 // Enable MFEM features based on the Ginkgo library.
 // #define MFEM_USE_GINKGO

@@ -36,6 +36,9 @@ MFEM_USE_SUPERLU       = @MFEM_USE_SUPERLU@
 MFEM_USE_SUPERLU5      = @MFEM_USE_SUPERLU5@
 MFEM_USE_MUMPS         = @MFEM_USE_MUMPS@
 MFEM_USE_STRUMPACK     = @MFEM_USE_STRUMPACK@
+MFEM_USE_CUDSS         = @MFEM_USE_CUDSS@
+MFEM_CUDSS_COMM_LIB    = @MFEM_CUDSS_COMM_LIB@
+MFEM_CUDSS_THREADING_LIB = @MFEM_CUDSS_THREADING_LIB@
 MFEM_USE_GINKGO        = @MFEM_USE_GINKGO@
 MFEM_USE_AMGX          = @MFEM_USE_AMGX@
 MFEM_USE_MAGMA         = @MFEM_USE_MAGMA@
@@ -38,6 +38,7 @@ option(MFEM_USE_SUPERLU "Enable SuperLU_DIST usage" OFF)
 option(MFEM_USE_SUPERLU5 "Use the old SuperLU_DIST 5.1 version" OFF)
 option(MFEM_USE_MUMPS "Enable MUMPS usage" OFF)
 option(MFEM_USE_STRUMPACK "Enable STRUMPACK usage" OFF)
+option(MFEM_USE_CUDSS "Enable cuDSS usage" OFF)
 option(MFEM_USE_GINKGO "Enable Ginkgo usage" OFF)
 option(MFEM_USE_AMGX "Enable AmgX usage" OFF)
 option(MFEM_USE_MAGMA "Enable MAGMA usage" OFF)
@@ -153,6 +153,7 @@ MFEM_USE_SUPERLU       = NO
 MFEM_USE_SUPERLU5      = NO
 MFEM_USE_MUMPS         = NO
 MFEM_USE_STRUMPACK     = NO
+MFEM_USE_CUDSS         = NO
 MFEM_USE_GINKGO        = NO
 MFEM_USE_AMGX          = NO
 MFEM_USE_MAGMA         = NO
@@ -368,6 +369,19 @@ STRUMPACK_OPT = -I$(STRUMPACK_DIR)/include $(SCOTCH_OPT)
 STRUMPACK_LIB = -L$(STRUMPACK_DIR)/lib -lstrumpack $(MPI_FORTRAN_LIB)\
 $(SCOTCH_LIB) $(SCALAPACK_LIB)

+# CUDSS library configuration
+CUDSS_DIR         = @MFEM_DIR@/../cudss
+CUDSS_INCLUDE_DIR = $(CUDSS_DIR)/include
+CUDSS_LIBRARY_DIR = $(CUDSS_DIR)/lib
+CUDSS_OPT         = -I$(CUDSS_INCLUDE_DIR)
+CUDSS_LIB         = \
+ $(XLINKER)-rpath,$(CUDSS_LIBRARY_DIR) -L$(CUDSS_LIBRARY_DIR) -lcudss
+# The cuDSS communication and threading libraries. 
+MFEM_CUDSS_COMM_LIB = $(abspath $(wildcard $(or $(CUDSS_COMM_LIB),\
+   $(subst @MFEM_DIR@,$(MFEM_DIR), $(CUDSS_LIBRARY_DIR)/libcudss_commlayer_openmpi.so))))
+MFEM_CUDSS_THREADING_LIB = $(abspath $(wildcard $(or $(CUDSS_THREADING_LIB),\
+   $(subst @MFEM_DIR@,$(MFEM_DIR),$(CUDSS_LIBRARY_DIR)/libcudss_mtlayer_gomp.so))))
+
 # Ginkgo library configuration
 GINKGO_DIR = @MFEM_DIR@/../ginkgo/install
 GINKGO_SEARCH_DIR = $(subst @MFEM_DIR@,$(MFEM_DIR),$(GINKGO_DIR))
@@ -621,7 +635,7 @@ PARELAG_LIB = -L$(PARELAG_DIR)/build/src -lParELAG
 AXOM_DIR = @MFEM_DIR@/../axom
 TRIBOL_DIR = @MFEM_DIR@/../tribol
 TRIBOL_OPT = -I$(TRIBOL_DIR)/include -I$(AXOM_DIR)/include
-TRIBOL_LIB = -L$(TRIBOL_DIR)/lib -ltribol -lredecomp -L$(AXOM_DIR)/lib -laxom_mint\
+TRIBOL_LIB = -L$(TRIBOL_DIR)/lib -ltribol -ltribol_shared -lredecomp -L$(AXOM_DIR)/lib -laxom_mint\
   -laxom_slam -laxom_slic -laxom_core

 # Enzyme configuration
@@ -224,17 +224,29 @@ int main(int argc, char *argv[])
   // 11. Solve the linear system A X = B.
   if (!pa)
   {
-#ifndef MFEM_USE_SUITESPARSE
-      // Use a simple symmetric Gauss-Seidel preconditioner with PCG.
-      GSSmoother M((SparseMatrix&)(*A));
-      PCG(*A, M, B, X, 1, 200, 1e-12, 0.0);
-#else
-      // If MFEM was compiled with SuiteSparse, use UMFPACK to solve the system.
-      UMFPackSolver umf_solver;
-      umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
-      umf_solver.SetOperator(*A);
-      umf_solver.Mult(B, X);
+#ifdef MFEM_USE_CUDSS
+      if (Device::Allows(Backend::CUDA_MASK))
+      {
+         // Use cuDSS to solve the system.
+         CuDSSSolver cudss_solver;
+         cudss_solver.SetOperator(*A);
+         cudss_solver.Mult(B, X);
+      }
+      else
 #endif
+      {
+#ifndef MFEM_USE_SUITESPARSE
+         // Use a simple symmetric Gauss-Seidel preconditioner with PCG.
+         GSSmoother M((SparseMatrix&)(*A));
+         PCG(*A, M, B, X, 1, 200, 1e-12, 0.0);
+#else
+         // If MFEM was compiled with SuiteSparse, use UMFPACK to solve the system.
+         UMFPackSolver umf_solver;
+         umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
+         umf_solver.SetOperator(*A);
+         umf_solver.Mult(B, X);
+#endif
+      }
   }
   else
   {
@@ -273,7 +285,7 @@ int main(int argc, char *argv[])
   if (visualization)
   {
      char vishost[] = "localhost";
-      int  visport   = 19916;
+      int visport = 19916;
      socketstream sol_sock(vishost, visport);
      sol_sock.precision(8);
      sol_sock << "solution\n" << mesh << x << flush;
@@ -83,6 +83,9 @@ int main(int argc, char *argv[])
   const char *device_config = "cpu";
   bool visualization = true;
   bool algebraic_ceed = false;
+#ifdef MFEM_USE_CUDSS
+   bool cudss_solver = false;
+#endif

   OptionsParser args(argc, argv);
   args.AddOption(&mesh_file, "-m", "--mesh",
@@ -102,6 +105,10 @@ int main(int argc, char *argv[])
   args.AddOption(&algebraic_ceed, "-a", "--algebraic",
                  "-no-a", "--no-algebraic",
                  "Use algebraic Ceed solver");
+#endif
+#ifdef MFEM_USE_CUDSS
+   args.AddOption(&cudss_solver, "-cudss", "--cudss-solver", "-no-cudss",
+                  "--no-cudss-solver", "Use the cuDSS Solver.");
 #endif
   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                  "--no-visualization",
@@ -248,33 +255,51 @@ int main(int argc, char *argv[])
   // 13. Solve the linear system A X = B.
   //     * With full assembly, use the BoomerAMG preconditioner from hypre.
   //     * With partial assembly, use Jacobi smoothing, for now.
-   Solver *prec = NULL;
-   if (pa)
+#ifdef MFEM_USE_CUDSS
+   if (!pa && (Device::Allows(Backend::CUDA_MASK) && cudss_solver))
   {
-      if (UsesTensorBasis(fespace))
-      {
-         if (algebraic_ceed)
-         {
-            prec = new ceed::AlgebraicSolver(a, ess_tdof_list);
-         }
-         else
-         {
-            prec = new OperatorJacobiSmoother(a, ess_tdof_list);
-         }
-      }
+      // Solve using a direct solver with cuDSS
+      CuDSSSolver cudss_solver(MPI_COMM_WORLD);
+      cudss_solver.SetMatrixSymType(
+         CuDSSSolver::SYMMETRIC_POSITIVE_DEFINITE);
+      cudss_solver.SetMatrixViewType(CuDSSSolver::UPPER);
+      cudss_solver.SetOperator(*A);
+      cudss_solver.Mult(B, X);
   }
   else
+#endif
   {
-      prec = new HypreBoomerAMG;
+      Solver *prec = NULL;
+      if (pa)
+      {
+         if (UsesTensorBasis(fespace))
+         {
+            if (algebraic_ceed)
+            {
+               prec = new ceed::AlgebraicSolver(a, ess_tdof_list);
+            }
+            else
+            {
+               prec = new OperatorJacobiSmoother(a, ess_tdof_list);
+            }
+         }
+      }
+      else
+      {
+         prec = new HypreBoomerAMG;
+      }
+      CGSolver cg(MPI_COMM_WORLD);
+      cg.SetRelTol(1e-12);
+      cg.SetMaxIter(2000);
+      cg.SetPrintLevel(1);
+      if (prec)
+      {
+         cg.SetPreconditioner(*prec);
+      }
+      cg.SetOperator(*A);
+      cg.Mult(B, X);
+      delete prec;
   }
-   CGSolver cg(MPI_COMM_WORLD);
-   cg.SetRelTol(1e-12);
-   cg.SetMaxIter(2000);
-   cg.SetPrintLevel(1);
-   if (prec) { cg.SetPreconditioner(*prec); }
-   cg.SetOperator(*A);
-   cg.Mult(B, X);
-   delete prec;

   // 14. Recover the parallel grid function corresponding to X. This is the
   //     local finite element solution on each processor.
@@ -830,15 +830,9 @@ ParComplexGridFunction::ParComplexGridFunction(ParMesh *m, std::istream &input)
      int vsize = pfes->GetVSize();
      Vector::Load(input, 2*vsize);

-      real_t *data_  = const_cast<real_t*>(HostRead());
-      for (int i = 0; i < vsize; i++)
-      {
-         if (pfes->GetDofSign(i) < 0)
-         {
-            data_[i] = -data_[i];
-            data_[i+vsize] = -data_[i+vsize];
-         }
-      }
+      real_t *h_data = HostReadWrite();
+      pfes->ApplyDofSigns(h_data);
+      pfes->ApplyDofSigns(h_data + vsize);


      // if the mesh is a legacy (v1.1) NC mesh, it has old vertex ordering
@@ -1051,15 +1045,14 @@ void ParComplexGridFunction::Save(std::ostream &os) const
   os << '\n';

   int vsize = pfes->GetVSize();
-   real_t *data_  = const_cast<real_t*>(HostRead());
-   for (int i = 0; i < vsize; i++)
-   {
-      if (pfes->GetDofSign(i) < 0)
-      {
-         data_[i] = -data_[i];
-         data_[i+vsize] = -data_[i+vsize];
-      }
-   }
+   // We use const_cast + HostRead (instead of HostReadWrite) because we only
+   // need to change the host data temporarily and this way we do not invalidate
+   // the data if it is on device. If we use HostReadWrite here, later calls to
+   // Read or ReadWrite will need to copy the data from host to device. With the
+   // approach used here, the host-to-device copy is avoided.
+   real_t *h_data = const_cast<real_t*>(HostRead());
+   pfes->ApplyDofSigns(h_data);
+   pfes->ApplyDofSigns(h_data + vsize);

   if (pfes->GetOrdering() == Ordering::byNODES)
   {
@@ -1070,14 +1063,8 @@ void ParComplexGridFunction::Save(std::ostream &os) const
      Vector::Print(os, pfes->GetVDim());
   }

-   for (int i = 0; i < vsize; i++)
-   {
-      if (pfes->GetDofSign(i) < 0)
-      {
-         data_[i] = -data_[i];
-         data_[i+vsize] = -data_[i+vsize];
-      }
-   }
+   pfes->ApplyDofSigns(h_data);
+   pfes->ApplyDofSigns(h_data + vsize);

   os.flush();
 }
@@ -114,6 +114,10 @@ void ConduitDataCollection::Save()
                                         n_mesh["fields"][name]);
   }

+   // TODO: in parallel, we need to call ParFiniteElementSpace::ApplyDofSigns
+   // for all ParGridFunction objects before and after saving, see
+   // ParGridFunction::Save.
+
   // save mesh data
   SaveMeshAndFields(myid,
                     n_mesh,
@@ -3363,18 +3363,36 @@ void FindPointsGSLIB::Interpolate(const GridFunction &field_in,
      int gf_order_h1 = std::max(gf_order, 1); // H1 should be at least order 1
      H1_FECollection fec(gf_order_h1, dim);
      const int ncomp = field_in.FESpace()->GetVDim();
-      FiniteElementSpace fes(mesh, &fec, ncomp,
-                             field_in.FESpace()->GetOrdering());
-      GridFunction field_in_h1(&fes);
-      field_in_h1.UseDevice(false);
+
+      std::unique_ptr<FiniteElementSpace> fes;
+      std::unique_ptr<GridFunction> field_in_h1;
+#ifdef MFEM_USE_MPI
+      if (auto *pmesh = dynamic_cast<ParMesh*>(mesh))
+      {
+         fes = std::make_unique<ParFiniteElementSpace>(
+                  pmesh, &fec, ncomp, field_in.FESpace()->GetOrdering());
+
+         field_in_h1 = std::make_unique<ParGridFunction>(
+                          static_cast<ParFiniteElementSpace*>(fes.get()));
+      }
+      else
+#endif
+      {
+         fes = std::make_unique<FiniteElementSpace>(
+                  mesh, &fec, ncomp, field_in.FESpace()->GetOrdering());
+
+         field_in_h1 = std::make_unique<GridFunction>(fes.get());
+      }
+
+      field_in_h1->UseDevice(false);

      if (avgtype == AvgType::ARITHMETIC)
      {
-         field_in_h1.ProjectDiscCoefficient(field_in_dg, GridFunction::ARITHMETIC);
+         field_in_h1->ProjectDiscCoefficient(field_in_dg, GridFunction::ARITHMETIC);
      }
      else if (avgtype == AvgType::HARMONIC)
      {
-         field_in_h1.ProjectDiscCoefficient(field_in_dg, GridFunction::HARMONIC);
+         field_in_h1->ProjectDiscCoefficient(field_in_dg, GridFunction::HARMONIC);
      }
      else
      {
@@ -3383,11 +3401,11 @@ void FindPointsGSLIB::Interpolate(const GridFunction &field_in,

      if (gf_order_h1 == mesh_order) // basis is GaussLobatto by default
      {
-         InterpolateH1(field_in_h1, field_out_l2, field_out_ordering);
+         InterpolateH1(*field_in_h1, field_out_l2, field_out_ordering);
      }
      else
      {
-         InterpolateGeneral(field_in_h1, field_out_l2, field_out_ordering);
+         InterpolateGeneral(*field_in_h1, field_out_l2, field_out_ordering);
      }

      // Copy interpolated values for the points on element border
@@ -284,12 +284,12 @@ GeometricMultigrid::GeometricMultigrid(
   ownedProlongations.SetSize(nlevels - 1);
   ownedProlongations = have_ess_bdr;

-   if (have_ess_bdr)
+   essentialTrueDofs.SetSize(nlevels);
+   for (int level = 0; level < nlevels; ++level)
   {
-      essentialTrueDofs.SetSize(nlevels);
-      for (int level = 0; level < nlevels; ++level)
+      essentialTrueDofs[level] = new Array<int>;
+      if (have_ess_bdr)
      {
-         essentialTrueDofs[level] = new Array<int>;
         fespaces.GetFESpaceAtLevel(level).GetEssentialTrueDofs(
            ess_bdr, *essentialTrueDofs[level]);
      }
@@ -187,8 +187,7 @@ public:
   /// mesh boundary element attributes that define the essential DOFs.
   ///
   /// If @a ess_bdr is empty, or all its entries are 0, then no essential
-   /// boundary conditions are imposed and the protected array essentialTrueDofs
-   /// remains empty.
+   /// boundary conditions are imposed.
   GeometricMultigrid(const FiniteElementSpaceHierarchy& fespaces_,
                      const Array<int> &ess_bdr);

@@ -349,6 +349,7 @@ void ParFiniteElementSpace::GetGroupComm(
      }
   }

+   bool have_sign_flips = false;
   if (g_ldof_sign)
   {
      g_ldof_sign->SetSize(GetNDofs());
@@ -428,6 +429,7 @@ void ParFiniteElementSpace::GetGroupComm(
                  if (g_ldof_sign)
                  {
                     (*g_ldof_sign)[dofs[l]] = -1;
+                     have_sign_flips = true;
                  }
               }
               else
@@ -466,6 +468,7 @@ void ParFiniteElementSpace::GetGroupComm(
                  if (g_ldof_sign)
                  {
                     (*g_ldof_sign)[dofs[l]] = -1;
+                     have_sign_flips = true;
                  }
               }
               else
@@ -504,6 +507,7 @@ void ParFiniteElementSpace::GetGroupComm(
                  if (g_ldof_sign)
                  {
                     (*g_ldof_sign)[dofs[l]] = -1;
+                     have_sign_flips = true;
                  }
               }
               else
@@ -527,12 +531,18 @@ void ParFiniteElementSpace::GetGroupComm(
      group_ldof.GetI()[gr+1] = group_ldof_counter;
   }

+   if (g_ldof_sign && have_sign_flips == false)
+   {
+      g_ldof_sign->DeleteAll();
+   }
+
   gc.Finalize();
 }

 void ParFiniteElementSpace::ApplyLDofSigns(Array<int> &dofs) const
 {
   MFEM_ASSERT(Conforming(), "wrong code path");
+   if (!HaveDofSigns()) { return; }

   for (int i = 0; i < dofs.Size(); i++)
   {
@@ -559,6 +569,24 @@ void ParFiniteElementSpace::ApplyLDofSigns(Table &el_dof) const
   ApplyLDofSigns(all_dofs);
 }

+void ParFiniteElementSpace::ApplyDofSigns(real_t *h_data) const
+{
+   if (!HaveDofSigns()) { return; }
+
+   const bool byvdim = (ordering == Ordering::byVDIM);
+   for (int i = 0; i < ndofs; i++)
+   {
+      if (ldof_sign[i] < 0)
+      {
+         for (int d = 0; d < vdim; d++)
+         {
+            const int idx = byvdim ? d+vdim*i : i+ndofs*d;
+            h_data[idx] = -h_data[idx];
+         }
+      }
+   }
+}
+
 void ParFiniteElementSpace::GetElementDofs(int i, Array<int> &dofs,
                                           DofTransformation &doftrans) const
 {
@@ -1193,15 +1221,15 @@ void ParFiniteElementSpace::GetEssentialTrueDofsVar(const Array<int>
   MFEM_VERIFY(IsVariableOrder() && R,
               "GetEssentialTrueDofsVar is only for variable-order spaces");

-   true_ess_dofs.SetSize(R->Height(), Device::GetDeviceMemoryType());
+   true_ess_dofs.SetSize(R->Height());
+   true_ess_dofs.HostWrite();
+   true_ess_dofs = 0;

   const int ntdofs = tdof2ldof.Size();
   MFEM_VERIFY(vdim * ntdofs == R->NumRows() &&
               vdim * ntdofs == true_ess_dofs.Size(), "");
   MFEM_VERIFY(ldof_ltdof.Size() == ndofs && ess_dofs.Size() == vdim * ndofs, "");

-   true_ess_dofs = 0;
-
   const bool bynodes = (ordering == Ordering::byNODES);
   const int vdim_factor = bynodes ? 1 : vdim;
   const int num_true_dofs = R->NumRows() / vdim;
@@ -340,8 +340,20 @@ public:

   inline ParMesh *GetParMesh() const { return pmesh; }

-   int GetDofSign(int i)
-   { return NURBSext || Nonconforming() ? 1 : ldof_sign[VDofToDof(i)]; }
+   /** @brief Return true if the parallel FE space has DOFs with signs opposite
+       of the DOFs in the respective serial FE space. */
+   bool HaveDofSigns() const { return ldof_sign.Size() != 0; }
+
+   /** @brief Apply the DOF signs to the given host data @a h_data which must be
+       of size GetVSize() if HaveDofSigns() is true. If HaveDofSigns() is false,
+       this method is no-op and returns immediately. */
+   void ApplyDofSigns(real_t *h_data) const;
+
+   /** @brief Return -1 if the given (vector) DOF @a i has a sign opposite of
+       the DOF in the respecive serial FE space. Otherwise, return 1. */
+   int GetDofSign(int i) const
+   { return !HaveDofSigns() ? 1 : ldof_sign[VDofToDof(i)]; }
+
   HYPRE_BigInt *GetDofOffsets()     const { return dof_offsets; }
   HYPRE_BigInt *GetTrueDofOffsets() const { return tdof_offsets; }
   HYPRE_BigInt GlobalVSize() const
@@ -80,6 +80,8 @@ ParGridFunction::ParGridFunction(ParMesh *pmesh, std::istream &input)
                                    fes->GetOrdering());
   delete fes;
   fes = pfes;
+
+   pfes->ApplyDofSigns(HostReadWrite());
 }

 void ParGridFunction::Update()
@@ -1082,18 +1084,17 @@ real_t ParGridFunction::ComputeDGFaceJumpError(Coefficient *exsol,

 void ParGridFunction::Save(std::ostream &os) const
 {
-   real_t *data_  = const_cast<real_t*>(HostRead());
-   for (int i = 0; i < size; i++)
-   {
-      if (pfes->GetDofSign(i) < 0) { data_[i] = -data_[i]; }
-   }
+   // We use const_cast + HostRead (instead of HostReadWrite) because we only
+   // need to change the host data temporarily and this way we do not invalidate
+   // the data if it is on device. If we use HostReadWrite here, later calls to
+   // Read or ReadWrite will need to copy the data from host to device. With the
+   // approach used here, the host-to-device copy is avoided.
+   real_t *h_data = const_cast<real_t*>(HostRead());
+   pfes->ApplyDofSigns(h_data);

   GridFunction::Save(os);

-   for (int i = 0; i < size; i++)
-   {
-      if (pfes->GetDofSign(i) < 0) { data_[i] = -data_[i]; }
-   }
+   pfes->ApplyDofSigns(h_data);
 }

 void ParGridFunction::Save(const char *fname, int precision) const
@@ -1264,7 +1265,13 @@ void ParGridFunction::SaveAsOne(std::ostream &os) const
   int *nfdofs = new int[NRanks];
   int *nrdofs = new int[NRanks];

+   // We use const_cast + HostRead (instead of HostReadWrite) because we only
+   // need to change the host data temporarily and this way we do not invalidate
+   // the data if it is on device. If we use HostReadWrite here, later calls to
+   // Read or ReadWrite will need to copy the data from host to device. With the
+   // approach used here, the host-to-device copy is avoided.
   real_t * h_data = const_cast<real_t *>(this->HostRead());
+   pfes->ApplyDofSigns(h_data);  // temporarily flip the dof signs

   values[0] = h_data;
   nv[0]     = pfes -> GetVSize();
@@ -1371,6 +1378,8 @@ void ParGridFunction::SaveAsOne(std::ostream &os) const
      MPI_Send(h_data, nv[0], MPITypeMap<real_t>::mpi_type, 0, 460, MyComm);
   }

+   pfes->ApplyDofSigns(h_data);  // restore the original h_data
+
   delete [] values;
   delete [] nv;
   delete [] nvdofs;
@@ -148,6 +148,12 @@ if (MFEM_USE_MKL_PARDISO)
  list(APPEND HDRS pardiso.hpp)
 endif()

+# cudss solver
+if (MFEM_USE_CUDSS)
+  list(APPEND SRCS cudss.cpp)
+  list(APPEND HDRS cudss.hpp)
+endif()
+
 convert_filenames_to_full_paths(SRCS)
 convert_filenames_to_full_paths(HDRS)

@@ -0,0 +1,409 @@
+// 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.
+
+#include "cudss.hpp"
+#include "../general/communication.hpp"
+#include <string>
+
+#ifdef MFEM_USE_CUDSS
+
+#ifdef MFEM_USE_SINGLE
+#define CUDA_REAL_T CUDA_R_32F
+#else
+#define CUDA_REAL_T CUDA_R_64F
+#endif
+
+// Define a cuDSS error check macro, MFEM_CUDSS_CHECK(x), where x returns/is of
+// type 'cudssStatus_t'. This macro evaluates 'x' and raises an error if the
+// result is not CUDSS_STATUS_SUCCESS.
+#define MFEM_CUDSS_CHECK(x)                                          \
+  do {                                                               \
+    cudssStatus_t mfem_err_internal_var_name = (x);                  \
+    if (mfem_err_internal_var_name != CUDSS_STATUS_SUCCESS) {        \
+      ::mfem::mfem_cudss_error(mfem_err_internal_var_name, #x,       \
+                               _MFEM_FUNC_NAME, __FILE__, __LINE__); \
+    }                                                                \
+  } while (0)
+
+namespace mfem
+{
+// Function used by the macro MFEM_CUDSS_CHECK.
+void mfem_cudss_error(cudssStatus_t status, const char *expr, const char *func,
+                      const char *file, int line)
+{
+   mfem::err << "\n\nCUDSS error: (" << expr << ") failed with error:\n --> "
+             << "CUDSS call ended unsuccessfully"
+             << " [code: " << static_cast<int>(status) << ']'
+             << "\n ... in function: " << func << "\n ... in file: " << file
+             << ':' << line << '\n';
+   mfem_error();
+}
+
+CuDSSSolver::CuDSSSolver() { InitCuDSS(); }
+
+#ifdef MFEM_USE_MPI
+CuDSSSolver::CuDSSSolver(MPI_Comm comm_) : mpi_comm(comm_)
+{
+   InitCuDSS();
+
+   // NOTE: Set the communication layer to NULL so that cuDSS picks it
+   // from the environment variable "CUDSS_COMM_LIB"
+   const char* comm_lib = GetEnv("CUDSS_COMM_LIB");
+#ifdef MFEM_CUDSS_COMM_LIB
+   if (comm_lib == nullptr)
+   {
+      comm_lib = MFEM_CUDSS_COMM_LIB;
+   }
+#endif
+   MFEM_CUDSS_CHECK(cudssSetCommLayer(handle, comm_lib));
+
+   MFEM_CUDSS_CHECK(cudssDataSet(handle, solverData, CUDSS_DATA_COMM,
+                                 &mpi_comm, sizeof(MPI_Comm *)));
+}
+#endif // MFEM_USE_MPI
+
+CuDSSSolver::~CuDSSSolver()
+{
+   // Destroy the system Matrix, RHS vector and solution vector
+   if (Ac)
+   {
+      MFEM_CUDSS_CHECK(cudssMatrixDestroy(*Ac));
+      MFEM_CUDSS_CHECK(cudssMatrixDestroy(xc));
+      MFEM_CUDSS_CHECK(cudssMatrixDestroy(yc));
+   }
+
+   // Destroy the cuDSS handle, solver config and solver data
+   MFEM_CUDSS_CHECK(cudssDataDestroy(handle, solverData));
+   MFEM_CUDSS_CHECK(cudssConfigDestroy(solverConfig));
+
+
+   MFEM_CUDSS_CHECK(cudssDestroy(handle));
+   handle = nullptr;
+
+
+   if (csr_offsets_d != NULL)
+   {
+      CuMemFree(csr_offsets_d);
+   }
+
+   if (csr_columns_d != NULL)
+   {
+      CuMemFree(csr_columns_d);
+   }
+
+   if (csr_values_d != NULL)
+   {
+      CuMemFree(csr_values_d);
+   }
+}
+
+void CuDSSSolver::InitCuDSS()
+{
+   // Create the cuDSS handle
+   MFEM_CUDSS_CHECK(cudssCreate(&handle));
+
+#ifdef MFEM_USE_OPENMP
+   // NOTE: Set the threading layer library name to NULL so that cuDSS picks
+   // it from the environment variable "CUDSS_THREADING_LIB"
+   const char* threading_lib = GetEnv("CUDSS_THREADING_LIB");
+#ifdef MFEM_CUDSS_THREADING_LIB
+   if (threading_lib == nullptr)
+   {
+      threading_lib = MFEM_CUDSS_THREADING_LIB;
+   }
+#endif
+   MFEM_CUDSS_CHECK(cudssSetThreadingLayer(handle, threading_lib));
+#endif  // MFEM_USE_OPENMP
+
+   // Create the solver configuration and data objects
+   MFEM_CUDSS_CHECK(cudssConfigCreate(&solverConfig));
+   MFEM_CUDSS_CHECK(cudssDataCreate(handle, &solverData));
+}
+
+void CuDSSSolver::SetMatrixSymType(MatType mtype_)
+{
+   switch (mtype_)
+   {
+      case MatType::SYMMETRIC_INDEFINITE:
+         mat_type = CUDSS_MTYPE_SYMMETRIC;
+         break;
+      case MatType::SYMMETRIC_POSITIVE_DEFINITE:
+         mat_type = CUDSS_MTYPE_SPD;
+         break;
+      default:
+         mat_type = CUDSS_MTYPE_GENERAL;
+         break;
+   }
+}
+
+void CuDSSSolver::SetMatrixViewType(MatViewType mvtype_)
+{
+   // If the MatType is NONSYMMETRIC, the matrix view type must be FULL.
+   if (mat_type == CUDSS_MTYPE_GENERAL)
+   {
+      mview = CUDSS_MVIEW_FULL;
+      return;
+   }
+
+   // If the matrix is symmetric, the following view type will be optional.
+   switch (mvtype_)
+   {
+      case MatViewType::LOWER:
+         mview = CUDSS_MVIEW_LOWER;
+         break;
+      case MatViewType::UPPER:
+         mview = CUDSS_MVIEW_UPPER;
+         break;
+      default:
+         mview = CUDSS_MVIEW_FULL;
+         break;
+   }
+}
+
+void CuDSSSolver::SetReorderingReuse(bool reuse)
+{
+   MFEM_VERIFY(Ac == nullptr,
+               "Set reordering reuse before setting the operator!");
+   reorder_reuse = reuse;
+}
+
+#ifdef MFEM_USE_MPI
+void CuDSSSolver::SetMatrix(const HypreParMatrix &op)
+{
+   bool cuDSSObjectInitialized = (Ac != nullptr);
+
+   hypre_ParCSRMatrix *parcsr_op = op;
+   op.HypreRead();
+   hypre_CSRMatrix *csr_op = hypre_MergeDiagAndOffd(parcsr_op);
+   op.HypreRead();
+#if MFEM_HYPRE_VERSION >= 21600
+   hypre_CSRMatrixBigJtoJ(csr_op);
+#endif
+
+   // Parameters of the Operator
+   n_loc = height;  // Equal to the csr_op->num_rows
+   n_global = internal::to_int(parcsr_op->global_num_rows);
+   row_start = parcsr_op->first_row_index;
+   row_end = row_start + n_loc - 1;
+   MFEM_VERIFY(!cuDSSObjectInitialized || !reorder_reuse ||
+               (reorder_reuse && (nnz == csr_op->num_nonzeros)),
+               "Inconsistent new matrix pattern!");
+   nnz = csr_op->num_nonzeros;
+
+   SetMatrixCuDSS(csr_op->i, csr_op->j, csr_op->data);
+   hypre_CSRMatrixDestroy(csr_op);
+}
+#endif // MFEM_USE_MPI
+
+void CuDSSSolver::SetMatrix(const SparseMatrix &op)
+{
+   bool cuDSSObjectInitialized = (Ac != nullptr);
+
+   // Parameters of the Operator
+   MFEM_VERIFY(!cuDSSObjectInitialized || !reorder_reuse ||
+               (reorder_reuse && (nnz == op.NumNonZeroElems())),
+               "Inconsistent new matrix pattern!");
+
+   SparseMatrix *A = const_cast<SparseMatrix *>(&op);
+
+   nnz = A->NumNonZeroElems();
+   n_global = height;  // Equal to the height in serial
+   n_loc = height;     // Equal to the height in serial
+
+   int *csr_offsets = const_cast<int *>(A->ReadI());
+   int *csr_columns = const_cast<int *>(A->ReadJ());
+   real_t *csr_values = const_cast<real_t *>(A->ReadData());
+
+   SetMatrixCuDSS(csr_offsets, csr_columns, csr_values);
+}
+
+void CuDSSSolver::SetMatrixCuDSS(int *csr_offsets, int *csr_columns,
+                                 real_t *csr_values)
+{
+   bool cuDSSObjectInitialized = (Ac != nullptr);
+   // Initial the cudssMatrix objects
+   if (!cuDSSObjectInitialized)
+   {
+      // Set the cudssMatrix object of csr operator
+      Ac = std::make_unique<cudssMatrix_t>();
+      // Create empty RHS and solution vectors
+      SetNumRHS(1);
+      // Allocate device memory for csr values
+      CuMemAlloc(&csr_values_d, nnz * sizeof(real_t));
+   }
+
+   CuMemcpyDtoD(csr_values_d, csr_values, nnz * sizeof(real_t));
+
+   // New cuDSS CSR matrix object and analysis or reuse the one from a previous
+   // matrix
+   if (!cuDSSObjectInitialized || !reorder_reuse)
+   {
+      if (reorder_reuse)  // !cuDSSObjectInitialized && reorder_reuse
+      {
+         // NOTE: For CuDSS solver to reuse the reordering (skipping analysis
+         // phase), it needs to access the I and J arrays of the **initial**
+         // matrix. Therefore, we need to copy and keep I and J in device memory.
+         CuMemAlloc(&csr_offsets_d, (n_loc + 1) * sizeof(int));
+         CuMemAlloc(&csr_columns_d, nnz * sizeof(int));
+
+         CuMemcpyDtoD(csr_offsets_d, csr_offsets, (n_loc + 1) * sizeof(int));
+         CuMemcpyDtoD(csr_columns_d, csr_columns, nnz * sizeof(int));
+
+         MFEM_CUDSS_CHECK(
+            cudssMatrixCreateCsr(
+               Ac.get(), n_global, n_global, nnz, csr_offsets_d, NULL,
+               csr_columns_d, csr_values_d, CUDA_R_32I, CUDA_REAL_T, mat_type, mview,
+               CUDSS_BASE_ZERO));
+      }
+      else    // !reorder_reuse
+      {
+         if (cuDSSObjectInitialized)
+         {
+            MFEM_CUDSS_CHECK(cudssMatrixDestroy(*Ac));
+         }
+         MFEM_CUDSS_CHECK(
+            cudssMatrixCreateCsr(
+               Ac.get(), n_global, n_global, nnz, csr_offsets, NULL, csr_columns,
+               csr_values_d, CUDA_R_32I, CUDA_REAL_T, mat_type, mview,
+               CUDSS_BASE_ZERO));
+      }
+#ifdef MFEM_USE_MPI
+      if (Mpi::IsInitialized())
+      {
+         MFEM_CUDSS_CHECK(cudssMatrixSetDistributionRow1d(*Ac, row_start, row_end));
+      }
+#endif
+      // Analysis
+      MFEM_CUDSS_CHECK(cudssExecute(handle, CUDSS_PHASE_ANALYSIS, solverConfig,
+                                    solverData, *Ac, yc, xc));
+   }
+   else    // cuDSSObjectInitialized && reorder_reuse
+   {
+      // NOTE: When reusing analysis result, we only update the Data array,
+      // without changing the I and J arrays.
+      MFEM_CUDSS_CHECK(cudssMatrixSetValues(*Ac, csr_values_d));
+   }
+
+   // Factorization
+   MFEM_CUDSS_CHECK(cudssExecute(handle, CUDSS_PHASE_FACTORIZATION, solverConfig,
+                                 solverData, *Ac, yc, xc));
+}
+
+void CuDSSSolver::SetOperator(const Operator &op)
+{
+   bool cuDSSObjectInitialized = (Ac != nullptr);
+   MFEM_VERIFY(
+      !cuDSSObjectInitialized || (height == op.Height() && width == op.Width()),
+      "Inconsistent new matrix size!");
+   height = op.Height();
+   width = op.Width();
+   if (const SparseMatrix *A = dynamic_cast<const SparseMatrix *>(&op))
+   {
+      SetMatrix(*A);
+   }
+#ifdef MFEM_USE_MPI
+   else if (const HypreParMatrix *A =
+               dynamic_cast<const HypreParMatrix *>(&op))
+   {
+      SetMatrix(*A);
+   }
+#endif // MFEM_USE_MPI
+   else
+   {
+      MFEM_ABORT("Unsupported Operator Type \n");
+   }
+}
+
+void CuDSSSolver::SetNumRHS(int nrhs_) const
+{
+   if (nrhs != nrhs_)
+   {
+      if (nrhs > 0)
+      {
+         // Destroy the previous RHS vector and solution vector
+         MFEM_CUDSS_CHECK(cudssMatrixDestroy(xc));
+         MFEM_CUDSS_CHECK(cudssMatrixDestroy(yc));
+      }
+      // Create empty RHS and solution vectors
+      MFEM_CUDSS_CHECK(cudssMatrixCreateDn(&xc, n_global, nrhs_, n_global, NULL,
+                                           CUDA_REAL_T, CUDSS_LAYOUT_COL_MAJOR));
+
+      MFEM_CUDSS_CHECK(cudssMatrixCreateDn(&yc, n_global, nrhs_, n_global, NULL,
+                                           CUDA_REAL_T, CUDSS_LAYOUT_COL_MAJOR));
+
+#ifdef MFEM_USE_MPI
+      MFEM_CUDSS_CHECK(cudssMatrixSetDistributionRow1d(xc, row_start, row_end));
+      MFEM_CUDSS_CHECK(cudssMatrixSetDistributionRow1d(yc, row_start, row_end));
+#endif // MFEM_USE_MPI
+   }
+   nrhs = nrhs_;
+}
+
+void CuDSSSolver::Mult(const Vector &x, Vector &y) const
+{
+   Array<const Vector *> X(1);
+   Array<Vector *> Y(1);
+   X[0] = &x;
+   Y[0] = &y;
+   ArrayMult(X, Y);
+}
+
+void CuDSSSolver::ArrayMult(const Array<const Vector *> &X,
+                            Array<Vector *> &Y) const
+{
+   SetNumRHS(X.Size());
+
+   Vector RHS, SOL;
+
+   if (nrhs == 1)
+   {
+      RHS.MakeRef(*(const_cast<Vector *>(X[0])), 0, X[0]->Size());
+      SOL.MakeRef(*Y[0], 0, Y[0]->Size());
+   }
+   else
+   {
+      // NOTE: RHS must have **global** num_rows and nrhs columns
+      RHS.SetSize(nrhs * n_global, *X[0]);
+      for (int i = 0; i < nrhs; i++)
+      {
+         Vector s(RHS, i * n_global, n_loc);
+         s = *X[i];
+      }
+
+      // NOTE: SOL must have **global** num_rows and nrhs columns
+      SOL.SetSize(nrhs * n_global, *Y[0]);
+   }
+
+   MFEM_CUDSS_CHECK(cudssMatrixSetValues(xc, const_cast<real_t *>(RHS.Read())));
+   MFEM_CUDSS_CHECK(cudssMatrixSetValues(yc, SOL.Write()));
+
+   // Solve
+   MFEM_CUDSS_CHECK(cudssExecute(handle, CUDSS_PHASE_SOLVE, solverConfig,
+                                 solverData, *Ac, yc, xc));
+
+   if (nrhs == 1)
+   {
+      SOL.SyncAliasMemory(*Y[0]);
+   }
+
+   if (nrhs > 1)
+   {
+      // Get solution for each right-hand side
+      for (int i = 0; i < nrhs; i++)
+      {
+         Vector s(SOL, i * n_global, n_loc);
+         *Y[i] = s;
+      }
+   }
+}
+
+} // namespace mfem
+#endif // MFEM_USE_CUDSS
@@ -0,0 +1,224 @@
+// 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.
+
+#ifndef MFEM_CUDSS
+#define MFEM_CUDSS
+
+#include "../config/config.hpp"
+
+#ifdef MFEM_USE_CUDSS
+
+#include "cudss.h"
+#include <memory>
+
+#ifdef MFEM_USE_MPI
+#include <mpi.h>
+#include "hypre.hpp"
+#else
+#include "operator.hpp"
+#include "sparsemat.hpp"
+#endif
+
+namespace mfem
+{
+/**
+ * @brief cuDSS: A high-performance CUDA Library for Direct Sparse Solvers
+ *
+ * Interface for the distributed cuDSS solver
+ */
+class CuDSSSolver : public Solver
+{
+public:
+   /// Specify the type of matrix we are applying the solver to
+   enum MatType
+   {
+      /// CUDSS_MTYPE_GENERAL: General matrix [default].
+      NONSYMMETRIC = 0,
+      /// CUDSS_MTYPE_SYMMETRIC: Real symmetric matrix.
+      SYMMETRIC_INDEFINITE = 1,
+      /// CUDSS_MTYPE_SPD: Symmetric positive-definite matrix.
+      SYMMETRIC_POSITIVE_DEFINITE = 2,
+   };
+
+   /// Specify the view type of matrix we are applying the solver to
+   enum MatViewType
+   {
+      /// CUDSS_MVIEW_FULL: Full matrix [default]
+      FULL = 0,
+      /// CUDSS_MVIEW_LOWER: Lower-triangular matrix (including the diagonal).
+      LOWER = 1,
+      /// CUDSS_MVIEW_UPPER: Upper-triangular matrix (including the diagonal).
+      UPPER = 2,
+   };
+
+   /**
+    * @brief Constructor.
+    */
+   CuDSSSolver();
+
+#ifdef MFEM_USE_MPI
+   /**
+    * @brief Constructor with MPI_Comm parameter.
+    */
+   CuDSSSolver(MPI_Comm comm);
+#endif
+
+   // Note: CuDSSSolver disables the move copy constructor and move assignment
+   // operator
+   CuDSSSolver(CuDSSSolver &&) = delete;
+   CuDSSSolver &operator=(CuDSSSolver &&) = delete;
+
+   /**
+    * @brief Set the matrix type
+    *
+    * Supported matrix types:
+    *          CuDSSSolver::NONSYMMETRIC,
+    *          CuDSSSolver::SYMMETRIC_INDEFINITE,
+    *      and CuDSSSolver::SYMMETRIC_POSITIVE_DEFINITE
+    *
+    * @param mtype_ Matrix type
+    *
+    * @note This method has to be called before SetOperator
+    */
+   void SetMatrixSymType(MatType mtype_);
+
+   /**
+    * @brief Set the matrix view type
+    *
+    * Supported matrix types:
+    *          CuDSSSolver::FULL,
+    *          CuDSSSolver::LOWER,
+    *      and CuDSSSolver::UPPER
+    *
+    * @param mvtype Matrix view type
+    *
+    * @note This method has to be called before SetOperator
+    */
+   void SetMatrixViewType(MatViewType mvtype);
+
+   /**
+    * @brief Set the flag controlling reuse of the symbolic factorization
+    * for multiple operators
+    *
+    * @param reuse Flag to reuse symbolic factorization
+    *
+    * @note This method has to be called before repeated calls to SetOperator
+    */
+   void SetReorderingReuse(bool reuse);
+
+   void SetOperator(const Operator &op) override;
+
+   /**
+    * @brief Solve $ y = Op^{-1} x $
+    *
+    * @param x RHS vector
+    * @param y Solution vector
+    */
+   void Mult(const Vector &x, Vector &y) const override;
+
+   /**
+    * @brief Solve $ Y_i = Op^{-1} X_i $
+    *
+    * @param X Array of RHS vectors
+    * @param Y Array of Solution vectors
+    */
+   void ArrayMult(const Array<const Vector *> &X,
+                  Array<Vector *> &Y) const override;
+
+   ~CuDSSSolver();
+
+private:
+#ifdef MFEM_USE_MPI
+   // MPI_Comm
+   MPI_Comm mpi_comm = MPI_COMM_NULL;
+
+   int row_start = 0;  // the first row index in CSR matrix operator
+   int row_end = 0;    // the end row index in CSR matrix operator
+#endif
+
+   // Parameter controlling whether or not to reuse the symbolic factorization
+   // for multiple calls to SetOperator
+   bool reorder_reuse = false;
+
+   // Parameter controlling the matrix type
+   cudssMatrixType_t mat_type = CUDSS_MTYPE_GENERAL;
+
+   int n_global = 0;      // global number of rows
+   int n_loc = 0;         // the number of the rows in CSR matrix operator
+
+   mutable int nrhs = 0;  // the number of the RHSs
+   int nnz = 0;        // the number of non zeros
+
+   // copy and keep the I and J arrays in device memory when skipping analysis
+   // phase
+   void *csr_offsets_d = NULL; // copy and keep I in device
+   void *csr_columns_d = NULL; // copy and keep J in device
+   void *csr_values_d = NULL;  // copy and keep csr data in device
+
+   // cuDSS object specifies available matrix types for sparse matrices
+   cudssMatrixViewType_t mview = CUDSS_MVIEW_FULL;
+
+   // cuDSS objects storage for sparse matrix Ac, RHS yc and solution xc
+   std::unique_ptr<cudssMatrix_t> Ac;
+   mutable cudssMatrix_t xc, yc;
+
+   // common for all cuDSS solver instances.
+   // cuDSS object holds the cuDSS library context
+   cudssHandle_t handle;
+
+   // cuDSS object stores configuration settings for the solver
+   mutable cudssConfig_t solverConfig;
+   // cuDSS object holds internal data
+   mutable cudssData_t solverData;
+
+   /// Method for configuring storage for distributed/centralized RHS and
+   /// solution
+   void SetNumRHS(int nrhs_) const;
+
+#ifdef MFEM_USE_MPI
+   /**
+    * @brief Set the HypreParMatrix object
+    *
+    * @param op HypreParMatrix object
+    *
+    * @note This method is called inside SetOperator
+    */
+   void SetMatrix(const HypreParMatrix &op);
+#endif
+
+   /**
+    * @brief Set the SparseMatrix object
+    *
+    * @param op SparseMatrix object
+    *
+    * @note This method is called inside SetOperator
+    */
+   void SetMatrix(const SparseMatrix &op);
+
+   /**
+    * @brief Set the matrix values for cuDSS
+    *
+    * @param csr_offsets Row offsets of the CSR matrix
+    * @param csr_columns Column indices of the CSR matrix
+    * @param csr_values Non-zero values of the CSR matrix
+    *
+    * @note This method is called inside SetMatrix.
+   */
+   void SetMatrixCuDSS(int* csr_offsets, int* csr_columns, real_t* csr_values);
+
+   /// Method for initializing the cuDSS library and creating the cuDSS handle
+   void InitCuDSS();
+};
+
+} // namespace mfem
+
+#endif // MFEM_USE_CUDSS
+#endif // MFEM_CUDSS
@@ -91,4 +91,8 @@

 #endif // MFEM_USE_MPI

+#ifdef MFEM_USE_CUDSS
+#include "cudss.hpp"
+#endif
+
 #endif
@@ -20,9 +20,7 @@
 #include "../fem/fem.hpp"

 #include "petsc.h"
-#if defined(PETSC_HAVE_HYPRE)
 #include "petscmathypre.h"
-#endif

 // Backward compatibility
 #if PETSC_VERSION_LT(3,11,0)
@@ -105,14 +103,6 @@ static PetscErrorCode MakeShellPC(PC,mfem::Solver&,bool);
 static PetscErrorCode MakeShellPCWithFactory(PC,
                                             mfem::PetscPreconditionerFactory*);

-// Equivalent functions are present in PETSc source code
-// if PETSc has been compiled with hypre support
-// We provide them here in case PETSC_HAVE_HYPRE is not defined
-#if !defined(PETSC_HAVE_HYPRE)
-static PetscErrorCode MatConvert_hypreParCSR_AIJ(hypre_ParCSRMatrix*,Mat*);
-static PetscErrorCode MatConvert_hypreParCSR_IS(hypre_ParCSRMatrix*,Mat*);
-#endif
-
 #if PETSC_VERSION_GE(3,15,0) && defined(PETSC_HAVE_DEVICE)
 #if defined(MFEM_USE_CUDA) && defined(PETSC_HAVE_CUDA)
 #define _USE_DEVICE
@@ -1149,11 +1139,10 @@ PetscParMatrix& PetscParMatrix::operator=(const HypreParMatrix& B)
   }
   height = B.Height();
   width  = B.Width();
-#if defined(PETSC_HAVE_HYPRE)
+
   ierr = MatCreateFromParCSR(B,MATAIJ,PETSC_USE_POINTER,&A);
-#else
-   ierr = MatConvert_hypreParCSR_AIJ(B,&A); CCHKERRQ(B.GetComm(),ierr);
-#endif
+   CCHKERRQ(B.GetComm(),ierr);
+
   SetUpForDevice();
   return *this;
 }
@@ -1529,7 +1518,6 @@ void PetscParMatrix::ConvertOperator(MPI_Comm comm, const Operator &op, Mat* A,
      }
      else if (tid == PETSC_MATHYPRE)
      {
-#if defined(PETSC_HAVE_HYPRE)
         if (istrans)
         {
            Mat B;
@@ -1541,9 +1529,6 @@ void PetscParMatrix::ConvertOperator(MPI_Comm comm, const Operator &op, Mat* A,
         {
            ierr = MatConvert(pA->A,MATHYPRE,MAT_INITIAL_MATRIX,A); PCHKERRQ(pA->A,ierr);
         }
-#else
-         MFEM_ABORT("Reconfigure PETSc with --download-hypre or --with-hypre")
-#endif
      }
      else if (tid == PETSC_MATSHELL)
      {
@@ -1558,33 +1543,21 @@ void PetscParMatrix::ConvertOperator(MPI_Comm comm, const Operator &op, Mat* A,
   {
      if (tid == PETSC_MATAIJ)
      {
-#if defined(PETSC_HAVE_HYPRE)
         ierr = MatCreateFromParCSR(const_cast<HypreParMatrix&>(*pH),MATAIJ,
                                    PETSC_USE_POINTER,A);
-#else
-         ierr = MatConvert_hypreParCSR_AIJ(const_cast<HypreParMatrix&>(*pH),A);
-#endif
         CCHKERRQ(pH->GetComm(),ierr);
      }
      else if (tid == PETSC_MATIS)
      {
-#if defined(PETSC_HAVE_HYPRE)
         ierr = MatCreateFromParCSR(const_cast<HypreParMatrix&>(*pH),MATIS,
                                    PETSC_USE_POINTER,A);
-#else
-         ierr = MatConvert_hypreParCSR_IS(const_cast<HypreParMatrix&>(*pH),A);
-#endif
         CCHKERRQ(pH->GetComm(),ierr);
      }
      else if (tid == PETSC_MATHYPRE || tid == ANY_TYPE)
      {
-#if defined(PETSC_HAVE_HYPRE)
         ierr = MatCreateFromParCSR(const_cast<HypreParMatrix&>(*pH),MATHYPRE,
                                    PETSC_USE_POINTER,A);
         CCHKERRQ(pH->GetComm(),ierr);
-#else
-         MFEM_ABORT("Reconfigure PETSc with --download-hypre or --with-hypre")
-#endif
      }
      else if (tid == PETSC_MATSHELL)
      {
@@ -2238,11 +2211,7 @@ PetscParMatrix * RAP(PetscParMatrix *A, PetscParMatrix *P)
 PetscParMatrix * RAP(HypreParMatrix *hA, PetscParMatrix *P)
 {
   PetscParMatrix *out,*A;
-#if defined(PETSC_HAVE_HYPRE)
   A = new PetscParMatrix(hA,Operator::PETSC_MATHYPRE);
-#else
-   A = new PetscParMatrix(hA);
-#endif
   out = RAP(P,A,P);
   delete A;
   return out;
@@ -2349,10 +2318,8 @@ Operator::Type PetscParMatrix::GetType() const
   if (ok == PETSC_TRUE) { return PETSC_MATSHELL; }
   ierr = PetscObjectTypeCompare(oA, MATNEST, &ok); PCHKERRQ(A,ierr);
   if (ok == PETSC_TRUE) { return PETSC_MATNEST; }
-#if defined(PETSC_HAVE_HYPRE)
   ierr = PetscObjectTypeCompare(oA, MATHYPRE, &ok); PCHKERRQ(A,ierr);
   if (ok == PETSC_TRUE) { return PETSC_MATHYPRE; }
-#endif
   return PETSC_MATGENERIC;
 }

@@ -5527,207 +5494,6 @@ static PetscErrorCode Convert_Vmarks_IS(MPI_Comm comm,
   PetscFunctionReturn(PETSC_SUCCESS);
 }

-#if !defined(PETSC_HAVE_HYPRE)
-
-#if defined(HYPRE_MIXEDINT)
-#error "HYPRE_MIXEDINT not supported"
-#endif
-
-#include "_hypre_parcsr_mv.h"
-static PetscErrorCode MatConvert_hypreParCSR_AIJ(hypre_ParCSRMatrix* hA,Mat* pA)
-{
-   MPI_Comm        comm = hypre_ParCSRMatrixComm(hA);
-   hypre_CSRMatrix *hdiag,*hoffd;
-   PetscScalar     *da,*oa,*aptr;
-   PetscInt        *dii,*djj,*oii,*ojj,*iptr;
-   PetscInt        i,dnnz,onnz,m,n;
-   PetscMPIInt     size;
-   PetscErrorCode  ierr;
-
-   PetscFunctionBeginUser;
-   hdiag = hypre_ParCSRMatrixDiag(hA);
-   hoffd = hypre_ParCSRMatrixOffd(hA);
-   m     = hypre_CSRMatrixNumRows(hdiag);
-   n     = hypre_CSRMatrixNumCols(hdiag);
-   dnnz  = hypre_CSRMatrixNumNonzeros(hdiag);
-   onnz  = hypre_CSRMatrixNumNonzeros(hoffd);
-   ierr  = PetscMalloc1(m+1,&dii); CHKERRQ(ierr);
-   ierr  = PetscMalloc1(dnnz,&djj); CHKERRQ(ierr);
-   ierr  = PetscMalloc1(dnnz,&da); CHKERRQ(ierr);
-   ierr  = PetscMemcpy(dii,hypre_CSRMatrixI(hdiag),(m+1)*sizeof(PetscInt));
-   CHKERRQ(ierr);
-   ierr  = PetscMemcpy(djj,hypre_CSRMatrixJ(hdiag),dnnz*sizeof(PetscInt));
-   CHKERRQ(ierr);
-   ierr  = PetscMemcpy(da,hypre_CSRMatrixData(hdiag),dnnz*sizeof(PetscScalar));
-   CHKERRQ(ierr);
-   iptr  = djj;
-   aptr  = da;
-   for (i=0; i<m; i++)
-   {
-      PetscInt nc = dii[i+1]-dii[i];
-      ierr = PetscSortIntWithScalarArray(nc,iptr,aptr); CHKERRQ(ierr);
-      iptr += nc;
-      aptr += nc;
-   }
-   ierr = MPI_Comm_size(comm,&size); CHKERRQ(ierr);
-   if (size > 1)
-   {
-      PetscInt *offdj,*coffd;
-
-      ierr  = PetscMalloc1(m+1,&oii); CHKERRQ(ierr);
-      ierr  = PetscMalloc1(onnz,&ojj); CHKERRQ(ierr);
-      ierr  = PetscMalloc1(onnz,&oa); CHKERRQ(ierr);
-      ierr  = PetscMemcpy(oii,hypre_CSRMatrixI(hoffd),(m+1)*sizeof(PetscInt));
-      CHKERRQ(ierr);
-      offdj = hypre_CSRMatrixJ(hoffd);
-      coffd = hypre_ParCSRMatrixColMapOffd(hA);
-      for (i=0; i<onnz; i++) { ojj[i] = coffd[offdj[i]]; }
-      ierr  = PetscMemcpy(oa,hypre_CSRMatrixData(hoffd),onnz*sizeof(PetscScalar));
-      CHKERRQ(ierr);
-      iptr  = ojj;
-      aptr  = oa;
-      for (i=0; i<m; i++)
-      {
-         PetscInt nc = oii[i+1]-oii[i];
-         ierr = PetscSortIntWithScalarArray(nc,iptr,aptr); CHKERRQ(ierr);
-         iptr += nc;
-         aptr += nc;
-      }
-      ierr = MatCreateMPIAIJWithSplitArrays(comm,m,n,PETSC_DECIDE,PETSC_DECIDE,dii,
-                                            djj,da,oii,ojj,oa,pA); CHKERRQ(ierr);
-   }
-   else
-   {
-      oii = ojj = NULL;
-      oa = NULL;
-      ierr = MatCreateSeqAIJWithArrays(comm,m,n,dii,djj,da,pA); CHKERRQ(ierr);
-   }
-   /* We are responsible to free the CSR arrays.  However, since we can take
-      references of a PetscParMatrix but we cannot take reference of PETSc
-      arrays, we need to create a PetscContainer object to take reference of
-      these arrays in reference objects */
-   void *ptrs[6] = {dii,djj,da,oii,ojj,oa};
-   const char *names[6] = {"_mfem_csr_dii",
-                           "_mfem_csr_djj",
-                           "_mfem_csr_da",
-                           "_mfem_csr_oii",
-                           "_mfem_csr_ojj",
-                           "_mfem_csr_oa"
-                          };
-   for (i=0; i<6; i++)
-   {
-      PetscContainer c;
-
-      ierr = PetscContainerCreate(comm,&c); CHKERRQ(ierr);
-      ierr = PetscContainerSetPointer(c,ptrs[i]); CHKERRQ(ierr);
-      ierr = PetscContainerSetCtxDestroy(c,__mfem_array_container_destroy);
-      CHKERRQ(ierr);
-      ierr = PetscObjectCompose((PetscObject)(*pA),names[i],(PetscObject)c);
-      CHKERRQ(ierr);
-      ierr = PetscContainerDestroy(&c); CHKERRQ(ierr);
-   }
-   PetscFunctionReturn(PETSC_SUCCESS);
-}
-
-static PetscErrorCode MatConvert_hypreParCSR_IS(hypre_ParCSRMatrix* hA,Mat* pA)
-{
-   Mat                    lA;
-   ISLocalToGlobalMapping rl2g,cl2g;
-   IS                     is;
-   hypre_CSRMatrix        *hdiag,*hoffd;
-   MPI_Comm               comm = hypre_ParCSRMatrixComm(hA);
-   void                   *ptrs[2];
-   const char             *names[2] = {"_mfem_csr_aux",
-                                       "_mfem_csr_data"
-                                      };
-   PetscScalar            *hdd,*hod,*aa,*data;
-   PetscInt               *col_map_offd,*hdi,*hdj,*hoi,*hoj;
-   PetscInt               *aux,*ii,*jj;
-   PetscInt               cum,dr,dc,oc,str,stc,nnz,i,jd,jo;
-   PetscErrorCode         ierr;
-
-   PetscFunctionBeginUser;
-   /* access relevant information in ParCSR */
-   str   = hypre_ParCSRMatrixFirstRowIndex(hA);
-   stc   = hypre_ParCSRMatrixFirstColDiag(hA);
-   hdiag = hypre_ParCSRMatrixDiag(hA);
-   hoffd = hypre_ParCSRMatrixOffd(hA);
-   dr    = hypre_CSRMatrixNumRows(hdiag);
-   dc    = hypre_CSRMatrixNumCols(hdiag);
-   nnz   = hypre_CSRMatrixNumNonzeros(hdiag);
-   hdi   = hypre_CSRMatrixI(hdiag);
-   hdj   = hypre_CSRMatrixJ(hdiag);
-   hdd   = hypre_CSRMatrixData(hdiag);
-   oc    = hypre_CSRMatrixNumCols(hoffd);
-   nnz  += hypre_CSRMatrixNumNonzeros(hoffd);
-   hoi   = hypre_CSRMatrixI(hoffd);
-   hoj   = hypre_CSRMatrixJ(hoffd);
-   hod   = hypre_CSRMatrixData(hoffd);
-
-   /* generate l2g maps for rows and cols */
-   ierr = ISCreateStride(comm,dr,str,1,&is); CHKERRQ(ierr);
-   ierr = ISLocalToGlobalMappingCreateIS(is,&rl2g); CHKERRQ(ierr);
-   ierr = ISDestroy(&is); CHKERRQ(ierr);
-   col_map_offd = hypre_ParCSRMatrixColMapOffd(hA);
-   ierr = PetscMalloc1(dc+oc,&aux); CHKERRQ(ierr);
-   for (i=0; i<dc; i++) { aux[i]    = i+stc; }
-   for (i=0; i<oc; i++) { aux[i+dc] = col_map_offd[i]; }
-   ierr = ISCreateGeneral(comm,dc+oc,aux,PETSC_OWN_POINTER,&is); CHKERRQ(ierr);
-   ierr = ISLocalToGlobalMappingCreateIS(is,&cl2g); CHKERRQ(ierr);
-   ierr = ISDestroy(&is); CHKERRQ(ierr);
-
-   /* create MATIS object */
-   ierr = MatCreate(comm,pA); CHKERRQ(ierr);
-   ierr = MatSetSizes(*pA,dr,dc,PETSC_DECIDE,PETSC_DECIDE); CHKERRQ(ierr);
-   ierr = MatSetType(*pA,MATIS); CHKERRQ(ierr);
-   ierr = MatSetLocalToGlobalMapping(*pA,rl2g,cl2g); CHKERRQ(ierr);
-   ierr = ISLocalToGlobalMappingDestroy(&rl2g); CHKERRQ(ierr);
-   ierr = ISLocalToGlobalMappingDestroy(&cl2g); CHKERRQ(ierr);
-
-   /* merge local matrices */
-   ierr = PetscMalloc1(nnz+dr+1,&aux); CHKERRQ(ierr);
-   ierr = PetscMalloc1(nnz,&data); CHKERRQ(ierr);
-   ii   = aux;
-   jj   = aux+dr+1;
-   aa   = data;
-   *ii  = *(hdi++) + *(hoi++);
-   for (jd=0,jo=0,cum=0; *ii<nnz; cum++)
-   {
-      PetscScalar *aold = aa;
-      PetscInt    *jold = jj,nc = jd+jo;
-      for (; jd<*hdi; jd++) { *jj++ = *hdj++;      *aa++ = *hdd++; }
-      for (; jo<*hoi; jo++) { *jj++ = *hoj++ + dc; *aa++ = *hod++; }
-      *(++ii) = *(hdi++) + *(hoi++);
-      ierr = PetscSortIntWithScalarArray(jd+jo-nc,jold,aold); CHKERRQ(ierr);
-   }
-   for (; cum<dr; cum++) { *(++ii) = nnz; }
-   ii   = aux;
-   jj   = aux+dr+1;
-   aa   = data;
-   ierr = MatCreateSeqAIJWithArrays(PETSC_COMM_SELF,dr,dc+oc,ii,jj,aa,&lA);
-   CHKERRQ(ierr);
-   ptrs[0] = aux;
-   ptrs[1] = data;
-   for (i=0; i<2; i++)
-   {
-      PetscContainer c;
-
-      ierr = PetscContainerCreate(PETSC_COMM_SELF,&c); CHKERRQ(ierr);
-      ierr = PetscContainerSetPointer(c,ptrs[i]); CHKERRQ(ierr);
-      ierr = PetscContainerSetCtxDestroy(c,__mfem_array_container_destroy);
-      CHKERRQ(ierr);
-      ierr = PetscObjectCompose((PetscObject)lA,names[i],(PetscObject)c);
-      CHKERRQ(ierr);
-      ierr = PetscContainerDestroy(&c); CHKERRQ(ierr);
-   }
-   ierr = MatISSetLocalMat(*pA,lA); CHKERRQ(ierr);
-   ierr = MatDestroy(&lA); CHKERRQ(ierr);
-   ierr = MatAssemblyBegin(*pA,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
-   ierr = MatAssemblyEnd(*pA,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
-   PetscFunctionReturn(PETSC_SUCCESS);
-}
-#endif
-
 #include <petsc/private/matimpl.h>

 static PetscErrorCode __mfem_MatCreateDummy(MPI_Comm comm, PetscInt m,
@@ -42,6 +42,9 @@
 #if !defined(PETSC_USE_64BIT_INDICES) && (defined(HYPRE_BIGINT) || defined(HYPRE_MIXEDINT))
 #error "Mismatch between HYPRE (64bit) and PETSc (32bit) integer types"
 #endif
+#if !defined(PETSC_HAVE_HYPRE)
+#error "MFEM requires PETSc built with HYPRE support"
+#endif

 #include "petscversion.h"
 #if PETSC_VERSION_GE(3,12,0)
@@ -307,7 +307,7 @@ endif
 MFEM_REQ_LIB_DEPS = SUPERLU MUMPS METIS FMS CONDUIT SIDRE LAPACK SUNDIALS\
 SUITESPARSE STRUMPACK GINKGO GNUTLS HDF5 NETCDF SLEPC PETSC MPFR PUMI HIOP\
 GSLIB OCCA CEED RAJA UMPIRE MKL_CPARDISO MKL_PARDISO AMGX MAGMA CALIPER PARELAG\
- TRIBOL BENCHMARK MOONOLITH ALGOIM
+ TRIBOL BENCHMARK MOONOLITH ALGOIM CUDSS


 PETSC_ERROR_MSG = $(if $(PETSC_FOUND),,. PETSC config not found: $(PETSC_VARS))
@@ -376,7 +376,8 @@ MFEM_DEFINES = MFEM_VERSION MFEM_VERSION_STRING MFEM_GIT_STRING MFEM_USE_MPI\
 MFEM_USE_SIMD MFEM_USE_ADIOS2 MFEM_USE_MKL_CPARDISO MFEM_USE_MKL_PARDISO MFEM_USE_AMGX\
 MFEM_USE_MAGMA MFEM_USE_MUMPS MFEM_USE_ADFORWARD MFEM_USE_CODIPACK MFEM_USE_CALIPER\
 MFEM_USE_BENCHMARK MFEM_USE_PARELAG MFEM_USE_TRIBOL MFEM_USE_ALGOIM MFEM_USE_ENZYME\
- MFEM_SOURCE_DIR MFEM_INSTALL_DIR MFEM_SHARED_BUILD MFEM_USE_DOUBLE MFEM_USE_SINGLE
+ MFEM_SOURCE_DIR MFEM_INSTALL_DIR MFEM_SHARED_BUILD MFEM_USE_DOUBLE MFEM_USE_SINGLE\
+ MFEM_USE_CUDSS MFEM_CUDSS_COMM_LIB MFEM_CUDSS_THREADING_LIB

 # List of makefile variables that will be written to config.mk:
 MFEM_CONFIG_VARS = MFEM_CXX MFEM_HOST_CXX MFEM_CPPFLAGS MFEM_CXXFLAGS\
@@ -411,7 +412,7 @@ MFEM_INSTALL_DIR = $(abspath $(MFEM_PREFIX))
 # If we have 'config' target, export variables used by config/makefile
 ifneq (,$(filter config,$(MAKECMDGOALS)))
   export $(MFEM_DEFINES) MFEM_DEFINES $(MFEM_CONFIG_VARS) MFEM_CONFIG_VARS
-   export VERBOSE HYPRE_OPT PUMI_DIR MUMPS_OPT GSLIB_OPT
+   export VERBOSE HYPRE_OPT PUMI_DIR MUMPS_OPT GSLIB_OPT CUDSS_OPT
 endif

 # If we have 'install' target, export variables used by config/makefile
@@ -742,6 +743,7 @@ status info:
 	$(info MFEM_USE_SUPERLU5      = $(MFEM_USE_SUPERLU5))
 	$(info MFEM_USE_MUMPS         = $(MFEM_USE_MUMPS))
 	$(info MFEM_USE_STRUMPACK     = $(MFEM_USE_STRUMPACK))
+	$(info MFEM_USE_CUDSS         = $(MFEM_USE_CUDSS))
 	$(info MFEM_USE_GINKGO        = $(MFEM_USE_GINKGO))
 	$(info MFEM_USE_AMGX          = $(MFEM_USE_AMGX))
 	$(info MFEM_USE_MAGMA         = $(MFEM_USE_MAGMA))
@@ -203,6 +203,69 @@ void Mesh::GetBoundingBox(Vector &min, Vector &max, int ref)
   }
 }

+void Mesh::GetElementBoundingBoxes(Vector &elmin, Vector &elmax, int nref)
+{
+   int nel = GetNE();
+   elmin.SetSize(spaceDim * nel);
+   elmax.SetSize(spaceDim * nel);
+   elmin = numeric_limits<real_t>::max();
+   elmax = -numeric_limits<real_t>::max();
+   if (Nodes == NULL)
+   {
+      Array<int> verts;
+      real_t *coord;
+      // create bounding boxes from vertex coordinates
+      for (int e = 0; e < nel; e++)
+      {
+         GetElementVertices(e, verts);
+         for (int v = 0; v < verts.Size(); v++)
+         {
+            coord = GetVertex(verts[v]);
+            for (int d = 0; d < spaceDim; d++)
+            {
+               elmin(d*nel + e) = min(elmin(d*nel + e), coord[d]);
+               elmax(d*nel + e) = max(elmax(d*nel + e), coord[d]);
+            }
+         }
+      }
+   }
+   else
+   {
+      Nodes->GetElementBounds(elmin, elmax, nref);
+   }
+}
+
+void Mesh::GetElementBoundingBox(int elem, Vector &elmin, Vector &elmax,
+                                 int nref)
+{
+   elmin.SetSize(spaceDim);
+   elmax.SetSize(spaceDim);
+   elmin = numeric_limits<real_t>::max();
+   elmax = -numeric_limits<real_t>::max();
+   if (Nodes == NULL)
+   {
+      Array<int> verts;
+      real_t *coord;
+      // create bounding boxes from vertex coordinates
+      GetElementVertices(elem, verts);
+      for (int v = 0; v < verts.Size(); v++)
+      {
+         coord = GetVertex(verts[v]);
+         for (int d = 0; d < spaceDim; d++)
+         {
+            elmin(d) = min(elmin(d), coord[d]);
+            elmax(d) = max(elmax(d), coord[d]);
+         }
+      }
+   }
+   else
+   {
+      const int max_order = Nodes->FESpace()->GetMaxElementOrder();
+      PLBound plb(Nodes->FESpace(), nref*(max_order+1));
+      Nodes->GetElementBounds(elem, plb, elmin, elmax);
+   }
+}
+
 void Mesh::GetCharacteristics(real_t &h_min, real_t &h_max,
                              real_t &kappa_min, real_t &kappa_max,
                              Vector *Vh, Vector *Vk)
@@ -1366,6 +1366,11 @@ public:
   /** For high-order meshes, the geometry is first refined @a ref times. */
   void GetBoundingBox(Vector &min, Vector &max, int ref = 2);

+   /// Returns the minimum and maximum corners of the element bounding box.
+   void GetElementBoundingBoxes(Vector &elmin, Vector &elmax, int nref=2);
+
+   void GetElementBoundingBox(int elem, Vector &elmin, Vector &elmax, int nref=2);
+
   void GetCharacteristics(real_t &h_min, real_t &h_max,
                           real_t &kappa_min, real_t &kappa_max,
                           Vector *Vh = NULL, Vector *Vk = NULL);
@@ -4838,6 +4838,8 @@ void ParMesh::Print(std::ostream &os, const std::string &comments) const
 {
   int shared_bdr_attr;
   Array<int> nc_shared_faces;
+   Array<int> interface_faces;
+   int interface_bdr_attr = 0;

   if (NURBSext)
   {
@@ -4882,6 +4884,19 @@ void ParMesh::Print(std::ostream &os, const std::string &comments) const
   const bool set_names = attribute_sets.SetsExist() ||
                          bdr_attribute_sets.SetsExist();

+   // Add material interfaces as boundary elements for visualization. We build a
+   // list of local faces to print as extra boundary elements. This does not
+   // modify the ParMesh object, only the printed mesh.
+   if (print_interfaces && Dim > 1)
+   {
+      FindInterface(interface_faces);
+
+      // Choose a boundary attribute that does not collide with existing ones,
+      // including those introduced by print_shared.
+      const int max_bdr_attr = bdr_attributes.Size() ? bdr_attributes.Max() : 0;
+      interface_bdr_attr = max_bdr_attr + 1 + (print_shared ? NRanks : 0);
+   }
+
   os << (!set_names ? "MFEM mesh v1.0\n" : "MFEM mesh v1.3\n");

   if (!comments.empty()) { os << '\n' << comments << '\n'; }
@@ -4916,6 +4931,11 @@ void ParMesh::Print(std::ostream &os, const std::string &comments) const
   {
      num_bdr_elems += s2l_face->Size();
   }
+   if (print_interfaces && Dim > 1)
+   {
+      // in 3D we print two oriented copies for each material interface face
+      num_bdr_elems += (Dim == 3 ? 2 : 1) * interface_faces.Size();
+   }
   os << "\nboundary\n" << num_bdr_elems << '\n';
   for (int i = 0; i < NumOfBdrElements; i++)
   {
@@ -4940,6 +4960,51 @@ void ParMesh::Print(std::ostream &os, const std::string &comments) const
      }
   }

+   // Print interface faces as additional boundary elements. In 3D we print two
+   // copies of each face, with opposite orientation, so material subdomains can
+   // be pulled apart with F11/F12 in GLVis.
+   if (print_interfaces && Dim > 1)
+   {
+      for (int i = 0; i < interface_faces.Size(); i++)
+      {
+         const int f = interface_faces[i];
+         const int *fv = faces[f]->GetVertices();
+         if (Dim == 2)
+         {
+            Segment seg(fv, interface_bdr_attr);
+            PrintElement(&seg, os);
+         }
+         else // Dim == 3
+         {
+            const Geometry::Type geom = faces[f]->GetGeometryType();
+            if (geom == Geometry::TRIANGLE)
+            {
+               int v0[3] = { fv[0], fv[1], fv[2] };
+               int v1[3] = { fv[0], fv[2], fv[1] };
+               Triangle t0(v0, interface_bdr_attr);
+               Triangle t1(v1, interface_bdr_attr);
+               PrintElement(&t0, os);
+               PrintElement(&t1, os);
+            }
+            else if (geom == Geometry::SQUARE)
+            {
+               int v0[4] = { fv[0], fv[1], fv[2], fv[3] };
+               int v1[4] = { fv[0], fv[3], fv[2], fv[1] };
+               Quadrilateral q0(v0, interface_bdr_attr);
+               Quadrilateral q1(v1, interface_bdr_attr);
+               PrintElement(&q0, os);
+               PrintElement(&q1, os);
+            }
+            else
+            {
+               MFEM_ABORT("unsupported 3D face geometry type '"
+                          << Geometry::Name[geom]
+                          << "' while printing interface boundaries.");
+            }
+         }
+      }
+   }
+
   if (set_names)
   {
      os << "\nbdr_attribute_sets\n";
@@ -4989,8 +5054,9 @@ void ParMesh::Print(adios2stream &os) const
 }
 #endif

-static void dump_element(const Element* elem, Array<int> &data)
+static void dump_element_with_attr(const Element* elem, Array<int> &data)
 {
+   data.Append(elem->GetAttribute());
   data.Append(elem->GetGeometryType());

   int nv = elem->GetNVertices();
@@ -5001,13 +5067,51 @@ static void dump_element(const Element* elem, Array<int> &data)
   }
 }

+void ParMesh::FindInterface(Array<int> &interface) const
+{
+   // We need face neighbor elements to determine if shared faces in parallel
+   // are on material interfaces.
+   const_cast<ParMesh*>(this)->ExchangeFaceNbrData();
+
+   interface.SetSize(0); // clear 'interface' since we Append to it
+   const int nf = GetNumFaces();
+   for (int f = 0; f < nf; f++)
+   {
+      if (!FaceIsTrueInterior(f)) { continue; } // skip true boundary
+
+      const int e1 = faces_info[f].Elem1No;
+      if (e1 < 0) { continue; }
+      const int a1 = elements[e1]->GetAttribute();
+
+      int a2 = a1;
+      if (faces_info[f].Elem2No >= 0)
+      {
+         a2 = elements[faces_info[f].Elem2No]->GetAttribute();
+      }
+      else
+      {
+         // Shared face: element 2 is a face-neighbor element with index
+         // -1-Elem2No, i.e. FlipIndexSign(Elem2No) (see Mesh::FaceInfo).
+         const int nbr_el = FlipIndexSign(faces_info[f].Elem2No);
+         MFEM_ASSERT(0 <= nbr_el && nbr_el < face_nbr_elements.Size(),
+                     "invalid face-neighbor index");
+         a2 = face_nbr_elements[nbr_el]->GetAttribute();
+      }
+
+      if (a1 != a2) { interface.Append(f); }
+   }
+}
+
 void ParMesh::PrintAsOne(std::ostream &os, const std::string &comments) const
 {
-   int i, j, k, p, nv_ne[2], &nv = nv_ne[0], &ne = nv_ne[1], vc;
+   int i, j, k, p, nv_ne[2], &nv = nv_ne[0], &ne = nv_ne[1];
+   long long vc; // global vertex offset
   const int *v;
   MPI_Status status;
   Array<real_t> vert;
   Array<int> ints;
+   Array<int> interface_faces;
+   int interface_bdr_attr = 0;

   if (MyRank == 0)
   {
@@ -5030,15 +5134,20 @@ void ParMesh::PrintAsOne(std::ostream &os, const std::string &comments) const
      os << "\ndimension\n" << Dim;
   }

-   nv = NumOfElements;
-   MPI_Reduce(&nv, &ne, 1, MPI_INT, MPI_SUM, 0, MyComm);
+   long long loc_ne = NumOfElements, glob_ne = 0;
+   MPI_Reduce(&loc_ne, &glob_ne, 1, MPI_LONG_LONG, MPI_SUM, 0, MyComm);
   if (MyRank == 0)
   {
-      os << "\n\nelements\n" << ne << '\n';
+      MFEM_VERIFY(static_cast<int>(glob_ne) == glob_ne,
+                  "integer overflow detected!");
+      os << "\n\nelements\n" << glob_ne << '\n';
      for (i = 0; i < NumOfElements; i++)
      {
-         // processor number + 1 as attribute and geometry type
-         os << 1 << ' ' << elements[i]->GetGeometryType();
+         // Print attribute + geometry:
+         // * if print_shared != 0, use processor number + 1 as attribute
+         // * otherwise, use the real attribute
+         os << (print_shared ? 1 : GetAttribute(i))
+            << ' ' << elements[i]->GetGeometryType();
         // vertices
         nv = elements[i]->GetNVertices();
         v  = elements[i]->GetVertices();
@@ -5048,19 +5157,26 @@ void ParMesh::PrintAsOne(std::ostream &os, const std::string &comments) const
         }
         os << '\n';
      }
-      vc = NumOfVertices;
+      vc = NumOfVertices;  // global offset for vertex indices
      for (p = 1; p < NRanks; p++)
      {
         MPI_Recv(nv_ne, 2, MPI_INT, p, 444, MyComm, &status);
         ints.SetSize(ne);
         if (ne)
         {
+            // Receive an array that contains attribute + geometry + vertices
+            // for each element.
            MPI_Recv(&ints[0], ne, MPI_INT, p, 445, MyComm, &status);
         }
+         MFEM_VERIFY(static_cast<int>(vc + nv) == (vc + nv),
+                     "integer overflow detected!");
         for (i = 0; i < ne; )
         {
-            // processor number + 1 as attribute and geometry type
-            os << p+1 << ' ' << ints[i];
+            // Print attribute + geometry:
+            // * if print_shared != 0, use processor number + 1 as attribute
+            // * otherwise, use the real attribute
+            os << (print_shared ? p+1 : ints[i]) << ' ' << ints[i+1];
+            i++;
            // vertices
            k = Geometries.GetVertices(ints[i++])->GetNPoints();
            for (j = 0; j < k; j++)
@@ -5074,11 +5190,11 @@ void ParMesh::PrintAsOne(std::ostream &os, const std::string &comments) const
   }
   else
   {
-      // for each element send its geometry type and its vertices
+      // for each element send its attribute, geometry type and vertices
      ne = 0;
      for (i = 0; i < NumOfElements; i++)
      {
-         ne += 1 + elements[i]->GetNVertices();
+         ne += 2 + elements[i]->GetNVertices(); // attribute + geom + vertices
      }
      nv = NumOfVertices;
      MPI_Send(nv_ne, 2, MPI_INT, 0, 444, MyComm);
@@ -5087,7 +5203,7 @@ void ParMesh::PrintAsOne(std::ostream &os, const std::string &comments) const
      ints.SetSize(0);
      for (i = 0; i < NumOfElements; i++)
      {
-         dump_element(elements[i], ints);
+         dump_element_with_attr(elements[i], ints);
      }
      MFEM_ASSERT(ints.Size() == ne, "");
      if (ne)
@@ -5096,36 +5212,57 @@ void ParMesh::PrintAsOne(std::ostream &os, const std::string &comments) const
      }
   }

-   // boundary + shared boundary
+   // Add material interfaces as boundary elements for visualization. We build a
+   // list of local faces to print as extra boundary elements. This does not
+   // modify the ParMesh object, only the printed mesh.
+   if (print_interfaces && Dim > 1)
+   {
+      FindInterface(interface_faces);
+
+      // Choose a boundary attribute that does not collide with existing ones.
+      const int max_bdr_attr = bdr_attributes.Size() ? bdr_attributes.Max() : 0;
+      // If print_shared is enabled, this attribute will be replaced:
+      interface_bdr_attr = max_bdr_attr + 1;
+   }
+
+   // boundary + (optionally) shared faces + (optionally) interface faces
   ne = NumOfBdrElements;
-   if (!pncmesh)
+   if (print_shared && !pncmesh)
   {
      ne += GetNSharedFaces();
   }
-   else if (Dim > 1)
+   if (print_shared && pncmesh && Dim > 1)
   {
      const NCMesh::NCList &list = pncmesh->GetSharedList(Dim - 1);
      ne += list.conforming.Size() + list.masters.Size() + list.slaves.Size();
      // In addition to the number returned by GetNSharedFaces(), include the
      // the master shared faces as well.
   }
-   ints.Reserve(ne * (1 + 2*(Dim-1))); // just an upper bound
+   if (print_interfaces && Dim > 1)
+   {
+      // In 3D we print two oriented copies for each material interface face:
+      ne += (Dim == 3 ? 2 : 1) * interface_faces.Size();
+   }
+   ints.Reserve(ne * (2 + (1 << Dim))); // just an upper bound
   ints.SetSize(0);

-   // for each boundary and shared boundary element send its geometry type
-   // and its vertices
+   // For each boundary, (optionally) shared face, and (optionally) interface
+   // face send its attribute, geometry type, and vertices.
   ne = 0;
-   for (i = j = 0; i < NumOfBdrElements; i++)
+   for (i = 0; i < NumOfBdrElements; i++)
   {
-      dump_element(boundary[i], ints); ne++;
+      dump_element_with_attr(boundary[i], ints); ne++;
   }
-   if (!pncmesh)
+   if (print_shared && !pncmesh)
   {
+      // Attribute for visualized shared faces on parallel interfaces:
+      constexpr int shared_attribute = 1;
      switch (Dim)
      {
         case 1:
            for (i = 0; i < svert_lvert.Size(); i++)
            {
+               ints.Append(shared_attribute);
               ints.Append(Geometry::POINT);
               ints.Append(svert_lvert[i]);
               ne++;
@@ -5135,19 +5272,21 @@ void ParMesh::PrintAsOne(std::ostream &os, const std::string &comments) const
         case 2:
            for (i = 0; i < shared_edges.Size(); i++)
            {
-               dump_element(shared_edges[i], ints); ne++;
+               dump_element_with_attr(shared_edges[i], ints); ne++;
            }
            break;

         case 3:
            for (i = 0; i < shared_trias.Size(); i++)
            {
+               ints.Append(shared_attribute);
               ints.Append(Geometry::TRIANGLE);
               ints.Append(shared_trias[i].v, 3);
               ne++;
            }
            for (i = 0; i < shared_quads.Size(); i++)
            {
+               ints.Append(shared_attribute);
               ints.Append(Geometry::SQUARE);
               ints.Append(shared_quads[i].v, 4);
               ne++;
@@ -5158,32 +5297,85 @@ void ParMesh::PrintAsOne(std::ostream &os, const std::string &comments) const
            MFEM_ABORT("invalid dimension: " << Dim);
      }
   }
-   else if (Dim > 1)
+   if (print_shared && pncmesh && Dim > 1)
   {
      const NCMesh::NCList &list = pncmesh->GetSharedList(Dim - 1);
      const int nfaces = GetNumFaces();
      for (i = 0; i < list.conforming.Size(); i++)
      {
         int index = list.conforming[i].index;
-         if (index < nfaces) { dump_element(faces[index], ints); ne++; }
+         if (index < nfaces)
+         {
+            dump_element_with_attr(faces[index], ints); ne++;
+         }
      }
      for (i = 0; i < list.masters.Size(); i++)
      {
         int index = list.masters[i].index;
-         if (index < nfaces) { dump_element(faces[index], ints); ne++; }
+         if (index < nfaces)
+         {
+            dump_element_with_attr(faces[index], ints); ne++;
+         }
      }
      for (i = 0; i < list.slaves.Size(); i++)
      {
         int index = list.slaves[i].index;
-         if (index < nfaces) { dump_element(faces[index], ints); ne++; }
+         if (index < nfaces)
+         {
+            dump_element_with_attr(faces[index], ints); ne++;
+         }
+      }
+   }
+   if (print_interfaces && Dim > 1)
+   {
+      for (i = 0; i < interface_faces.Size(); i++)
+      {
+         const int f = interface_faces[i];
+         const int *fv = faces[f]->GetVertices();
+         if (Dim == 2)
+         {
+            Segment seg(fv, interface_bdr_attr);
+            dump_element_with_attr(&seg, ints); ne++;
+         }
+         else // Dim == 3
+         {
+            const Geometry::Type geom = faces[f]->GetGeometryType();
+            if (geom == Geometry::TRIANGLE)
+            {
+               int v0[3] = { fv[0], fv[1], fv[2] };
+               int v1[3] = { fv[0], fv[2], fv[1] };
+               Triangle t0(v0, interface_bdr_attr);
+               Triangle t1(v1, interface_bdr_attr);
+               dump_element_with_attr(&t0, ints); ne++;
+               dump_element_with_attr(&t1, ints); ne++;
+            }
+            else if (geom == Geometry::SQUARE)
+            {
+               int v0[4] = { fv[0], fv[1], fv[2], fv[3] };
+               int v1[4] = { fv[0], fv[3], fv[2], fv[1] };
+               Quadrilateral q0(v0, interface_bdr_attr);
+               Quadrilateral q1(v1, interface_bdr_attr);
+               dump_element_with_attr(&q0, ints); ne++;
+               dump_element_with_attr(&q1, ints); ne++;
+            }
+            else
+            {
+               MFEM_ABORT("unsupported 3D face geometry type '"
+                          << Geometry::Name[geom]
+                          << "' while printing interface boundaries.");
+            }
+         }
      }
   }

-   MPI_Reduce(&ne, &k, 1, MPI_INT, MPI_SUM, 0, MyComm);
+   long long loc_nb = ne, glob_nb = 0;
+   MPI_Reduce(&loc_nb, &glob_nb, 1, MPI_LONG_LONG, MPI_SUM, 0, MyComm);
   if (MyRank == 0)
   {
-      os << "\nboundary\n" << k << '\n';
-      vc = 0;
+      MFEM_VERIFY(static_cast<int>(glob_nb) == glob_nb,
+                  "integer overflow detected!");
+      os << "\nboundary\n" << glob_nb << '\n';
+      vc = 0;  // global vertex offset
      for (p = 0; p < NRanks; p++)
      {
         if (p)
@@ -5202,8 +5394,11 @@ void ParMesh::PrintAsOne(std::ostream &os, const std::string &comments) const
         }
         for (i = 0; i < ne; )
         {
-            // processor number + 1 as bdr. attr. and bdr. geometry type
-            os << p+1 << ' ' << ints[i];
+            // Print bdr attribute + bdr geometry:
+            // * if print_shared != 0, use processor number + 1 as bdr attribute
+            // * otherwise, use the real bdr attribute
+            os << (print_shared ? p+1 : ints[i]) << ' ' << ints[i+1];
+            i++;
            k = Geometries.NumVerts[ints[i++]];
            // vertices
            for (j = 0; j < k; j++)
@@ -5212,7 +5407,7 @@ void ParMesh::PrintAsOne(std::ostream &os, const std::string &comments) const
            }
            os << '\n';
         }
-         vc += nv;
+         vc += nv;  // checked for overflow above, when printing the elements
      }
   }
   else
@@ -5227,6 +5422,7 @@ void ParMesh::PrintAsOne(std::ostream &os, const std::string &comments) const
   }

   // vertices / nodes
+   // checked for overflow above, when printing the elements:
   MPI_Reduce(const_cast<int*>(&NumOfVertices), &nv, 1, MPI_INT, MPI_SUM, 0,
              MyComm);
   if (MyRank == 0)
@@ -5253,8 +5449,8 @@ void ParMesh::PrintAsOne(std::ostream &os, const std::string &comments) const
            vert.SetSize(nv*spaceDim);
            if (nv)
            {
-               MPI_Recv(&vert[0], nv*spaceDim, MPITypeMap<real_t>::mpi_type, p, 449, MyComm,
-                        &status);
+               MPI_Recv(&vert[0], nv*spaceDim, MPITypeMap<real_t>::mpi_type, p,
+                        449, MyComm, &status);
            }
            for (i = 0; i < nv; i++)
            {
@@ -5281,8 +5477,8 @@ void ParMesh::PrintAsOne(std::ostream &os, const std::string &comments) const
         }
         if (NumOfVertices)
         {
-            MPI_Send(&vert[0], NumOfVertices*spaceDim, MPITypeMap<real_t>::mpi_type, 0, 449,
-                     MyComm);
+            MPI_Send(&vert[0], NumOfVertices*spaceDim,
+                     MPITypeMap<real_t>::mpi_type, 0, 449, MyComm);
         }
      }
   }
@@ -6879,7 +7075,9 @@ void ParMesh::Swap(ParMesh &other)
   // Nodes, NCMesh, and NURBSExtension are taken care of by Mesh::Swap
   mfem::Swap(pncmesh, other.pncmesh);

-   print_shared = other.print_shared;
+   // Keep Print() behavior consistent after move/swap operations.
+   mfem::Swap(print_shared, other.print_shared);
+   mfem::Swap(print_interfaces, other.print_interfaces);
 }

 void ParMesh::Destroy()
@@ -96,10 +96,12 @@ protected:
   mutable long glob_offset_sequence;
   void ComputeGlobalElementOffset() const;

-   // Enable Print() to add the parallel interface as boundary (typically used
-   // for visualization purposes)
+   // See SetPrintShared()
   bool print_shared = true;

+   // See SetPrintInterfaces()
+   bool print_interfaces = false;
+
   /// Create from a nonconforming mesh.
   ParMesh(const ParNCMesh &pncmesh);

@@ -303,6 +305,9 @@ protected:
   void GetSharedTriCommunicator(int ordering,
                                 GroupCommunicator& stria_comm) const;

+   // Optionally called by Print() and PrintAsOne() (so this needs to be const)
+   void FindInterface(Array<int> &interface) const;
+
   // Similar to Mesh::GetFacesTable()
   STable3D *GetSharedFacesTable();

@@ -720,13 +725,35 @@ public:
       begin with '#'. */
   void ParPrint(std::ostream &out, const std::string &comments = "") const;

-   // Enable Print() to add the parallel interface as boundary (typically used
-   // for visualization purposes)
+   /** @brief Enable Print() and PrintAsOne() to add the parallel interface as
+       boundary (typically used for visualization purposes).
+
+       In PrintAsOne(), this setting also controls what element and boundary
+       attributes are printed:
+       - if @a print == false, use the real element and boundary attributes,
+       - otherwise, processor rank + 1 is used for both, the element and
+         boundary attributes.
+
+      The ParMesh object itself is not modified, this only affects file output
+      for visualization.
+
+      The default value of this flag is true. */
   void SetPrintShared(bool print) { print_shared = print; }

+   /** @brief Enable Print() and PrintAsOne() to add material interfaces, i.e.
+       intefaces between different mesh element attributes, as boundary
+       (typically used for visualization purposes).
+
+       The ParMesh object itself is not modified, this only affects file output
+       for visualization.
+
+       The default value of this flag is false. */
+   void SetPrintInterfaces(bool print) { print_interfaces = print; }
+
   /** Print the part of the mesh in the calling processor using the mfem v1.0
-       format. Depending on SetPrintShared(), the parallel interface can be
-       added as boundary for visualization (true by default). If @a comments is
+       format. Depending on SetPrintShared() and SetPrintInterfaces(), the
+       parallel interface and/or material interfaces can be added as boundary
+       for visualization (true/false respectively by default). If @a comments is
       non-empty, it will be printed after the first line of the file, and each
       line should begin with '#'. */
   void Print(std::ostream &out = mfem::out,
@@ -748,12 +775,18 @@ public:
       as boundary (for visualization purposes) using Netgen/Truegrid format .*/
   void PrintXG(std::ostream &out = mfem::out) const override;

-   /** Write the mesh to the stream 'out' on Process 0 in a form suitable for
-       visualization: the mesh is written as a disjoint mesh and the shared
-       boundary is added to the actual boundary; both the element and boundary
-       attributes are set to the processor number. If @a comments is non-empty,
-       it will be printed after the first line of the file, and each line should
-       begin with '#'. */
+   /** @brief Write the mesh to the stream 'out' on Process 0 in a form suitable
+       for visualization.
+
+       The mesh is written as a disjoint mesh. If SetPrintShared() is enabled,
+       the shared boundary is added to the actual boundary and both the element
+       and boundary attributes are set to the processor number + 1.
+
+       If SetPrintInterfaces() is enabled, material interfaces, i.e. interfaces
+       between different mesh element attributes, are added as boundary as well.
+
+       If @a comments is non-empty, it will be printed after the first line of
+       the file, and each line should begin with '#'. */
   void PrintAsOne(std::ostream &out = mfem::out,
                   const std::string &comments = "") const;

@@ -98,7 +98,7 @@ point Tribol to them. The steps are as follows:
   AXOM_DIR = @MFEM_DIR@/../axom-repo/axom
   TRIBOL_DIR = @MFEM_DIR@/../tribol-repo/tribol
   TRIBOL_OPT = -I$(TRIBOL_DIR)/include -I$(AXOM_DIR)/include
-   TRIBOL_LIB = -L$(TRIBOL_DIR)/lib -ltribol -lredecomp -L$(AXOM_DIR)/lib \
+   TRIBOL_LIB = -L$(TRIBOL_DIR)/lib -ltribol -ltribol_shared -lredecomp -L$(AXOM_DIR)/lib \
   -laxom_mint -laxom_slam -laxom_slic -laxom_core
   ```
 3. [**Axom:**](https://github.com/LLNL/axom.git) Starting from the MFEM root
Author	SHA1	Message	Date
Ketan Mittal	eeaf33d17c	Merge branch 'master' of https://github.com/mfem/mfem into mesh-dist	2026-06-09 10:00:05 -07:00
Ketan Mittal	70cc72c462	initial commit	2026-06-09 09:59:57 -07:00
Tzanio Kolev	c6ea42a681	Merge pull request #5045 from mfem/stefanozampini/petsc-requires-hypre PETSc: Remove ad-hoc code to support PETSc without hypre	2026-06-06 16:52:11 -07:00
Tzanio Kolev	61bb7755a3	Merge pull request #5324 from mfem/multigrid-bcs-fix Bugfix for GeometricMultigrid with no essential BCs	2026-06-05 16:17:51 -07:00
Tzanio Kolev	a7fa61464c	Merge pull request #5267 from mfem/print-interfaces-dev Optional output of material interfaces in parallel	2026-06-05 16:13:18 -07:00
Tzanio Kolev	93329a7ba0	Merge pull request #5124 from yuyangdai/cudss-dev Add support for parallel NVIDIA's GPU-accelerated direct sparse solver cuDSS solver[cudss-dev]	2026-06-05 16:10:03 -07:00
yuyangdai	fb93e21a0c	Update INSTALL file to correct cuDSS library options	2026-06-05 09:16:00 +08:00
Stefano Zampini	19d60f95ec	PETSc: Remove ad-hoc code to support PETSc without hypre	2026-06-04 11:27:36 +01:00
Will Pazner	861d318162	Fix Doxygen comment in GeometricMultigrid	2026-06-03 15:48:39 -07:00
Tzanio Kolev	d29e914146	Merge branch 'master' into multigrid-bcs-fix	2026-06-02 18:20:52 -07:00
yuyangdai	b5e99d8d83	Fix CMake script to correctly locate cuDSS library and allocate device memory for csr values when set the cudssMatrix.	2026-06-02 13:31:25 +08:00
daiyuyangandAndrew Ho	91f58c293d	Update linalg/cudss.cpp Co-authored-by: Andrew Ho <ho37@llnl.gov>	2026-06-02 13:17:09 +08:00
Tzanio Kolev	6451e64637	Merge branch 'master' into cudss-dev	2026-06-01 15:55:15 -07:00
Tzanio Kolev	b7ac1963bb	Merge pull request #5340 from mfem/update-contact-miniapp-readme Update Tribol shared libraries	2026-06-01 15:53:50 -07:00
Veselin Dobrev	baefb786fb	Correction in CHANGELOG	2026-06-01 08:58:54 -07:00
Veselin Dobrev	99523b3a97	Address reviewer feedback: add/improve Doxygen documentation.	2026-06-01 08:43:12 -07:00
Tzanio Kolev	59d6820ff3	Merge branch 'master' into cudss-dev	2026-05-28 14:21:30 -07:00
Tzanio Kolev	9ec379f176	Merge branch 'master' into update-contact-miniapp-readme	2026-05-28 14:19:45 -07:00
Tzanio Kolev	72cc503eb9	Updated CHANGELOG	2026-05-28 12:53:25 -07:00
Veselin Dobrev	5aff935c98	Minor future-proofing suggested by Copilot. MFEM_ASSERT does not need to explicitly print the name of the function because that is already done automatically.	2026-05-27 10:18:49 -07:00
Tzanio KolevandCopilot Autofix powered by AI	d56298ba17	Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com>	2026-05-26 19:02:02 -07:00
Tzanio KolevandCopilot Autofix powered by AI	b9d67dec34	Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com>	2026-05-26 19:00:55 -07:00
Tzanio Kolev	5b664e393d	Merge branch 'master' into print-interfaces-dev	2026-05-26 18:43:14 -07:00
Veselin Dobrev	c582282084	Merge pull request #5345 from mfem/copilot-dev Updates in developer docs + instructions for the GitHub Copilot reviews	2026-05-26 18:27:49 -07:00
Tzanio Kolev	944ece4090	TODO item for ConduitDataCollection::Save()	2026-05-26 17:35:13 -07:00
Veselin DobrevandTzanio Kolev	f225d0c3ef	Apply suggestions from code review Remove FIXME comments -- no actions needed. Co-authored-by: Tzanio Kolev <tzanio@llnl.gov>	2026-05-26 16:58:36 -07:00
Tzanio Kolev	a9c98c2e3a	Merge branch 'master' into print-interfaces-dev	2026-05-26 10:35:15 -07:00
Tzanio KolevandCopilot Autofix powered by AI	ac8c2948e7	Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com>	2026-05-26 10:02:50 -07:00
Tzanio KolevandCopilot Autofix powered by AI	f75b4c10d2	Potential fix for pull request finding Co-authored-by: Copilot Autofix powered by AI <175728472+Copilot@users.noreply.github.com>	2026-05-26 10:02:10 -07:00
Veselin Dobrev	ecf6b0b44f	Address some reviewer suggestions and comments	2026-05-26 08:41:10 -07:00
Tzanio Kolev	5a52c676e1	Small changes in Copilot instructions	2026-05-25 15:49:19 -07:00
Tzanio Kolev	5be9693235	Initial Copilot instructions	2026-05-25 15:34:45 -07:00
Tzanio Kolev	5421abc4c4	Fixes and updates in CONTRIBUTING.md	2026-05-25 12:58:57 -07:00
Tzanio Kolev	f2d3fb45dc	Fixes and updates in INSTALL	2026-05-25 12:04:26 -07:00
John Camier	239c83a742	Merge branch 'master' into cudss-dev	2026-05-24 20:14:04 -07:00
maxpaik16	c315028a09	Merge branch 'master' into update-contact-miniapp-readme	2026-05-24 19:03:34 -07:00
yuyangdai	99b66f4f45	Update CHANGELOG	2026-05-25 09:20:19 +08:00
Tzanio Kolev	e7be50eb91	Merge pull request #5333 from nmnobre/gslib FindPointsGSLIB: use parallel-aware ProjectDiscCoefficient	2026-05-24 10:58:16 -07:00
Tzanio Kolev	2dd3915ad5	Merge branch 'master' into gslib	2026-05-23 10:15:21 -07:00
Tzanio Kolev	29b572819a	Merge branch 'master' into cudss-dev	2026-05-23 10:14:37 -07:00
Tzanio Kolev	0a0acfda66	Merge pull request #5329 from nmnobre/host Fix GetEssentialTrueDofsVar memory allocation	2026-05-22 20:15:51 -07:00
maxpaik16	7dbce44472	Update tribol libraries	2026-05-22 16:36:51 -07:00
maxpaik16	c67bd21219	Update tribol shared libraries	2026-05-22 16:35:29 -07:00
Tzanio Kolev	38f5b93520	Merge branch 'master' into print-interfaces-dev	2026-05-21 09:45:23 -07:00
Tzanio Kolev	a078dfd59e	Reviewer comments	2026-05-21 09:45:02 -07:00
Tzanio Kolev	66818f5525	Merge branch 'master' into gslib	2026-05-21 07:41:18 -07:00
Tzanio Kolev	b810a5e540	Merge branch 'master' into host	2026-05-21 07:41:14 -07:00
Tzanio Kolev	e63e421343	Merge branch 'master' into cudss-dev	2026-05-21 07:40:59 -07:00
Nuno Nobre	95acb1f85f	Remove unnecessary namespace qualification	2026-05-19 17:15:03 +01:00
Nuno Nobre	1c7261db07	Fence code using parallel objs w/ MPI-conditional directive	2026-05-18 23:53:35 +01:00
Nuno Nobre	e1b678664a	FindPointsGSLIB: use parallel-aware ProjectDiscCoefficient	2026-05-18 19:02:30 +01:00
Nuno NobreandAndrew Ho	8f06539b6e	Do not assume true_ess_dofs is empty nor host-allocated Co-authored-by: Andrew Ho <ho37@llnl.gov>	2026-05-16 00:32:25 +01:00
Nuno Nobre	30a8eb5ccf	Fix GetEssentialTrueDofsVar memory allocation	2026-05-14 23:58:41 +01:00
John Camier	078e59a33b	Merge branch 'master' into cudss-dev	2026-05-14 13:20:48 -07:00
Tzanio KolevandVeselin Dobrev	da4e1b5137	Update makefile Co-authored-by: Veselin Dobrev <v-dobrev@users.noreply.github.com>	2026-05-12 13:54:55 -07:00
Tzanio Kolev	b5a1660c2d	Merge branch 'master' into cudss-dev	2026-05-09 10:41:10 -07:00
Will Pazner	5ffc2ef502	Bugfix for GeometricMultigrid with no essential BCs	2026-05-07 18:04:46 -07:00
yuyangdai	59cef5f9e3	Add a link for communication layer library in cuDSS	2026-05-06 10:07:31 +08:00
John Camier	05be944a86	Merge branch 'master' into cudss-dev	2026-05-05 15:33:09 -07:00
John Camier	c103cfa84a	Merge branch 'master' into cudss-dev	2026-05-05 06:34:53 -07:00
John Camier	29caf08098	Merge branch 'master' into cudss-dev	2026-04-29 17:05:01 -07:00
daiyuyang	e217864f16	Merge branch 'master' into cudss-dev	2026-04-28 15:17:58 +08:00
yuyangdai	de8aacddce	Replace enum class with enum	2026-04-28 14:59:09 +08:00
Veselin Dobrev	f1174bfbf5	Added new methods in class ParFiniteElementSpace: HaveDofSigns and ApplyDofSigns. Used the new methods to fix bugs in: * the ParGridFunction constructor that reads input from a stream * the method ParGridFunction::SaveAsOne	2026-04-17 09:08:19 -07:00
Veselin Dobrev	65acd08d38	Small tweak in ParMesh::Print	2026-04-10 19:50:10 -07:00
Veselin Dobrev	0acc85d962	In ParMesh::PrintAsOne() fix the interface attribute to make it different from all real boundary attributes.	2026-04-10 17:41:26 -07:00
Veselin Dobrev	4e35641c28	Fix GCC warning	2026-04-10 12:57:48 -07:00
Veselin Dobrev	8594867ab6	Modify ParMesh::PrintAsOne() to use the settings of SetPrintShared() and SetPrintInterfaces(). Added some suggestions/questions as FIXME comments.	2026-04-10 12:21:02 -07:00
John Camier	6e0fbbd3bf	Merge branch 'master' into cudss-dev	2026-04-09 06:43:58 -07:00
yuyangdai	a2b34ab650	Add conditional compilation for cudss_solver in ex1p	2026-04-08 15:05:50 +08:00
yuyangdai	f478f687ff	Update cuDSS solver integration: - Use the full name of default communication library and threading library. - Check the `CUDSS_COMM_LIB` and `CUDSS_THREADING_LIB` in environment first. - Add the `cudss-solver` option in ex1p	2026-04-08 13:48:56 +08:00
daiyuyangandWill Pazner	0bf998510f	Update config/defaults.mk Co-authored-by: Will Pazner <11493037+pazner@users.noreply.github.com>	2026-04-08 09:45:40 +08:00
yuyangdai	d916299a49	Update INSTALL documentation for CUDSS requirements and backend support	2026-04-03 14:33:59 +08:00
yuyangdai	10840ac6b4	Refactor CuDSSSolver initialization.	2026-04-03 13:56:20 +08:00
yuyangdai	b437016f6e	Refactor cuDSS library path resolution in defaults.mk for improved handling of missing libraries	2026-04-02 15:01:00 +08:00
yuyangdai	5170bbe010	Fix CUDSS library checks and improve status output in CMake and makefile	2026-04-01 16:16:01 +08:00
yuyangdai	f8fcdf6a97	Update cuDSS configuration and library paths in CMake and makefiles; refactor cuDSS solver integration. - Add `MFEM_CUDSS_COMM_LIB` for OpenMPI communication library - Add 'MFEM_CUDSS_THREADING_LIB' for threading library - Add SetMatrixCuDSS() to set matrix values for cudss - Remove SetMatrixSortRow() - Remove unused options and simplify conditions in ex1.cpp and ex1p.cpp.	2026-04-01 13:49:02 +08:00
John Camier	d10c908b38	Merge branch 'master' into cudss-dev	2026-03-29 17:51:58 -07:00
John Camier	89f85cee21	Merge branch 'master' into cudss-dev	2026-03-26 13:49:44 -07:00
camierjs	00da00b93a	Add Copyright headers	2026-03-24 08:42:34 -07:00
camierjs	4eca673111	Cosmetic trailing spaces	2026-03-24 08:33:06 -07:00
John Camier	4507a02249	Merge branch 'master' into cudss-dev	2026-03-24 08:01:00 -07:00
Tzanio Kolev	61c7187a86	Review comments	2026-03-12 11:17:14 -07:00
Tzanio KolevandCopilot	4fbefc6987	Update mesh/pmesh.cpp Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com>	2026-03-11 12:20:13 -07:00
Tzanio Kolev	8dd75d2548	typo	2026-03-11 11:32:06 -07:00
Tzanio Kolev	b0c2ec505f	Added optional output of material interfaces in parallel	2026-03-11 11:20:59 -07:00
yuyangdai	a3e229ef82	Make SetMatrix methods private and remove Init method - SetMatrix methods are called inside SetOperator - Init method has been removed as it is unnecessary	2026-02-09 16:40:39 +08:00
John Camier	00b6dcdd37	Merge branch 'master' into cudss-dev	2026-02-07 13:21:22 -08:00
Tzanio Kolev	19f444489f	Merge branch 'master' into cudss-dev	2026-02-05 11:03:38 -08:00
daiyuyang	10c637837c	Add dependency check for CUDSS in MFEMConfig.cmake.in	2026-01-22 16:56:17 +08:00
daiyuyangandAndrew Ho	665ba30f65	Update linalg/cudss.hpp Co-authored-by: Andrew Ho <ho37@llnl.gov>	2026-01-22 09:38:05 +08:00
daiyuyangandAndrew Ho	d6de2c1a1d	Update examples/ex1p.cpp Co-authored-by: Andrew Ho <ho37@llnl.gov>	2026-01-22 09:35:32 +08:00
daiyuyangandAndrew Ho	c4b4cfdc2f	Update examples/ex1p.cpp Co-authored-by: Andrew Ho <ho37@llnl.gov>	2026-01-22 09:35:18 +08:00
daiyuyangandAndrew Ho	3435475ae0	Update examples/ex1.cpp Co-authored-by: Andrew Ho <ho37@llnl.gov>	2026-01-22 09:34:32 +08:00
daiyuyangandAndrew Ho	e43b58fa02	Update config/cmake/modules/FindCUDSS.cmake Co-authored-by: Andrew Ho <ho37@llnl.gov>	2026-01-22 09:32:27 +08:00
John Camier	27c8412439	Merge branch 'master' into cudss-dev	2026-01-21 08:23:31 -08:00
daiyuyang	5aaef22cc0	Add conditional support for cuDSS solver in ex1.cpp	2026-01-19 15:16:31 +08:00
daiyuyang	f9df36a6de	Update parameter names in documentation for clarity in cudss.hpp	2026-01-19 14:21:19 +08:00
daiyuyang	b3a75a9295	Add cuDSS solver option in ex1.cpp	2026-01-19 13:42:00 +08:00
daiyuyang	b8c90d24f1	Make MPI optional and add OpenMP support in CuDSSSolver - Make the MPI optional: the CuDSSSolver now does not need MPI necessary - Add the OpenMP supports by cudssSetThreadingLayer() API	2026-01-19 13:40:56 +08:00
daiyuyang	a8eba594d2	Update linalg/cudss.hpp and linalg/cudss.cpp; Revert general/device.cpp - Move ` CUDA_REAL_T` to the cpp file. - Remove `CuDSSHandle` singleton and replace it with a static cudssHandle_t variable. - Delete `CuDSSHandle::Init()` from ex1p	2025-12-17 09:52:36 +08:00
yuyangdai	52291cadbf	Revert code docs; add constructor comment; update the `.gitignore`	2025-12-15 13:09:45 +08:00
yuyangdai	1865b430b0	Revert `makefile` and `example/CMakeLists.txt`; delete `examples/cudss`	2025-12-15 11:01:42 +08:00
yuyangdai	d8734b4b18	Update the examples/ex1p.cpp - Add the global cudss handle before using the cudss solver.	2025-12-15 10:27:47 +08:00
yuyangdai	eae1fa217b	Add `CuDSSHandle` singleton and update `CuDSSSolver` - Remove unused variable `myid`. - Move `MFEM_CUDSS_CHECK` and `mfem_cudss_error` into linalg/cudss.cpp. - Disable move copy constructor and move assignment. - Introduce `CuDSSHandle` singleton to manage cuDSS handle lifetime. - Rename `InitHandle` to `Init`	2025-12-15 10:27:22 +08:00
gengyan.zgy	7dc8b0d9fe	Fix CuDSSSolver::ArrayMult() for single RHS	2025-12-11 14:32:20 +08:00
yuyangdai	3f236b406f	Update examples/ex1p.cpp - Add cudss solver into the ex1p.cpp	2025-12-11 13:55:16 +08:00
yuyangdai	47dcdd8159	Update config/cmake/modules/FindCUDSS.cmake and examples/cudss/CMakeLists.txt - Add a newline at end of file	2025-12-11 13:53:49 +08:00
yuyangdai	d93d36d6d3	Update linalg/cudss.cpp and linalg/cudss.hpp - Define a cuDSS error check macro, MFEM_CUDSS_CHECK(x); - Rename the method from InitRhsSol to SetNumRHS; - Utilize CuMemAlloc/ CuMemcpyDtoD instead of cudaMalloc/cudaMemcpy; - Update MPI_Comm usage.	2025-12-11 13:52:57 +08:00
yuyangdai	6bcad07a5c	Update the `n_global` variable name for global number of rows	2025-12-05 13:50:36 +08:00
yuyangdai	d56493c2c5	Refactored enum classes MatType and MatViewType, updated variable names, and update some comments.	2025-12-05 11:01:29 +08:00
daiyuyang	ec67fe536e	Merge branch 'master' into cudss-dev	2025-12-02 14:04:17 +08:00
daiyuyang	eff793d6fa	Merge branch 'master' into cudss-dev	2025-11-28 13:53:59 +08:00
yuyangdai	e5e7280a17	Ignore the solution files of cudss sample run	2025-11-24 15:18:53 +08:00
yuyangdai	773dc57712	Using the `make style` to format the class CuDSSSolver and example/cudss/ex1p codes	2025-11-24 15:16:48 +08:00
yuyangdai	98a043b195	Update the doc files for CuDSSSolver	2025-11-24 15:15:22 +08:00
yuyangdai	c2a3c83099	Update the makefiles for CuDSSSolver	2025-11-24 15:14:51 +08:00
daiyuyang	3ea90aff27	Add cudss solver	2025-11-24 15:10:45 +08:00
gengyan.zgy	23a362d6ad	Update cmake files for cuDSS	2025-11-24 15:10:30 +08:00