documentation

* Update GitHub CI to handle Enzyme/LLVM changes in Homebrew.

* GitHub CI: try to allow the Enzyme job to fail

* GitHub CI: fix a typo in last commit

* GitHub CI: another try to allow the Enzyme job to fail.

* GitHub CI: another try to allow the Enzyme job to fail.

* GitHub CI: another try to allow the Enzyme job to fail.

* GitHub CI: another try to allow the Enzyme job to fail.

* GitHub CI: try without fail-fast == false.

* GitHub CI: try to fix the LLVM link issue in the Enzyme job.

.github/workflows/builds-and-tests.yml

@@ -132,12 +132,14 @@ jobs:
             hypre-target: int32
             precision: fp64
             enzyme: true
-            config-opts: MFEM_USE_ENZYME=YES ENZYME_DIR=$(brew --prefix enzyme)
+            config-opts: MFEM_USE_ENZYME=YES ENZYME_DIR=$(brew --prefix enzyme) LDFLAGS=-L$LLVM_PREFIX/lib/c++
 
     name: ${{ matrix.os }}-${{ matrix.build-system }}-${{ matrix.target }}-${{ matrix.mpi }}-${{ matrix.hypre-target }}-${{ matrix.precision }}${{ matrix.enzyme && '-enzyme' || '' }}
 
     runs-on: ${{ matrix.os }}
 
+    continue-on-error: ${{ matrix.enzyme && true || false }}
+
     steps:
       # Fix 'No space left on device' errors for Ubuntu builds.
       - name: Run Actions Cleaner
@@ -292,10 +294,12 @@ jobs:
         run: |
           export HOMEBREW_NO_INSTALL_CLEANUP=1
           brew update
-          brew install llvm@20 enzyme
-          echo "LLVM_PREFIX=$(brew --prefix llvm@20)" >> $GITHUB_ENV
-          echo "OMPI_CC=$(brew --prefix llvm@20)/bin/clang" >> $GITHUB_ENV
-          echo "OMPI_CXX=$(brew --prefix llvm@20)/bin/clang++" >> $GITHUB_ENV
+          brew install enzyme
+          ENZYME_LLVM=$(brew info enzyme | sed -n 's/^Required:.*\(llvm[^ ]*\).*/\1/p')
+          LLVM_PREFIX=$(brew --prefix $ENZYME_LLVM)
+          echo "LLVM_PREFIX=$LLVM_PREFIX" >> $GITHUB_ENV
+          echo "OMPI_CC=$LLVM_PREFIX/bin/clang" >> $GITHUB_ENV
+          echo "OMPI_CXX=$LLVM_PREFIX/bin/clang++" >> $GITHUB_ENV
 
       # MFEM build and test
       - name: build

fem/bilinearform.cpp

@@ -1275,22 +1275,6 @@ void BilinearForm::Update(FiniteElementSpace *nfes)
    height = width = fes->GetVSize();
 
    if (ext) { ext->Update(); }
-   for (int k = 0; k < domain_integs.Size(); ++k)
-   {
-      domain_integs[k]->Update();
-   }
-   for (int k = 0; k < boundary_integs.Size(); ++k)
-   {
-      boundary_integs[k]->Update();
-   }
-   for (int k = 0; k < interior_face_integs.Size(); ++k)
-   {
-      interior_face_integs[k]->Update();
-   }
-   for (int k = 0; k < boundary_integs.Size(); ++k)
-   {
-      boundary_face_integs[k]->Update();
-   }
 }
 
 void BilinearForm::SetDiagonalPolicy(DiagonalPolicy policy)
@@ -2353,31 +2337,6 @@ void MixedBilinearForm::Update()
    height = test_fes->GetVSize();
    width = trial_fes->GetVSize();
    if (ext) { ext->Update(); }
-
-   for (int k = 0; k < domain_integs.Size(); ++k)
-   {
-      domain_integs[k]->Update();
-   }
-   for (int k = 0; k < boundary_integs.Size(); ++k)
-   {
-      boundary_integs[k]->Update();
-   }
-   for (int k = 0; k < interior_face_integs.Size(); ++k)
-   {
-      interior_face_integs[k]->Update();
-   }
-   for (int k = 0; k < boundary_integs.Size(); ++k)
-   {
-      boundary_face_integs[k]->Update();
-   }
-   for (int k = 0; k < trace_face_integs.Size(); ++k)
-   {
-      trace_face_integs[k]->Update();
-   }
-   for (int k = 0; k < boundary_trace_face_integs.Size(); ++k)
-   {
-      boundary_trace_face_integs[k]->Update();
-   }
 }
 
 MixedBilinearForm::~MixedBilinearForm()

fem/bilininteg.cpp

@@ -21,11 +21,6 @@ using namespace std;
 namespace mfem
 {
 
-void BilinearFormIntegrator::Update()
-{
-   // default no-op
-}
-
 void BilinearFormIntegrator::AssemblePA(const FiniteElementSpace&)
 {
    MFEM_ABORT("BilinearFormIntegrator::AssemblePA(fes)\n"
@@ -3465,12 +3460,6 @@ real_t ElasticityIntegrator::ComputeFluxEnergy(const FiniteElement &fluxelem,
    return energy;
 }
 
-void DGTraceIntegrator::Update()
-{
-   qspace[0].reset();
-   qspace[1].reset();
-}
-
 void DGTraceIntegrator::AssembleFaceMatrix(const FiniteElement &el1,
                                            const FiniteElement &el2,
                                            FaceElementTransformations &Trans,

fem/bilininteg.hpp

@@ -23,8 +23,6 @@
 
 namespace mfem
 {
-class QuadratureSpace;
-class FaceQuadratureSpace;
 
 /// Abstract base class BilinearFormIntegrator
 class BilinearFormIntegrator : public NonlinearFormIntegrator
@@ -46,10 +44,6 @@ public:
    // make sense for the action of the nonlinear operator (but they all make
    // sense for its Jacobian).
 
-   /// Signal this integrator that something about either the trial or test space has changed.
-   virtual void Update();
-
-
    /// Method defining partial assembly.
    /** The result of the partial assembly is stored internally so that it can be
        used later in the methods AddMultPA() and AddMultTransposePA(). */
@@ -3317,7 +3311,6 @@ protected:
    VectorCoefficient *u;
    real_t alpha, beta;
    // PA extension
-   std::unique_ptr<FaceQuadratureSpace> qspace[2];
    Vector pa_data;
    const DofToQuad *maps;             ///< Not owned
    const FaceGeometricFactors *geom;  ///< Not owned
@@ -3340,8 +3333,6 @@ public:
                      real_t a, real_t b)
    { rho = &rho_; u = &u_; alpha = a; beta = b; }
 
-   void Update() override;
-
    using BilinearFormIntegrator::AssembleFaceMatrix;
    void AssembleFaceMatrix(const FiniteElement &el1,
                            const FiniteElement &el2,

fem/fespace.hpp

@@ -922,6 +922,9 @@ public:
    { return mesh->GetBdrElementType(i); }
 
    /// Returns ElementTransformation for the @a i-th element.
+   /// @note The returned pointer references an object owned by the associated
+   /// @a Mesh that will be modified by other calls to `GetElementTransformation`.
+   /// As such, this pointer should @b not be deleted by the caller.
    ElementTransformation *GetElementTransformation(int i) const
    { return mesh->GetElementTransformation(i); }
 

fem/integ/bilininteg_dgtrace_pa.cpp

@@ -139,6 +139,8 @@ void DGTraceIntegrator::SetupPA(const FiniteElementSpace &fes, FaceType type)
    const MemoryType mt = (pa_mt == MemoryType::DEFAULT) ?
                          Device::GetDeviceMemoryType() : pa_mt;
 
+   nf = fes.GetNFbyType(type);
+   if (nf==0) { return; }
    // Assumes tensor-product elements
    Mesh *mesh = fes.GetMesh();
    const FiniteElement &el = *fes.GetTypicalTraceElement();
@@ -146,17 +148,6 @@ void DGTraceIntegrator::SetupPA(const FiniteElementSpace &fes, FaceType type)
                                IntRule:
                                &GetRule(el.GetGeomType(), el.GetOrder(),
                                         *mesh->GetTypicalElementTransformation());
-
-   if (!qspace[static_cast<int>(type)])
-   {
-      qspace[static_cast<int>(type)].reset(
-         new FaceQuadratureSpace(*mesh, *ir, type));
-   }
-
-   FaceQuadratureSpace& qs = *qspace[static_cast<int>(type)];
-   nf = qs.GetNumFaces();
-   if (nf==0) { return; }
-
    const int symmDims = 4;
    nq = ir->GetNPoints();
    dim = mesh->Dimension();
@@ -168,6 +159,8 @@ void DGTraceIntegrator::SetupPA(const FiniteElementSpace &fes, FaceType type)
    dofs1D = maps->ndof;
    quad1D = maps->nqpt;
    pa_data.SetSize(symmDims * nq * nf, Device::GetMemoryType());
+
+   FaceQuadratureSpace qs(*mesh, *ir, type);
    CoefficientVector vel(*u, qs, CoefficientStorage::COMPRESSED);
 
    CoefficientVector r(qs, CoefficientStorage::COMPRESSED);

fem/qspace.cpp

@@ -17,9 +17,8 @@ namespace mfem
 {
 
 QuadratureSpaceBase::QuadratureSpaceBase(Mesh &mesh_, Geometry::Type geom,
-                                         const IntegrationRule &ir,
-                                         QSpaceStorage storage)
-   : mesh(mesh_), order(ir.GetOrder()), storage(storage)
+                                         const IntegrationRule &ir)
+   : mesh(mesh_), order(ir.GetOrder())
 {
    for (int g = 0; g < Geometry::NumGeom; g++)
    {
@@ -97,10 +96,10 @@ void QuadratureSpaceBase::Integrate(VectorCoefficient &coeff,
 
 void QuadratureSpace::ConstructOffsets()
 {
-   const int num_elem = ne;
+   const int num_elem = mesh.GetNE();
+   ne = num_elem;
 
-   if (storage == QSpaceStorage::COMPRESSED &&
-       mesh.GetNumGeometries(mesh.Dimension()) == 1)
+   if (mesh.GetNumGeometries(mesh.Dimension()) == 1)
    {
       Array<Geometry::Type> geoms;
       mesh.GetGeometries(mesh.Dimension(), geoms);
@@ -125,9 +124,14 @@ void QuadratureSpace::ConstructOffsets()
    }
 }
 
-QuadratureSpace::QuadratureSpace(Mesh *mesh_, std::istream &in,
-                                 QSpaceStorage storage)
-   : QuadratureSpaceBase(*mesh_, 0, storage)
+void QuadratureSpace::Construct()
+{
+   ConstructIntRules(mesh.Dimension());
+   ConstructOffsets();
+}
+
+QuadratureSpace::QuadratureSpace(Mesh *mesh_, std::istream &in)
+   : QuadratureSpaceBase(*mesh_)
 {
    const char *msg = "invalid input stream";
    std::string ident;
@@ -146,24 +150,15 @@ QuadratureSpace::QuadratureSpace(Mesh *mesh_, std::istream &in,
       return;
    }
 
-   ne = mesh.GetNE();
-   ConstructIntRules(mesh.Dimension());
+   Construct();
 }
 
-QuadratureSpace::QuadratureSpace(Mesh *mesh_, int order_, QSpaceStorage storage)
-   : QuadratureSpaceBase(*mesh_, order_, storage)
-{
-   ne = mesh.GetNE();
-   ConstructIntRules(mesh.Dimension());
-}
-
-QuadratureSpace::QuadratureSpace(Mesh &mesh_, const IntegrationRule &ir,
-                                 QSpaceStorage storage)
-   : QuadratureSpaceBase(mesh_, mesh_.GetTypicalElementGeometry(), ir, storage)
+QuadratureSpace::QuadratureSpace(Mesh &mesh_, const IntegrationRule &ir)
+   : QuadratureSpaceBase(mesh_, mesh_.GetTypicalElementGeometry(), ir)
 {
    MFEM_VERIFY(mesh.GetNumGeometries(mesh.Dimension()) <= 1,
                "Constructor not valid for mixed meshes");
-   ne = mesh.GetNE();
+   ConstructOffsets();
 }
 
 void QuadratureSpace::Save(std::ostream &os) const
@@ -185,53 +180,55 @@ const Vector &QuadratureSpace::GetGeometricFactorWeights() const
 }
 
 FaceQuadratureSpace::FaceQuadratureSpace(Mesh &mesh_, int order_,
-                                         FaceType face_type_,
-                                         QSpaceStorage storage)
-   : QuadratureSpaceBase(mesh_, order_, storage), face_type(face_type_),
-     face_indices(mesh.GetFaceIndices(face_type_)),
-     face_indices_inv(mesh.GetInvFaceIndices(face_type_))
+                                         FaceType face_type_)
+   : QuadratureSpaceBase(mesh_, order_),
+     face_type(face_type_),
+     num_faces(mesh.GetNFbyType(face_type))
 {
-   ne = face_indices.Size();
-   ConstructIntRules(mesh.Dimension() - 1);
+   Construct();
 }
 
 FaceQuadratureSpace::FaceQuadratureSpace(Mesh &mesh_, const IntegrationRule &ir,
-                                         FaceType face_type_,
-                                         QSpaceStorage storage)
-   : QuadratureSpaceBase(mesh_, mesh_.GetTypicalFaceGeometry(), ir, storage),
-     face_type(face_type_), face_indices(mesh.GetFaceIndices(face_type_)),
-     face_indices_inv(mesh.GetInvFaceIndices(face_type_))
+                                         FaceType face_type_)
+   : QuadratureSpaceBase(mesh_, mesh_.GetTypicalFaceGeometry(), ir),
+     face_type(face_type_),
+     num_faces(mesh.GetNFbyType(face_type))
 {
    MFEM_VERIFY(mesh.GetNumGeometries(mesh.Dimension() - 1) <= 1,
                "Constructor not valid for mixed meshes");
-   ne = face_indices.Size();
+   ConstructOffsets();
 }
 
 void FaceQuadratureSpace::ConstructOffsets()
 {
-   if (storage == QSpaceStorage::COMPRESSED &&
-       mesh.GetNumGeometries(mesh.Dimension() - 1) == 1)
+   face_indices.SetSize(num_faces);
+   offsets.SetSize(num_faces + 1);
+   ne = num_faces;
+   int offset = 0;
+   int f_idx = 0;
+   for (int i = 0; i < mesh.GetNumFacesWithGhost(); i++)
    {
-      Array<Geometry::Type> geoms;
-      mesh.GetGeometries(mesh.Dimension() - 1, geoms);
-      offsets.SetSize(1);
-      offsets.HostWrite();
-      offsets[0] = int_rule[geoms[0]]->GetNPoints();
-      size = ne * offsets[0];
-   }
-   else
-   {
-      offsets.SetSize(face_indices.Size() + 1);
-      int offset = 0;
-      for (int i = 0; i < mesh.GetNFbyType(face_type); ++i)
+      const Mesh::FaceInformation face = mesh.GetFaceInformation(i);
+      if (face.IsNonconformingCoarse() || !face.IsOfFaceType(face_type))
       {
-         offsets[i] = offset;
-         Geometry::Type geom = mesh.GetFaceGeometry(face_indices[i]);
-         MFEM_ASSERT(int_rule[geom] != NULL, "Missing integration rule");
-         offset += int_rule[geom]->GetNPoints();
+         continue;
       }
-      offsets[face_indices.Size()] = size = offset;
+      face_indices[f_idx] = i;
+      face_indices_inv[i] = f_idx;
+      offsets[f_idx] = offset;
+      Geometry::Type geom = mesh.GetFaceGeometry(i);
+      MFEM_ASSERT(int_rule[geom] != NULL, "Missing integration rule");
+      offset += int_rule[geom]->GetNPoints();
+
+      f_idx++;
    }
+   offsets[num_faces] = size = offset;
+}
+
+void FaceQuadratureSpace::Construct()
+{
+   ConstructIntRules(mesh.Dimension() - 1);
+   ConstructOffsets();
 }
 
 int FaceQuadratureSpace::GetPermutedIndex(int idx, int iq) const

fem/qspace.hpp

@@ -19,49 +19,39 @@
 namespace mfem
 {
 
-enum class QSpaceStorage
-{
-   FULL,
-   COMPRESSED
-};
-
 /// Abstract base class for QuadratureSpace and FaceQuadratureSpace.
 /** This class represents the storage layout for QuadratureFunction%s, that may
     be defined either on mesh elements or mesh faces. */
 class QuadratureSpaceBase
 {
 protected:
+   friend class QuadratureFunction; // Uses the offsets.
+
    Mesh &mesh; ///< The underlying mesh.
    int order; ///< The order of integration rule.
-   int size = -1; ///< Total number of quadrature points. -1 indicates
-   ///< offsets/size not computed yet.
+   int size; ///< Total number of quadrature points.
    int ne; ///< Actual number of entities
    mutable Vector weights; ///< Integration weights.
    mutable long nodes_sequence = 0; ///< Nodes counter for cache invalidation.
 
-   QSpaceStorage storage;
-
    /// @brief Entity quadrature point offset array.
    ///
    /// Supports a constant compression scheme for meshes which have a single
    /// geometry type. When compressed, will have a single value. The true offset
    /// can be computed as i * offsets[0], where i is the entity index. Otherwise
-   /// has size num_entities + 1. Lazily constructed.
+   /// has size num_entities + 1.
    ///
    Array<int> offsets;
    /// The quadrature rules used for each geometry type.
    const IntegrationRule *int_rule[Geometry::NumGeom];
 
    /// Protected constructor. Used by derived classes.
-   QuadratureSpaceBase(Mesh &mesh_, int order_ = 0,
-                       QSpaceStorage storage = QSpaceStorage::COMPRESSED)
-      : mesh(mesh_), order(order_), storage(storage)
-   {}
+   QuadratureSpaceBase(Mesh &mesh_, int order_ = 0)
+      : mesh(mesh_), order(order_) { }
 
    /// Protected constructor. Used by derived classes.
    QuadratureSpaceBase(Mesh &mesh_, Geometry::Type geom,
-                       const IntegrationRule &ir,
-                       QSpaceStorage storage = QSpaceStorage::COMPRESSED);
+                       const IntegrationRule &ir);
 
    /// Fill the @ref int_rule array for each geometry type using @ref order.
    void ConstructIntRules(int dim);
@@ -72,21 +62,13 @@ protected:
    /// Compute the integration weights.
    void ConstructWeights() const;
 
-   virtual void ConstructOffsets() = 0;
-
 public:
-   QSpaceStorage StorageType() const { return storage; }
-
    /// @brief Gets the offset for a given entity @a idx.
    ///
    /// The quadrature point values for entity i are stored in the indices
    /// between Offset(i) and Offset(i+1)
    int Offset(int idx) const
    {
-      if (size < 0)
-      {
-         const_cast<QuadratureSpaceBase *>(this)->ConstructOffsets();
-      }
       return (offsets.Size() == 1) ? (idx * offsets[0]) : offsets[idx];
    }
 
@@ -97,24 +79,10 @@ public:
    /// can be computed as i * offsets[0], where i is the entity index. Otherwise
    /// has size num_entities + 1.
    ///
-   const Array<int> &Offsets() const
-   {
-      if (size < 0)
-      {
-         const_cast<QuadratureSpaceBase *>(this)->ConstructOffsets();
-      }
-      return offsets;
-   }
+   const Array<int> &Offsets() const { return offsets; }
 
    /// Return the total number of quadrature points.
-   int GetSize() const
-   {
-      if (size < 0)
-      {
-         const_cast<QuadratureSpaceBase *>(this)->ConstructOffsets();
-      }
-      return size;
-   }
+   int GetSize() const { return size; }
 
    /// Return the order of the quadrature rule(s) used by all elements.
    int GetOrder() const { return order; }
@@ -174,20 +142,19 @@ class QuadratureSpace : public QuadratureSpaceBase
 {
 protected:
    const Vector &GetGeometricFactorWeights() const override;
-   void ConstructOffsets() override;
+   void ConstructOffsets();
+   void Construct();
 public:
    /// Create a QuadratureSpace based on the global rules from #IntRules.
-   QuadratureSpace(Mesh *mesh_, int order_,
-                   QSpaceStorage storage = QSpaceStorage::COMPRESSED);
+   QuadratureSpace(Mesh *mesh_, int order_)
+      : QuadratureSpaceBase(*mesh_, order_) { Construct(); }
 
    /// @brief Create a QuadratureSpace with an IntegrationRule, valid only when
    /// the mesh has one element type.
-   QuadratureSpace(Mesh &mesh_, const IntegrationRule &ir,
-                   QSpaceStorage storage = QSpaceStorage::COMPRESSED);
+   QuadratureSpace(Mesh &mesh_, const IntegrationRule &ir);
 
    /// Read a QuadratureSpace from the stream @a in.
-   QuadratureSpace(Mesh *mesh_, std::istream &in,
-                   QSpaceStorage storage = QSpaceStorage::COMPRESSED);
+   QuadratureSpace(Mesh *mesh_, std::istream &in);
 
    /// Returns number of elements in the mesh.
    inline int GetNE() const { return mesh.GetNE(); }
@@ -224,29 +191,29 @@ public:
 class FaceQuadratureSpace : public QuadratureSpaceBase
 {
    FaceType face_type; ///< Is the space defined on interior or boundary faces?
+   const int num_faces; ///< Number of faces.
 
    /// Map from boundary or interior face indices to mesh face indices.
-   const Array<int> &face_indices;
+   Array<int> face_indices;
 
    /// Inverse of the map @a face_indices.
-   const std::unordered_map<int,int> &face_indices_inv;
+   std::unordered_map<int,int> face_indices_inv;
 
    const Vector &GetGeometricFactorWeights() const override;
-   void ConstructOffsets() override;
+   void ConstructOffsets();
+   void Construct();
 
 public:
    /// Create a FaceQuadratureSpace based on the global rules from #IntRules.
-   FaceQuadratureSpace(Mesh &mesh_, int order_, FaceType face_type_,
-                       QSpaceStorage storage = QSpaceStorage::COMPRESSED);
+   FaceQuadratureSpace(Mesh &mesh_, int order_, FaceType face_type_);
 
    /// @brief Create a FaceQuadratureSpace with an IntegrationRule, valid only
    /// when the mesh has one type of face geometry.
    FaceQuadratureSpace(Mesh &mesh_, const IntegrationRule &ir,
-                       FaceType face_type_,
-                       QSpaceStorage storage = QSpaceStorage::COMPRESSED);
+                       FaceType face_type_);
 
    /// Returns number of faces in the mesh.
-   inline int GetNumFaces() const { return face_indices.Size(); }
+   inline int GetNumFaces() const { return num_faces; }
 
    /// Returns the face type (boundary or interior).
    FaceType GetFaceType() const { return face_type; }

general/array.hpp

@@ -211,9 +211,6 @@ public:
    /// Delete the first entry with value == 'el'.
    inline void DeleteFirst(const T &el);
 
-   /// Delete entries at @a indices, and resize.
-   inline void DeleteAt(const Array<int> &indices);
-
    /// Delete the whole array.
    inline void DeleteAll();
 
@@ -938,30 +935,6 @@ inline void Array<T>::DeleteFirst(const T &el)
    }
 }
 
-template <class T>
-inline void Array<T>::DeleteAt(const Array<int> &indices)
-{
-   // Make a copy of the indices, sorted.
-   Array<int> sorted_indices(indices);
-   sorted_indices.Sort();
-
-   int rm_count = 0;
-   for (int i = 0; i < size; i++)
-   {
-      if (rm_count < sorted_indices.Size() && i == sorted_indices[rm_count])
-      {
-         rm_count++;
-      }
-      else
-      {
-         data[i-rm_count] = data[i]; // shift data rm_count
-      }
-   }
-
-   // Resize to remove tail
-   size -= rm_count;
-}
-
 template <class T>
 inline void Array<T>::DeleteAll()
 {

general/scan.hpp

@@ -1,176 +0,0 @@
-// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
-// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
-// LICENSE and NOTICE for details. LLNL-CODE-806117.
-//
-// This file is part of the MFEM library. For more information and source code
-// availability visit https://mfem.org.
-//
-// MFEM is free software; you can redistribute it and/or modify it under the
-// terms of the BSD-3 license. We welcome feedback and contributions, see file
-// CONTRIBUTING.md for details.
-
-#ifndef MFEM_SCAN_HPP
-#define MFEM_SCAN_HPP
-
-#ifdef MFEM_USE_CUDA
-#include <cub/device/device_scan.cuh>
-#define MFEM_CUB_NAMESPACE cub
-#elif MFEM_USE_HIP
-#include <hipcub/device/device_scan.hpp>
-#define MFEM_CUB_NAMESPACE hipcub
-#endif
-
-#include <functional>
-#include <numeric>
-
-namespace mfem
-{
-/// Equivalent to InclusiveScan(use_dev, d_in, d_out, num_items, workspace,
-/// std::plus<>{})
-template <class InputIt, class OutputIt>
-void InclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
-                   Array<char> &workspace)
-{
-   // forward to InclusiveSum for potentially faster kernels
-#if defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP)
-   if (use_dev && mfem::Device::Allows(Backend::CUDA_MASK | Backend::HIP_MASK))
-   {
-      size_t bytes = workspace.Size();
-      if (bytes)
-      {
-         auto err = MFEM_CUB_NAMESPACE::DeviceScan::InclusiveSum(
-                       workspace.Write(), bytes, d_in, d_out, num_items);
-#if defined(MFEM_USE_CUDA)
-         if (err == cudaSuccess)
-         {
-            return;
-         }
-#elif defined(MFEM_USE_HIP)
-         if (err == hipSuccess)
-         {
-            return;
-         }
-#endif
-      }
-      // try allocating a larger buffer
-      bytes = 0;
-      MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::InclusiveSum(
-                        nullptr, bytes, d_in, d_out, num_items));
-      workspace.SetSize(bytes);
-      MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::InclusiveSum(
-                        workspace.Write(), bytes, d_in, d_out, num_items));
-      return;
-   }
-#endif
-   std::inclusive_scan(d_in, d_in + num_items, d_out);
-}
-
-/// Performs an inclusive scan of [d_in, d_in+num_items) -> [d_out,
-/// d_out+num_items). This call is potentially asynchronous on the device.
-/// @a d_in input start.
-/// @a d_out output start. Can perform in-place scans with d_out = d_in
-/// @a workspace temporary workspace used for device scans. TODO: replace with
-/// internal temporary workspace once that's added to the memory manager.
-/// @a scan_op binary scan functor. Must be associative. If only weakly
-/// associative (i.e. floating point addition) results are not deterministic. On
-/// device this must also be commutative.
-template <class InputIt, class OutputIt, class ScanOp>
-void InclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
-                   Array<char> &workspace, ScanOp scan_op)
-{
-#if defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP)
-   if (use_dev && mfem::Device::Allows(Backend::CUDA_MASK | Backend::HIP_MASK))
-   {
-      size_t bytes = workspace.Size();
-      if (bytes)
-      {
-         auto err = MFEM_CUB_NAMESPACE::DeviceScan::InclusiveScan(
-                       workspace.Write(), bytes, d_in, d_out, scan_op, num_items);
-#if defined(MFEM_USE_CUDA)
-         if (err == cudaSuccess)
-         {
-            return;
-         }
-#elif defined(MFEM_USE_HIP)
-         if (err == hipSuccess)
-         {
-            return;
-         }
-#endif
-      }
-      // try allocating a larger buffer
-      bytes = 0;
-      MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::InclusiveScan(
-                        nullptr, bytes, d_in, d_out, scan_op, num_items));
-      workspace.SetSize(bytes);
-      MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::InclusiveScan(
-                        workspace.Write(), bytes, d_in, d_out, scan_op, num_items));
-      return;
-   }
-#endif
-   std::inclusive_scan(d_in, d_in + num_items, d_out, scan_op);
-}
-
-/// Performs an exclusive scan of [d_in, d_in+num_items) -> [d_out,
-/// d_out+num_items). This call is potentially asynchronous on the device.
-/// @a d_in input start.
-/// @a d_out output start. Can perform in-place scans with d_out = d_in
-/// @a workspace temporary workspace used for device scans. TODO: replace with
-/// internal temporary workspace once that's added to the memory manager.
-/// @a scan_op binary scan functor. Must be associative. If only weakly
-/// associative (i.e. floating point addition) results are not deterministic. On
-/// device this must also be commutative.
-template <class InputIt, class OutputIt, class T, class ScanOp>
-void ExclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
-                   T init_value, Array<char> &workspace, ScanOp scan_op)
-{
-#if defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP)
-   if (use_dev && mfem::Device::Allows(Backend::CUDA_MASK | Backend::HIP_MASK))
-   {
-      size_t bytes = workspace.Size();
-      if (bytes)
-      {
-         auto err = MFEM_CUB_NAMESPACE::DeviceScan::ExclusiveScan(
-                       workspace.Write(), bytes, d_in, d_out, scan_op, init_value,
-                       num_items);
-#if defined(MFEM_USE_CUDA)
-         if (err == cudaSuccess)
-         {
-            return;
-         }
-#elif defined(MFEM_USE_HIP)
-         if (err == hipSuccess)
-         {
-            return;
-         }
-#endif
-      }
-      // try allocating a larger buffer
-      bytes = 0;
-      MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::ExclusiveScan(
-                        nullptr, bytes, d_in, d_out, scan_op, init_value, num_items));
-      workspace.SetSize(bytes);
-      MFEM_GPU_CHECK(MFEM_CUB_NAMESPACE::DeviceScan::ExclusiveScan(
-                        workspace.Write(), bytes, d_in, d_out, scan_op, init_value,
-                        num_items));
-      return;
-   }
-#endif
-   std::exclusive_scan(d_in, d_in + num_items, d_out, init_value, scan_op);
-}
-
-/// Equivalent to ExclusiveScan(use_dev, d_in, d_out, num_items, init_value,
-/// workspace, std::plus<>{})
-template <class InputIt, class OutputIt, class T>
-void ExclusiveScan(bool use_dev, InputIt d_in, OutputIt d_out, size_t num_items,
-                   T init_value, Array<char> &workspace)
-{
-   ExclusiveScan(use_dev, d_in, d_out, num_items, init_value, workspace,
-                 std::plus<> {});
-}
-
-} // namespace mfem
-
-#undef MFEM_CUB_NAMESPACE
-
-#endif

linalg/CMakeLists.txt

@@ -23,6 +23,7 @@ list(APPEND SRCS
   complex_operator.cpp
   constraints.cpp
   densemat.cpp
+  functional.cpp
   symmat.cpp
   handle.cpp
   matrix.cpp
@@ -49,6 +50,7 @@ list(APPEND HDRS
   complex_operator.hpp
   constraints.hpp
   densemat.hpp
+  functional.hpp
   dinvariants.hpp
   symmat.hpp
   dtensor.hpp

linalg/functional.cpp

@@ -0,0 +1,68 @@
+#include "functional.hpp"
+
+namespace mfem
+{
+QuadraticFunctional::QuadraticFunctional(const Operator *A_,
+                                         const Vector *b_, const real_t beta_, const real_t c_)
+   : Functional(A_ ? A_->Width() : 0)
+   , A(A_), beta(beta_), b(b_), c(c_)
+   , aux(A_ ? A_->Width() : 0)
+{
+   MFEM_VERIFY(A_ != nullptr,
+               "QuadraticFunctional: A must not be nullptr. "
+               << "Use QuadraticFunctional() constructor to create an empty Quadratic functional.");
+   MFEM_VERIFY(A_->Width() == A_->Height(),
+               "QuadraticFunctional: A must be a square operator.");
+   MFEM_VERIFY(b_ == nullptr || A_->Width() == b_->Size(),
+               "QuadraticFunctional: A and b must have compatible sizes");
+}
+
+void QuadraticFunctional::SetOperator(const Operator &A_)
+{
+   MFEM_VERIFY(A_.Width() == A_.Height(),
+               "QuadraticFunctional: A must be a square operator.");
+   A = &A_;
+   width = A_.Width();
+   aux.SetSize(width);
+}
+
+void QuadraticFunctional::SetVector(const Vector &b_, const real_t beta_)
+{
+   MFEM_VERIFY(A != nullptr && A->Width() == b_.Size(),
+               "QuadraticFunctional: A and b must have compatible sizes.");
+   b = &b_;
+   beta = beta_;
+}
+void QuadraticFunctional::SetConstant(real_t c_) { c = c_; }
+
+void QuadraticFunctional::Mult(const Vector &x, Vector &y) const
+{
+   MFEM_ASSERT(A != nullptr, "QuadraticFunctional: A must not be nullptr");
+   A->Mult(x, aux); // aux = A(x)
+   if (b != nullptr) { aux.Add(2.0*beta, *b); } // aux = A(x) + 2*beta*b
+   real_t result = 0.0;
+#ifdef MFEM_USE_MPI
+   if (IsParallel())
+   {
+      result = InnerProduct(GetComm(), x, aux); // result = <A(x) + 2*beta*b, x>
+   }
+   else
+   {
+      result = InnerProduct(x, aux); // result = <A(x) + 2*beta*b, x>
+   }
+#else
+   result = InnerProduct(x, aux); // result  = <A(x) + 2*beta*b, x>
+#endif
+   y.SetSize(1);
+   y[0] = result*0.5 + c; // y = 0.5 * <A(x) + 2*beta*b, x> + c
+}
+
+void QuadraticFunctional::EvalGradient(const Vector &x,
+                                       Vector &y) const
+{
+   y.SetSize(A->Width());
+   A->Mult(x, y);
+   if (b) { y.Add(beta,*b); }
+}
+
+} // namespace mfem

linalg/functional.hpp

@@ -0,0 +1,744 @@
+#ifndef MFEM_FUNCTIONAL_HPP
+#define MFEM_FUNCTIONAL_HPP
+
+#include "../config/config.hpp"
+
+#ifdef MFEM_USE_MPI
+#include "../general/communication.hpp"
+#endif
+
+#include "operator.hpp"
+#include "blockvector.hpp"
+#include "solvers.hpp"
+#include <cxxabi.h>
+#include <vector>
+
+namespace mfem
+{
+/// @brief A base class for functionals F:R^n->R
+///
+/// This class provides an interface for evaluating
+/// $ F:R^n->R, \nabla F:R^n->R^n $, and $ \nabla^2 F:R^n x R^n ->R^n $.
+/// F.Mult(x, y) evaluates the functional at a point x, and stores the result in y[0]
+/// F.GetGradient() returns an operator that evaluates the gradient
+/// F.GetGradient().GetGradient(x) returns an Hessian action operator.
+///
+/// The usual Operator::GetGradient(const Vector &x) method for this method
+/// is deprecated as $ \nabla F $ only takes a single argument x.
+/// It is redundant to use F.GetGradient(x).Mult(x, y) to evaluate the gradient.
+/// Instead, use F.GetGradient().Mult(x, y) to evaluate the gradient at x.
+///
+/// The gradient and Hessian can be defined in two ways:
+/// 1. If the gradient is available as a sperate operator,
+///    then override Functional::GetGradient().
+///    In this case, Functional::HessianMult() will not be called.
+///
+/// 2. Otherwise, override Functional::EvalGradient() and Functional::HessianMult()
+///    to evaluate the gradient and Hessian action, respectively.
+///    The helper classes, GradientOperator and HessianActionOperator,
+///    will call these methods to evaluate the gradient and Hessian action.
+///    If Hessian is a seperate operator, then you can override
+///    The GradientOperator::GetGradient(x) will call Functional::GetHessian(x)
+///
+class Functional : public Operator
+{
+   Operator * riesz_map = nullptr; ///< Riesz map operator, if available
+public:
+   /// @brief Create a Functional with optional gradient and hessian
+   /// @param n number of variables
+   Functional(int n=0)
+      : Operator(1, n)
+      , grad_operator(*this)
+      , hessian_action_operator(*this)
+   { }
+
+#ifdef MFEM_USE_MPI
+   Functional(MPI_Comm comm, int n=0)
+      : Functional(n)
+   { SetComm(comm); }
+
+   void SetComm(MPI_Comm comm_)
+   {
+      parallel = comm_ != MPI_COMM_NULL;
+      comm = comm_;
+   }
+   MPI_Comm GetComm() const { return comm; }
+   bool IsParallel() const { return parallel; }
+#else
+   constexpr bool IsParallel() const { return false; }
+#endif
+
+   void SetRieszMap(Operator &op) { riesz_map = &op; }
+   /// @brief return the GradientOperator that evaluates the gradient
+   ///        input x is not used. Use GetGradient().Mult(x,y) to evaluate the gradient
+   ///        we recommend using GetGradient() instead of GetGradient(x)
+   /// Deprecated. See Functional::GetGradient()
+   MFEM_DEPRECATED
+   Operator &GetGradient(const Vector &dummy) const override final { return GetGradient(); }
+
+   /// @brief Return the GradientOperator that wraps Functional::EvalGradient() for Mult().
+   /// @note If the functional has a corresponding standalone gradient operator,
+   /// override this method to return the gradient operator.
+   virtual Operator &GetGradient() const { return grad_operator; }
+
+   /// @brief Evaluate the functional at a point x that will be called by GradientOperator::Mult()
+   /// @note This method is not meant to be called directly. See, GradientOperator
+   virtual void EvalGradient(const Vector &x, Vector &y) const
+   {
+      MFEM_ABORT("Functional::EvalGradient() not implemented");
+   }
+   /// @brief Evaluate the Hessian action at a point x and direction d
+   /// that will be called by Functional::GetGradient().GetHessian(x).Mult(d,y)
+   /// @note This method is not meant to be called directly. See, HessianActionOperator
+   virtual void HessianMult(const Vector &x, const Vector &d, Vector &y) const
+   {
+      MFEM_ABORT("Functional::HessianMult() not implemented.");
+   }
+   /// @brief Return the HessianActionOperator at evaluation point x
+   ///        that wraps Functional::HessianMult() for Mult().
+   ///        See, HessianActionOperator and Functional::HessianMult().
+   ///
+   /// @note If the Hessian is available as a seperate operator, override this method.
+   ///
+   /// @warning If GetGradient() is overridden, this method will not be used.
+   virtual Operator &GetHessian(const Vector &x) const
+   {
+      hessian_action_operator.SetX(x);
+      return hessian_action_operator;
+   }
+private:
+#ifdef MFEM_USE_MPI
+   bool parallel=false;
+#else
+   const static bool parallel=false;
+#endif
+#ifdef MFEM_USE_MPI
+   MPI_Comm comm;
+#endif
+   /// @brief A helper class to return an operator that evaluates the gradient
+   ///        using Functional::EvalGradient() method.
+   class GradientOperator : public Operator
+   {
+   private: const Functional &f; mutable Vector der;
+   public:
+      GradientOperator(const Functional &f) : Operator(f.Width()), f(f) {}
+      /// @brief Evaluate the gradient of Functional at a point x
+      void Mult(const Vector &x, Vector &y) const override final
+      {
+         if (f.riesz_map)
+         {
+            der.SetSize(f.Width());
+            f.EvalGradient(x, der);
+            f.riesz_map->Mult(der, y);
+         }
+         else
+         {
+            f.EvalGradient(x, y);
+         }
+      }
+      /// @brief Evaluate the Hessian of Functional at a point x
+      Operator &GetGradient(const Vector &x) const override final { return f.GetHessian(x); }
+   };
+   friend class GradientOperator;
+
+   /// @brief A helper class to return an operator that applies the Hessian action
+   ///        using Functional::HessianMult() method.
+   class HessianActionOperator : public Operator
+   {
+   private:
+      const Functional &f;
+      const Vector *x;
+   public:
+      HessianActionOperator(const Functional &f) : Operator(f.Width()), f(f) {}
+      void SetX(const Vector &new_x) { x = &new_x; }
+      void Mult(const Vector &d, Vector &y) const override { f.HessianMult(*x, d, y); }
+   };
+   friend class HessianActionOperator;
+
+   mutable GradientOperator grad_operator;
+   mutable HessianActionOperator hessian_action_operator;
+};
+
+/// @brief Stacked functioanl operator, [f1, ..., fk] where fi:R^n->R are functionals
+/*
+   Typical usage of this class is to provide a single operator for multiple constraints.
+   For example, consider a minimization problem with k constraints,
+   min f0(u) s.t. fi(u)=0, i=1,...,k.
+   The Lagrangian functional is
+   $ L(u, lambda) = F0(u) + sum_i lambda_i * fi(u) $
+   The first-order optimality conditions are
+   $ \nabla f0(u) + \sum_i lambda_i * grad fi(u) = 0 $
+   $ fi(u) = 0 $
+   where lambda_i are the Lagrange multipliers.
+   The StackedFunctional class can be used to represent the list of constraints fi(u).
+
+   StackedFunctional::Mult(u, y) will evaluate each functional y[i]=fi(u)
+
+   StackedFunctional::GetGradient(u) represents an operator, column-stacked gradient
+   That is, [grad f0(u), ..., grad fk(u)] in R^{n x k}
+   If you want to extract the gradient as a matrix, use
+   StackedFunctional::GetGradientMatrix(const Vector &x, DenseMatrix &grad)
+   As functionals are not assumed to return a sparse vector, the gradient is dense.
+
+   StackedFunctional::GetGradient(u).Mult(lambda, y) contract the gradients with the Lagrange multipliers
+   y = sum lambda_i * grad fi(u)
+   StackedFunctional::GetGradient(u).MultTranspose(d, y) return the directional derivative for each k
+   y[i] = <grad fi(u), d>
+
+   StackedFunctional::GetHessian(u, lambda).Mult(d, y) will return the contracted Hessian action
+   $ y = \sum_i \lambda_i * H_{fi}(u, d) $
+*/
+/// @warning Functionals should be all serial or all parallel.
+///
+class StackedFunctional : public Operator
+{
+public:
+   StackedFunctional(int n=0)
+      : Operator(0, n)
+      , funcs(0)
+      , grad_helper_op(*this)
+      , hessian_helper_op(*this)
+   {}
+   StackedFunctional(Functional &f)
+      : Operator(0, f.Width())
+      , grad_helper_op(*this)
+      , hessian_helper_op(*this)
+   { AddFunctional(f); }
+
+   StackedFunctional(const std::vector<Functional*> &funcs)
+      : Operator((int)funcs.size(), funcs[0]->Width())
+      , grad_helper_op(*this)
+      , hessian_helper_op(*this)
+   { for (auto &f : funcs) { AddFunctional(*f); } }
+
+   void AddFunctional(Functional &f)
+   {
+#ifdef MFEM_USE_MPI
+      if (funcs.empty()) { if (f.IsParallel()) { SetComm(f.GetComm()); } }
+#endif
+      MFEM_VERIFY(f.Width() == Width(),
+                  "StackedFunctional::AddFunctional: Functional width does not match with the operator.");
+      MFEM_VERIFY(parallel == f.IsParallel(),
+                  "StackedFunctional::AddFunctional: Parallelism mismatch.");
+      funcs.push_back(&f);
+      height++;
+   }
+
+   void Mult(const Vector &x, Vector &y) const override
+   {
+      y.SetSize(Height());
+      Vector yview;
+      for (int i=0; i<Height(); i++)
+      {
+         yview.MakeRef(y, i, 1);
+         funcs[i]->Mult(x, yview);
+      }
+   }
+
+   Operator &GetGradient(const Vector &x) const override
+   {
+      grad_helper_op.SetX(x);
+      return grad_helper_op;
+   }
+
+   void GetGradientMatrix(const Vector &x, DenseMatrix &grads) const
+   {
+      grads.SetSize(Width(), Height());
+      Vector grad;
+      for (int i=0; i<Height(); i++)
+      {
+         grads.GetColumnReference(i, grad);
+         funcs[i]->GetGradient().Mult(x, grad);
+      }
+   }
+   Functional &GetFunctional(int i) const
+   {
+      MFEM_VERIFY(i >= 0 && i < Height(),
+                  "StackedFunctional::GetFunctional: Index out of bounds.");
+      return *funcs[i];
+   }
+
+   Operator &GetHessian(const Vector &x, const Vector &lambda) const
+   {
+      hessian_helper_op.SetX(x, lambda);
+      return hessian_helper_op;
+   }
+
+   bool parallel;
+   bool IsParallel() const { return parallel; }
+#ifdef MFEM_USE_MPI
+   void SetComm(MPI_Comm comm_)
+   {
+      parallel = comm != MPI_COMM_NULL;
+      comm = comm_;
+   }
+   MPI_Comm GetComm() const { return comm; }
+#endif
+
+protected:
+#ifdef MFEM_USE_MPI
+   MPI_Comm comm;
+#endif
+   std::vector<Functional*> funcs;
+
+   class GradientOperator : public Operator
+   {
+   public:
+      GradientOperator(const StackedFunctional &op)
+         : Operator(op.Width(), op.Height())
+         , op(op)
+         , tmp_grad(op.Width())
+      {}
+      void SetX(const Vector &x) const { x_curr = &x; }
+      Operator &GetGradient(const Vector &lambda) const override
+      {
+         op.hessian_helper_op.SetX(*x_curr, lambda);
+         return op.hessian_helper_op;
+      }
+
+      void Mult(const Vector &lambda, Vector &y) const override
+      {
+         y.SetSize(op.Width());
+         y = 0.0;
+         for (int i=0; i<op.Height(); i++)
+         {
+            op.funcs[i]->GetGradient().Mult(*x_curr, tmp_grad);
+            y.Add(lambda[i], tmp_grad);
+         }
+      }
+      void MultTranspose(const Vector &x, Vector &y) const override
+      {
+         y.SetSize(op.Height());
+         for (int i=0; i<op.Height(); i++)
+         {
+            op.funcs[i]->GetGradient().Mult(x, tmp_grad);
+            y[i] = InnerProduct(tmp_grad, *x_curr);
+         }
+#ifdef MFEM_USE_MPI
+         if (op.IsParallel())
+         {
+            MPI_Allreduce(MPI_IN_PLACE, y.GetData(), op.Height(),
+                          MPITypeMap<real_t>::mpi_type, MPI_SUM,
+                          op.GetComm());
+         }
+#endif
+      }
+
+   private:
+      const StackedFunctional &op;
+      mutable const Vector *x_curr;
+      mutable Vector tmp_grad;
+
+   };
+   class HessianActionOperator : public Operator
+   {
+   public:
+      HessianActionOperator(const StackedFunctional &op)
+         : Operator(op.Width()), op(op)
+      {}
+      void SetX(const Vector &x, const Vector &lambda) const { x_curr = &x; lambda_curr = &lambda; }
+
+      void Mult(const Vector &d, Vector &y) const override
+      {
+         y.SetSize(op.Width());
+         y = 0.0;
+         for (int i=0; i<op.Height(); i++)
+         {
+            op.funcs[i]->GetGradient().GetGradient(*x_curr).Mult(d, tmp_hessian);
+            y.Add((*lambda_curr)[i], tmp_hessian);
+         }
+      }
+   private:
+      const StackedFunctional &op;
+      mutable Vector tmp_hessian;
+      mutable const Vector *x_curr;
+      mutable const Vector *lambda_curr;
+   };
+   friend class GradientOperator;
+   friend class HessianActionOperator;
+   mutable GradientOperator grad_helper_op;
+   mutable HessianActionOperator hessian_helper_op;
+private:
+};
+
+class ConstrainedOptimizationProblem : public Functional
+{
+public:
+   ConstrainedOptimizationProblem(Functional &objective_,
+                                  Operator *eq_constraints_=nullptr,
+                                  Operator *ineq_constraints_=nullptr)
+      : Functional(objective_.Width())
+      , objective(objective_)
+      , eq_constraints(eq_constraints_)
+      , ineq_constraints(ineq_constraints_)
+   {
+      // Check Size
+      MFEM_VERIFY((eq_constraints == nullptr ||
+                   eq_constraints->Width() == objective.Width()),
+                  "ConstrainedFunctional: Equality constraints width does not match with the objective.");
+      MFEM_VERIFY((ineq_constraints == nullptr ||
+                   ineq_constraints->Width() == objective.Width()),
+                  "ConstrainedFunctional: Inequality constraints width does not match with the objective.");
+#ifdef MFEM_USE_MPI
+      if (objective.IsParallel()) { SetComm(objective.GetComm()); }
+#endif
+   }
+
+   Functional &GetObjective() { return objective; }
+   const Functional &GetObjective() const { return objective; }
+   Operator *GetEqualityConstraints() { return eq_constraints; }
+   const Operator *GetEqualityConstraints() const { return eq_constraints; }
+   Operator *GetInequalityConstraints() { return ineq_constraints; }
+   const Operator *GetInequalityConstraints() const { return ineq_constraints; }
+protected:
+   Functional &objective;
+   Operator *eq_constraints;
+   Operator *ineq_constraints;
+};
+
+/// @brief A Lagrangian functional for
+/// min F(u)
+/// subject to C(u) = 0
+/// That is, L(u, lambda) = F(u) + <lambda, C(u)>
+///
+/// We assume that $ F:R^n -> R $ is a functional,
+///                $ C:R^n -> R^k $ is an equality constraint operator,
+/// C should return a residual. That is,
+/// C(u) = c, then C.Mult(u, y) should return y[i] = C_i(u) - c_i.
+///
+/// C.GetGradient(u):R^k -> R^n that takes lambda and returns the contracted gradient at x
+/// C.GetGradient(u).Mult(lambda, y) returns y = sum lambda_i * grad C_i(u)
+///
+/// C's gradient should support MultTranspose method
+/// That is, C.GetGradient(u).MultTranspose(d, y) returns y[i] = <grad C_i(u), d>
+///
+class LagrangianFunctional : public ConstrainedOptimizationProblem
+{
+private:
+   mutable Vector eq_residual;
+public:
+   LagrangianFunctional(Functional &objective,
+                        Operator &eq_constraints)
+      : ConstrainedOptimizationProblem(objective, &eq_constraints)
+      , eq_residual(eq_constraints.Height())
+   {
+      width = objective.Width() + eq_constraints.Height();
+
+      offsets.SetSize(3);
+      offsets[0] = 0;
+      offsets[1] = objective.Width();
+      offsets[2] = eq_constraints.Height();
+      offsets.PartialSum();
+   }
+
+   void Mult(const Vector &x, Vector &y) const override
+   {
+      const BlockVector input_block(const_cast<Vector&>(x), offsets);
+      const Vector &u = input_block.GetBlock(0);
+      const Vector &lambda = input_block.GetBlock(1);
+
+      y.SetSize(1);
+      y[0] = 0.0;
+      objective.Mult(u, y);
+      eq_constraints->Mult(u, eq_residual);
+      y[0] += InnerProduct(lambda, eq_residual);
+   }
+
+   void EvalGradient(const Vector &x, Vector &y) const override
+   {
+      const BlockVector input_block(const_cast<Vector&>(x), offsets);
+      const Vector &u = input_block.GetBlock(0);
+      const Vector &lambda = input_block.GetBlock(1);
+
+      y.SetSize(Width());
+      BlockVector output_block(y, offsets);
+      Vector &opt_residual = output_block.GetBlock(0);
+      eq_residual = output_block.GetBlock(1);
+      y = 0.0;
+      // grad F(u) + \sum_i lambda_i grad C_i(u)
+      objective.GetGradient().Mult(u, opt_residual);
+      eq_constraints->GetGradient(u).AddMult(lambda, opt_residual);
+      // grad C_i(u)^T
+      eq_constraints->GetGradient(u).MultTranspose(u, eq_residual);
+   }
+
+   /// @brief Evaluate the Hessian action at a point x=[u, lambda]
+   /// and direction d=[v, mu]
+   /// $ [H_F(u,d) + \sum_i lambda_i H_{C_i}(u, d) + <\nabla C_i(u), mu> $
+   void HessianMult(const Vector &x, const Vector &d, Vector &y) const override
+   {
+      const BlockVector input_block(const_cast<Vector&>(x), offsets);
+      const Vector &u = input_block.GetBlock(0);
+      const Vector &lambda = input_block.GetBlock(1);
+
+      const BlockVector direction_block(const_cast<Vector&>(x), offsets);
+      const Vector &v = direction_block.GetBlock(0);
+      const Vector &mu = direction_block.GetBlock(1);
+
+      y.SetSize(Width());
+      BlockVector output_block(y, offsets);
+      Vector &opt_H = output_block.GetBlock(0); // Optimality Hessian
+      Vector &eq_H = output_block.GetBlock(1); // Equality Hessian
+
+      // H_F(u,d) + \sum_i lambda_i H_{C_i}(u, d) + mu_i grad C_i(u)
+      objective.GetGradient().GetGradient(u).Mult(v, opt_H);
+      eq_constraints->GetGradient(u).GetGradient(lambda).AddMult(v, opt_H);
+      eq_constraints->GetGradient(u).Mult(mu, eq_H);
+      // <grad C_i(u), d>
+      eq_constraints->GetGradient(u).MultTranspose(d, eq_H);
+   }
+
+protected:
+   Array<int> offsets; // offsets for [x, lambda, mu]
+};
+
+/// @brief An augmented Lagrangian functional of the form
+/// F(u) + 0.5 mu * ||C(u)||^2 + <lambda, C(u)>
+/// where F is the objective functional,
+/// C is the equality constraint operator,
+/// lambda is the Lagrange multiplier vector (initialized to zero),
+/// mu is the penalty parameter (defaults to 1.0)
+///
+/// Currently, only equality constraints are supported.
+///
+/// AugLagrangianFunctional::Update() will update the penalty and Lagrange multiplier vectors
+/// By default, lambda <- lambda + mu * C(u)
+///             mu <- mu (no update)
+class AugLagrangianFunctional : public ConstrainedOptimizationProblem
+{
+public:
+   AugLagrangianFunctional(Functional &objective_,
+                           Operator &eq_constraints_)
+      : ConstrainedOptimizationProblem(objective_, &eq_constraints_)
+      , lambda(eq_constraints_.Height())
+      , mu(1.0)
+      , eq_residual(eq_constraints_.Height())
+      , eq_dir(eq_constraints_.Height())
+   {
+      lambda = 0.0;
+   }
+
+   void SetLambda(const Vector &lambda_)
+   {
+      MFEM_VERIFY(lambda_.Size() == eq_constraints->Height(),
+                  "AugLagrangianFunctional: Lambda size does not match with the equality constraints.");
+      lambda = lambda_;
+   }
+
+   void SetPenalty(real_t mu_)
+   {
+      MFEM_VERIFY(mu_ >= 0.0,
+                  "AugLagrangianFunctional: Penalty parameter mu must be non-negative.");
+      mu = mu_;
+   }
+
+   virtual void Update(const Vector &x)
+   {
+      // Update the Lagrange multipliers
+      eq_constraints->AddMult(x, lambda, mu);
+      // Update the penalty parameter
+      // Do nothing
+   }
+   const Vector &GetLambda() const { return lambda; }
+   real_t GetPenalty() const { return mu; }
+
+
+   void Mult(const Vector &x, Vector &y) const override
+   {
+      y.SetSize(1);
+      objective.Mult(x, y);
+      eq_constraints->Mult(x, eq_residual);
+      y[0] += lambda*eq_residual;
+      y[0] += 0.5 * mu * (eq_residual*eq_residual);
+   }
+
+   void EvalGradient(const Vector &x, Vector &y) const override
+   {
+      y.SetSize(Width());
+      // grad F(x) + \sum_i (lambda_i + mu * C_i(x)) grad C_i(x)
+      Vector curr_lambda = lambda; // store lambda + mu * C(x)
+      objective.GetGradient().Mult(x, y);
+      eq_constraints->Mult(x, eq_residual);
+      curr_lambda.Add(mu, eq_residual);
+      eq_constraints->GetGradient(x).AddMult(curr_lambda, y);
+   }
+
+   /// @brief Evaluate the Hessian action at a point x=[u, lambda]
+   /// and direction d=[v, mu]
+   /// $ H_F(u,d) + \sum_i \lambda_i H_{C_i}(u, d) + <\nabla C_i(u), mu> $
+   void HessianMult(const Vector &x, const Vector &d, Vector &y) const override
+   {
+      // H_F(u,d) + \sum_i lambda_i H_{C_i}(u, d) + mu_i grad C_i(u) <grad C_i(u), d>
+      objective.GetGradient().GetGradient(x).Mult(d, y);
+
+      Vector curr_lambda = lambda;
+      eq_constraints->Mult(x, eq_residual);
+      curr_lambda.Add(mu, eq_residual);
+
+      eq_constraints->GetGradient(x).GetGradient(curr_lambda).AddMult(d, y);
+      // eq_dir = <grad C_i(u), d>
+      eq_constraints->GetGradient(x).MultTranspose(d, eq_dir);
+      // mu_i <grad C_i(u), eq_dir>
+      eq_constraints->GetGradient(x).AddMult(eq_dir, y, mu);
+   }
+
+protected:
+   Vector lambda;
+   real_t mu;
+   mutable Vector eq_residual; // residual of the equality constraints, R^k
+   // directional derivative of the equality constraints, R^k
+   mutable Vector eq_dir;
+};
+
+/// @brief Quadratic functional of the form
+///        f(u) = 0.5 * <A u, u> + beta<b, u> + c
+///        where A is a square (possibly nonlinear) operator,
+///        beta is a scalar (defaults to 1.0, not used when b is nullptr),
+///        b is a vector (independent of u, optional),
+///        c is a constant (independent of u, optional).
+///        GetHessian() returns the operator A.
+///
+class QuadraticFunctional : public Functional
+{
+public:
+   QuadraticFunctional()
+      : Functional(0)
+      , A(nullptr), b(nullptr), c(0.0)
+   {}
+   QuadraticFunctional(const Operator *A_, const Vector *b_=nullptr,
+                       const real_t beta_=1.0,
+                       const real_t c_=0.0);
+#ifdef MFEM_USE_MPI
+   QuadraticFunctional(MPI_Comm comm_)
+      : QuadraticFunctional()
+   { SetComm(comm_); }
+
+   QuadraticFunctional(MPI_Comm comm_, const Operator *A_,
+                       const Vector *b_=nullptr,
+                       const real_t beta_=1.0, const real_t c_=0.0)
+      : QuadraticFunctional(A_, b_, beta_, c_)
+   { SetComm(comm_); }
+#endif
+
+   void SetOperator(const Operator &A_);
+   void SetVector(const Vector &b_, const real_t beta_=1.0);
+   void SetConstant(real_t c_);
+
+   void Mult(const Vector &x, Vector &y) const override;
+   void EvalGradient(const Vector &x, Vector &y) const override;
+protected:
+   const Operator *A;
+   real_t beta;
+   const Vector *b;
+   real_t c;
+   mutable Vector aux;
+protected:
+   /// @brief return the underlying Operator A
+   /// @warning Modifying the returned operator leads to undefined behavior.
+   Operator& GetHessian(const Vector &dummy) const override
+   {
+      return const_cast<Operator&>(*A);
+   }
+};
+
+class Optimizer : public IterativeSolver
+{
+public:
+   Optimizer() : IterativeSolver(), f(nullptr) { }
+#ifdef MFEM_USE_MPI
+   Optimizer(MPI_Comm comm) : IterativeSolver(comm), f(nullptr) { }
+#endif
+   // @brief Set the subproblem functional operator
+   // @param op the functional operator
+   // @note The functional will be stored in subproblem, and oper will be set to the gradient of the functional.
+   void SetOperator(const Functional &f_)
+   {
+      f = &f_;
+      IterativeSolver::SetOperator(f_.GetGradient());
+   }
+   virtual void SetLinearSolver(Solver &prec) { IterativeSolver::SetPreconditioner(prec); }
+
+   /// @brief This will abort. Should be called only with a Functional operator.
+   void SetOperator(const Operator &op) override
+   {
+      MFEM_ABORT("OptSolver::SetOperator() should not be called directly. Use SetFunctional() instead.");
+   }
+protected:
+   const Functional * f;
+};
+
+class NewtonOptimizer : public Optimizer
+{
+private:
+   real_t step_size = 1.0; // default step size
+public:
+   NewtonOptimizer() : Optimizer() { }
+#ifdef MFEM_USE_MPI
+   NewtonOptimizer(MPI_Comm comm) : Optimizer(comm) { }
+#endif
+   void SetStepSize(real_t step_size_) { step_size = step_size_; }
+
+   void Mult(const Vector &x, Vector &y) const override
+   {
+      dx.SetSize(x.Size());
+      y.SetSize(x.Size());
+      y = x;
+      MFEM_ASSERT(f != nullptr,
+                  "NewtonOptimizer::Mult() called without a functional operator.");
+      MFEM_ASSERT(prec != nullptr,
+                  "NewtonOptimizer::Mult() called without a linear solver.");
+      for (int i=0; i<max_iter; i++)
+      {
+         oper->Mult(y, grad);
+         Operator &hess = oper->GetGradient(y);
+         prec->SetOperator(hess);
+         prec->Mult(grad, dx);
+         y.Add(-step_size, dx);
+         if (Dot(dx, dx) < abs_tol*abs_tol)
+         {
+            break;
+         }
+      }
+   }
+private:
+   mutable Vector grad;
+   mutable Vector dx;
+};
+
+class GradientDescentOptimizer : public Optimizer
+{
+private:
+   real_t step_size = 1.0; // default step size
+public:
+   GradientDescentOptimizer() : Optimizer() { }
+#ifdef MFEM_USE_MPI
+   GradientDescentOptimizer(MPI_Comm comm) : Optimizer(comm) { }
+#endif
+   void SetStepSize(real_t step_size_) { step_size = step_size_; }
+
+   void Mult(const Vector &x, Vector &y) const override
+   {
+      grad.SetSize(x.Size());
+      y.SetSize(x.Size());
+      y = x;
+      MFEM_ASSERT(f != nullptr,
+                  "NewtonOptimizer::Mult() called without a functional operator.");
+      MFEM_ASSERT(prec != nullptr,
+                  "NewtonOptimizer::Mult() called without a linear solver.");
+      for (int i=0; i<max_iter; i++)
+      {
+         oper->Mult(y, grad);
+         y.Add(-step_size, grad);
+         if (Dot(grad, grad) < abs_tol*abs_tol)
+         {
+            break;
+         }
+      }
+   }
+private:
+   mutable Vector grad;
+};
+
+
+} // namespace mfem
+#endif // MFEM_FUNCTIONAL_HPP

linalg/vector.cpp

@@ -14,7 +14,6 @@
 #include "../general/forall.hpp"
 #include "../general/reducers.hpp"
 #include "../general/hash.hpp"
-#include "../general/scan.hpp"
 #include "vector.hpp"
 
 #ifdef MFEM_USE_OPENMP
@@ -1253,44 +1252,4 @@ real_t Vector::Sum() const
    return res;
 }
 
-void Vector::DeleteAt(const Array<int> &indices)
-{
-   const bool use_dev = UseDevice();
-
-   Array<int> flag(size);
-   const auto d_flag = flag.Write(use_dev);
-   mfem::forall_switch(use_dev, size, [=] MFEM_HOST_DEVICE (int i)
-   {
-      d_flag[i] = true;
-   });
-   const auto d_indices = indices.Read(use_dev);
-   mfem::forall_switch(use_dev, indices.Size(), [=] MFEM_HOST_DEVICE (int i)
-   {
-      d_flag[d_indices[i]] = false;
-   });
-
-   Array<int> out_idx(size);
-   auto d_out_idx = out_idx.Write(use_dev);
-   Array<char> workspace;
-   // Perform inclusive scan so that the last entry is the new size.
-   InclusiveScan(use_dev, d_flag, d_out_idx, size, workspace);
-
-   Vector copy(*this);
-   auto d_in = copy.Read(use_dev);
-   auto d_out = Write(use_dev);
-   mfem::forall_switch(use_dev, size, [=] MFEM_HOST_DEVICE (int i)
-   {
-      if (d_flag[i])
-      {
-         // Transform inclusive scan to exclusive by shifting.
-         const int j = (i > 0) ? d_out_idx[i - 1] : 0;
-         d_out[j] = d_in[i];
-      }
-   });
-
-   // Get the new size of the vector. Copy only the last entry.
-   Memory<int> submem(out_idx.GetMemory(), out_idx.Size() - 1, 1);
-   size = submem.Read(MemoryClass::HOST, 1)[0];
-}
-
 } // namespace mfem

linalg/vector.hpp

@@ -171,12 +171,6 @@ public:
    /// Resize the vector to size @a s using the MemoryType of @a v.
    void SetSize(int s, const Vector &v) { SetSize(s, v.GetMemory().GetMemoryType()); }
 
-   /// Update \ref Capacity() to @a res (if less than current), keeping existing entries.
-   void Reserve(int res);
-
-   /// Delete entries at @a indices and resize vector accordingly.
-   void DeleteAt(const Array<int> &indices);
-
    /// Set the Vector data.
    /// @warning This method should be called only when OwnsData() is false.
    void SetData(real_t *d) { data.Wrap(d, data.Capacity(), false); }
@@ -627,18 +621,6 @@ inline void Vector::SetSize(int s, MemoryType mt)
    data.UseDevice(use_dev);
 }
 
-inline void Vector::Reserve(int res)
-{
-   if (res > Capacity())
-   {
-      Memory<real_t> p(res, data.GetMemoryType());
-      p.CopyFrom(data, size);
-      p.UseDevice(data.UseDevice());
-      data.Delete();
-      data = p;
-   }
-}
-
 inline void Vector::NewMemoryAndSize(const Memory<real_t> &mem, int s,
                                      bool own_mem)
 {

mesh/mesh.cpp

@@ -979,46 +979,6 @@ const Array<int>& Mesh::GetElementAttributes() const
    return elem_attrs_cache;
 }
 
-void Mesh::ComputeFaceInfo(FaceType ftype) const
-{
-   auto &fidcs = face_indices[static_cast<int>(ftype)];
-   auto &ifidcs = inv_face_indices[static_cast<int>(ftype)];
-   fidcs.SetSize(GetNFbyType(ftype));
-   fidcs.HostWrite();
-   ifidcs.reserve(fidcs.Size());
-   int f_idx = 0;
-   for (int i = 0; i < GetNumFacesWithGhost(); ++i)
-   {
-      const FaceInformation face = GetFaceInformation(i);
-      if (face.IsNonconformingCoarse() || !face.IsOfFaceType(ftype))
-      {
-         continue;
-      }
-      fidcs[f_idx] = i;
-      ifidcs[i] = f_idx;
-      ++f_idx;
-   }
-}
-
-const Array<int> &Mesh::GetFaceIndices(FaceType ftype) const
-{
-   if (face_indices[static_cast<int>(ftype)].Size() == 0)
-   {
-      ComputeFaceInfo(ftype);
-   }
-   return face_indices[static_cast<int>(ftype)];
-}
-
-const std::unordered_map<int, int> &
-Mesh::GetInvFaceIndices(FaceType ftype) const
-{
-   if (inv_face_indices[static_cast<int>(ftype)].empty())
-   {
-      ComputeFaceInfo(ftype);
-   }
-   return inv_face_indices[static_cast<int>(ftype)];
-}
-
 void Mesh::DeleteGeometricFactors()
 {
    for (int i = 0; i < geom_factors.Size(); i++)
@@ -1906,11 +1866,6 @@ void Mesh::Destroy()
    bdr_face_attrs_cache.DeleteAll();
    attributes.DeleteAll();
    bdr_attributes.DeleteAll();
-
-   face_indices[0].DeleteAll();
-   face_indices[1].DeleteAll();
-   inv_face_indices[0] = std::unordered_map<int, int>();
-   inv_face_indices[1] = std::unordered_map<int, int>();
 }
 
 void Mesh::ResetLazyData()
@@ -8182,12 +8137,6 @@ void Mesh::GenerateFaces()
       FreeElement(f);
    }
 
-   // delete caches
-   face_indices[0].SetSize(0);
-   face_indices[1].SetSize(0);
-   inv_face_indices[0].clear();
-   inv_face_indices[1].clear();
-
    // (re)generate the interior faces and the info for them
    faces.SetSize(nfaces);
    faces_info.SetSize(nfaces);
@@ -10987,11 +10936,6 @@ void Mesh::Swap(Mesh& other, bool non_geometry)
    // copy attribute caches
    mfem::Swap(elem_attrs_cache, other.elem_attrs_cache);
    mfem::Swap(bdr_face_attrs_cache, other.bdr_face_attrs_cache);
-
-   mfem::Swap(face_indices[0], other.face_indices[0]);
-   mfem::Swap(face_indices[1], other.face_indices[1]);
-   inv_face_indices[0].swap(other.inv_face_indices[0]);
-   inv_face_indices[1].swap(other.inv_face_indices[1]);
 }
 
 void Mesh::GetElementData(const Array<Element*> &elem_array, int geom,

mesh/mesh.hpp

@@ -278,13 +278,6 @@ protected:
 
    // used during NC mesh initialization only
    Array<Triple<int, int, int> > tmp_vertex_parents;
-   /// cache for FaceIndices(ftype)
-   mutable Array<int> face_indices[2];
-   /// cache for FaceIndices(ftype)
-   mutable std::unordered_map<int, int> inv_face_indices[2];
-
-   /// compute face_indices[ftype] and inv_face_indices[type]
-   void ComputeFaceInfo(FaceType ftype) const;
 
 public:
    typedef Geometry::Constants<Geometry::SEGMENT>     seg_t;
@@ -319,11 +312,6 @@ public:
    // (true) is set in mesh_readers.cpp.
    static bool remove_unused_vertices;
 
-   /// Map from boundary or interior face indices to mesh face indices.
-   const Array<int>& GetFaceIndices(FaceType ftype) const;
-   /// Inverse of the map FaceIndices(ftype)
-   const std::unordered_map<int, int>& GetInvFaceIndices(FaceType ftype) const;
-
 protected:
    Operation last_operation;
 

tests/unit/CMakeLists.txt

@@ -22,7 +22,6 @@ set(UNIT_TESTS_SRCS
   dfem/test_mass.cpp
   general/test_array.cpp
   general/test_reduction.cpp
-  general/test_scan.cpp
   general/test_arrays_by_name.cpp
   general/test_error.cpp
   general/test_mem.cpp

tests/unit/general/test_array.cpp

@@ -124,19 +124,3 @@ TEST_CASE("Array stl-interactions", "[Array]")
       CHECK(x[i] == y[i]);
    }
 }
-
-TEST_CASE("Array delete at indices", "[Array]")
-{
-   Array<int>       test({0,1,2,3,4,5,6,7,8});
-   Array<int> rm_indices({0,    3,4,  6,  8});
-   Array<int>     result({  1,2,    5,  7  });
-
-   test.DeleteAt(rm_indices);
-
-   REQUIRE(test.Size() == result.Size());
-
-   for (int i = 0; i < test.Size(); i++)
-   {
-      CHECK(test[i] == result[i]);
-   }
-}

tests/unit/general/test_scan.cpp

@@ -1,102 +0,0 @@
-// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
-// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
-// LICENSE and NOTICE for details. LLNL-CODE-806117.
-//
-// This file is part of the MFEM library. For more information and source code
-// availability visit https://mfem.org.
-//
-// MFEM is free software; you can redistribute it and/or modify it under the
-// terms of the BSD-3 license. We welcome feedback and contributions, see file
-// CONTRIBUTING.md for details.
-
-#include <algorithm>
-#include <limits>
-
-#include "mfem.hpp"
-#include "unit_tests.hpp"
-
-// must be included after mfem.hpp
-#include "general/scan.hpp"
-
-using namespace mfem;
-
-TEST_CASE("Inclusive Scan", "[Scan],[GPU]")
-{
-   Array<char> workspace;
-   Array<int> a(10);
-
-   for (int use_dev = 0; use_dev < 2; ++use_dev)
-   {
-      CAPTURE(use_dev);
-      a.HostReadWrite();
-      for (int i = 0; i < a.Size(); ++i)
-      {
-         a[i] = i;
-      }
-      auto dptr = a.ReadWrite(use_dev);
-      InclusiveScan(use_dev, dptr, dptr, a.Size(), workspace);
-      a.HostRead();
-      for (int i = 0; i < a.Size(); ++i)
-      {
-         int expected = (i + 1) * i / 2;
-         CAPTURE(i);
-         REQUIRE(a[i] == expected);
-      }
-      a.HostReadWrite();
-      for (int i = 0; i < a.Size(); ++i)
-      {
-         a[i] = i + 1;
-      }
-      a.ReadWrite(use_dev);
-      InclusiveScan(use_dev, dptr, dptr, a.Size(), workspace, std::multiplies<> {});
-      a.HostRead();
-      int expected = 1;
-      for (int i = 0; i < a.Size(); ++i)
-      {
-         expected *= i + 1;
-         CAPTURE(i);
-         REQUIRE(a[i] == expected);
-      }
-   }
-}
-
-TEST_CASE("Exclusive Scan", "[Scan],[GPU]")
-{
-   Array<char> workspace;
-   Array<int> a(10);
-
-   for (int use_dev = 0; use_dev < 2; ++use_dev)
-   {
-      CAPTURE(use_dev);
-      a.HostReadWrite();
-      for (int i = 0; i < a.Size(); ++i)
-      {
-         a[i] = i;
-      }
-      auto dptr = a.ReadWrite(use_dev);
-      ExclusiveScan(use_dev, dptr, dptr, a.Size(), 5, workspace);
-      a.HostRead();
-      for (int i = 0; i < a.Size(); ++i)
-      {
-         int expected = (i + 1) * i / 2 - i + 5;
-         CAPTURE(i);
-         REQUIRE(a[i] == expected);
-      }
-      a.HostReadWrite();
-      for (int i = 0; i < a.Size(); ++i)
-      {
-         a[i] = i + 1;
-      }
-      a.ReadWrite(use_dev);
-      ExclusiveScan(use_dev, dptr, dptr, a.Size(), 5, workspace,
-                    std::multiplies<> {});
-      a.HostRead();
-      int expected = 5;
-      for (int i = 0; i < a.Size(); ++i)
-      {
-         CAPTURE(i);
-         REQUIRE(a[i] == expected);
-         expected *= i + 1;
-      }
-   }
-}

tests/unit/linalg/test_vector.cpp

@@ -247,21 +247,3 @@ TEST_CASE("Vector Sum", "[Vector],[GPU]")
 
    REQUIRE(sum_1 == MFEM_Approx(sum_2));
 }
-
-TEST_CASE("Vector delete at indices", "[Vector][GPU]")
-{
-   Vector           test({0,1,2,3,4,5,6,7,8});
-   Array<int> rm_indices({0,    3,4,  6,  8});
-   Vector         result({  1,2,    5,  7  });
-
-   test.UseDevice(true);
-   test.DeleteAt(rm_indices);
-
-   REQUIRE(test.Size() == result.Size());
-
-   test.HostReadWrite();
-   for (int i = 0; i < test.Size(); i++)
-   {
-      CHECK(test[i] == result[i]);
-   }
-}