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/dfem/doperator.hpp

@@ -231,53 +231,22 @@ public:
       const std::vector<FieldDescriptor> &parameters,
       const ParMesh &mesh);
 
-   /// MultLevel enum to indicate if the T->L Operators are used in the
-   /// Mult method.
-   enum MultLevel
-   {
-      TVECTOR,
-      LVECTOR
-   };
-
-   /// @brief Set the MultLevel mode for the DifferentiableOperator.
-   /// The default is TVECTOR, which means that the Operator will use
-   /// T->L before Mult and L->T Operators after.
-   void SetMultLevel(MultLevel level)
-   {
-      mult_level = level;
-   }
-
    /// @brief Compute the action of the operator on a given vector.
    ///
-   /// @param solutions_in The solution vector in which to compute the action.
-   /// This has to be a T-dof vector if MultLevel is set to TVECTOR, or L-dof
-   /// Vector if MultLevel is set to LVECTOR.
-   /// @param result_in Result vector of the action of the operator on
-   /// solutions. The result is a T-dof vector or L-dof vector depending on
-   /// the MultLevel.
-   void Mult(const Vector &solutions_in, Vector &result_in) const override
+   /// @param solutions_t The solution vector in which to compute the action.
+   /// This has to be a T-dof vector.
+   /// @param result_t Result vector of the action of the operator on
+   /// solutions_t. The result is a T-dof vector.
+   void Mult(const Vector &solutions_t, Vector &result_t) const override
    {
       MFEM_ASSERT(!action_callbacks.empty(), "no integrators have been set");
-
-      if (mult_level == MultLevel::LVECTOR)
+      prolongation(solutions, solutions_t, solutions_l);
+      residual_l = 0.0;
+      for (auto &action : action_callbacks)
       {
-         get_lvectors(solutions, solutions_in, solutions_l);
-         result_in = 0.0;
-         for (auto &action : action_callbacks)
-         {
-            action(solutions_l, parameters_l, result_in);
-         }
-      }
-      else
-      {
-         prolongation(solutions, solutions_in, solutions_l);
-         residual_l = 0.0;
-         for (auto &action : action_callbacks)
-         {
-            action(solutions_l, parameters_l, residual_l);
-         }
-         prolongation_transpose(residual_l, result_in);
+         action(solutions_l, parameters_l, residual_l);
       }
+      prolongation_transpose(residual_l, result_t);
    }
 
    /// @brief Add a domain integrator to the operator.
@@ -376,8 +345,6 @@ public:
 private:
    const ParMesh &mesh;
 
-   MultLevel mult_level = TVECTOR;
-
    std::vector<action_t> action_callbacks;
    std::map<size_t,
        std::vector<derivative_action_t>> derivative_action_callbacks;
@@ -387,6 +354,7 @@ private:
        std::vector<assemble_derivative_hypreparmatrix_callback_t>>
        assemble_derivative_hypreparmatrix_callbacks;
 
+
    std::vector<FieldDescriptor> solutions;
    std::vector<FieldDescriptor> parameters;
    // solutions and parameters

fem/dfem/util.hpp

@@ -1076,24 +1076,6 @@ void prolongation(const std::vector<FieldDescriptor> fields,
    }
 }
 
-inline
-void get_lvectors(const std::vector<FieldDescriptor> fields,
-                  const Vector &x,
-                  std::vector<Vector> &fields_l)
-{
-   int data_offset = 0;
-   for (std::size_t i = 0; i < fields.size(); i++)
-   {
-      const int sz = GetVSize(fields[i]);
-      fields_l[i].SetSize(sz);
-
-      const Vector x_i(const_cast<Vector&>(x), data_offset, sz);
-      fields_l[i] = x_i;
-
-      data_offset += sz;
-   }
-}
-
 /// @brief Get a transpose prolongation callback for a field descriptor.
 ///
 /// In the special case of a one field operator, the transpose prolongation

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); }
 

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

tests/unit/CMakeLists.txt

@@ -20,7 +20,6 @@ set(UNIT_TESTS_SRCS
   dfem/test_diffusion.cpp
   dfem/test_divergence.cpp
   dfem/test_mass.cpp
-  dfem/test_lvector_interface.cpp
   general/test_array.cpp
   general/test_reduction.cpp
   general/test_arrays_by_name.cpp
@@ -145,11 +144,11 @@ set(UNIT_TESTS_SRCS
 # SERIAL CPU TESTS: unit_tests
 #-----------------------------------------------------------
 if (MFEM_USE_CUDA)
-    set_property(SOURCE unit_test_main.cpp ${UNIT_TESTS_SRCS}
+   set_property(SOURCE unit_test_main.cpp ${UNIT_TESTS_SRCS}
                 PROPERTY LANGUAGE CUDA)
 endif()
 if (MFEM_USE_HIP)
-    set_property(SOURCE unit_test_main.cpp ${UNIT_TESTS_SRCS}
+  set_property(SOURCE unit_test_main.cpp ${UNIT_TESTS_SRCS}
                PROPERTY HIP_SOURCE_PROPERTY_FORMAT TRUE)
 endif()
 
@@ -176,7 +175,7 @@ COMMAND ${CMAKE_COMMAND} -E copy_directory
 #   make unit_tests
 #   ctest -R unit_tests [-V]
 if (MFEM_USE_DOUBLE) # otherwise returns MFEM_SKIP_RETURN_VALUE
-    add_test(NAME unit_tests COMMAND unit_tests)
+  add_test(NAME unit_tests COMMAND unit_tests)
 endif()
 
 #-----------------------------------------------------------
@@ -184,16 +183,16 @@ endif()
 #-----------------------------------------------------------
 # Create CUDA executable and test
 if (MFEM_USE_CUDA)
-    # gpu_unit_tests
-    set(GPU_UNIT_TESTS_SRCS gpu_unit_test_main.cpp)
-    set_property(SOURCE ${GPU_UNIT_TESTS_SRCS} PROPERTY LANGUAGE CUDA)
-    mfem_add_executable(gpu_unit_tests ${GPU_UNIT_TESTS_SRCS} ${UNIT_TESTS_SRCS})
-    target_link_libraries(gpu_unit_tests mfem)
-    add_dependencies(gpu_unit_tests copy_data)
-    add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} gpu_unit_tests)
-    if (MFEM_USE_DOUBLE) # otherwise returns MFEM_SKIP_RETURN_VALUE
-        add_test(NAME gpu_unit_tests COMMAND gpu_unit_tests)
-    endif()
+   # gpu_unit_tests
+   set(GPU_UNIT_TESTS_SRCS gpu_unit_test_main.cpp)
+   set_property(SOURCE ${GPU_UNIT_TESTS_SRCS} PROPERTY LANGUAGE CUDA)
+   mfem_add_executable(gpu_unit_tests ${GPU_UNIT_TESTS_SRCS} ${UNIT_TESTS_SRCS})
+   target_link_libraries(gpu_unit_tests mfem)
+   add_dependencies(gpu_unit_tests copy_data)
+   add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} gpu_unit_tests)
+   if (MFEM_USE_DOUBLE) # otherwise returns MFEM_SKIP_RETURN_VALUE
+      add_test(NAME gpu_unit_tests COMMAND gpu_unit_tests)
+   endif()
 endif()
 
 #-----------------------------------------------------------
@@ -201,15 +200,15 @@ endif()
 #-----------------------------------------------------------
 # Create HIP 'gpu_unit_tests' executable and test
 if (MFEM_USE_HIP)
-    # gpu_unit_tests
-    set(GPU_UNIT_TESTS_SRCS gpu_unit_test_main.cpp)
-    mfem_add_executable(gpu_unit_tests ${GPU_UNIT_TESTS_SRCS} ${UNIT_TESTS_SRCS})
-    target_link_libraries(gpu_unit_tests mfem)
-    add_dependencies(gpu_unit_tests copy_data)
-    add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} gpu_unit_tests)
-    if (MFEM_USE_DOUBLE) # otherwise returns MFEM_SKIP_RETURN_VALUE
-        add_test(NAME gpu_unit_tests COMMAND gpu_unit_tests)
-    endif()
+   # gpu_unit_tests
+   set(GPU_UNIT_TESTS_SRCS gpu_unit_test_main.cpp)
+   mfem_add_executable(gpu_unit_tests ${GPU_UNIT_TESTS_SRCS} ${UNIT_TESTS_SRCS})
+   target_link_libraries(gpu_unit_tests mfem)
+   add_dependencies(gpu_unit_tests copy_data)
+   add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} gpu_unit_tests)
+   if (MFEM_USE_DOUBLE) # otherwise returns MFEM_SKIP_RETURN_VALUE
+      add_test(NAME gpu_unit_tests COMMAND gpu_unit_tests)
+   endif()
 endif()
 
 #-----------------------------------------------------------
@@ -225,7 +224,7 @@ function(add_serial_miniapp_test name test_uvm)
 
     set(${NAME}_TESTS_SRCS miniapps/test_${name}.cpp)
     if (MFEM_USE_CUDA)
-        set_property(SOURCE ${${NAME}_TESTS_SRCS} PROPERTY LANGUAGE CUDA)
+       set_property(SOURCE ${${NAME}_TESTS_SRCS} PROPERTY LANGUAGE CUDA)
     endif(MFEM_USE_CUDA)
 
     mfem_add_executable(${name}_tests_cpu ${${NAME}_TESTS_SRCS})
@@ -245,25 +244,25 @@ function(add_serial_miniapp_test name test_uvm)
     endif()
 
     if (MFEM_USE_CUDA OR MFEM_USE_HIP)
-        mfem_add_executable(${name}_tests_gpu ${${NAME}_TESTS_SRCS})
-        target_compile_definitions(${name}_tests_gpu PUBLIC MFEM_${NAME}_DEVICE="gpu")
-        target_link_libraries(${name}_tests_gpu mfem)
-        add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} ${name}_tests_gpu)
-        if (MFEM_USE_DOUBLE) # otherwise returns MFEM_SKIP_RETURN_VALUE
-            add_test(NAME ${name}_tests_gpu COMMAND ${name}_tests_gpu)
-        endif()
+       mfem_add_executable(${name}_tests_gpu ${${NAME}_TESTS_SRCS})
+       target_compile_definitions(${name}_tests_gpu PUBLIC MFEM_${NAME}_DEVICE="gpu")
+       target_link_libraries(${name}_tests_gpu mfem)
+       add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} ${name}_tests_gpu)
+       if (MFEM_USE_DOUBLE) # otherwise returns MFEM_SKIP_RETURN_VALUE
+           add_test(NAME ${name}_tests_gpu COMMAND ${name}_tests_gpu)
+       endif()
 
-        if (test_uvm)
-            mfem_add_executable(${name}_tests_gpu_uvm ${${NAME}_TESTS_SRCS})
-            target_compile_definitions(${name}_tests_gpu_uvm PUBLIC
+       if (test_uvm)
+           mfem_add_executable(${name}_tests_gpu_uvm ${${NAME}_TESTS_SRCS})
+           target_compile_definitions(${name}_tests_gpu_uvm PUBLIC
                MFEM_${NAME}_DEVICE="gpu:uvm")
-            target_link_libraries(${name}_tests_gpu_uvm mfem)
-            add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME}
+           target_link_libraries(${name}_tests_gpu_uvm mfem)
+           add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME}
                ${name}_tests_gpu_uvm)
-            if (MFEM_USE_DOUBLE) # otherwise returns MFEM_SKIP_RETURN_VALUE
-                add_test(NAME ${name}_tests_gpu_uvm COMMAND ${name}_tests_gpu_uvm)
-            endif()
-        endif()
+           if (MFEM_USE_DOUBLE) # otherwise returns MFEM_SKIP_RETURN_VALUE
+               add_test(NAME ${name}_tests_gpu_uvm COMMAND ${name}_tests_gpu_uvm)
+           endif()
+       endif()
     endif()
 endfunction(add_serial_miniapp_test)
 
@@ -278,25 +277,25 @@ add_dependencies(tmop_pa_tests_cpu copy_miniapps_meshing_data)
 #-----------------------------------------------------------
 # Add 'ceed_tests' executable and test; add extra tests 'ceed_test_*'
 if (MFEM_USE_CEED)
-    set(CEED_TESTS_SRCS
+   set(CEED_TESTS_SRCS
       ceed/test_ceed.cpp
       ceed/test_ceed_main.cpp)
-    if (MFEM_USE_CUDA)
-        set_property(SOURCE ${CEED_TESTS_SRCS} PROPERTY LANGUAGE CUDA)
-    endif(MFEM_USE_CUDA)
-    mfem_add_executable(ceed_tests ${CEED_TESTS_SRCS})
-    target_link_libraries(ceed_tests mfem)
-    add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} ceed_tests)
-    # Add CEED tests
-    add_test(NAME ceed_tests COMMAND ceed_tests)
-    if (MFEM_USE_CUDA)
-        add_test(NAME ceed_tests_cuda_ref
+   if (MFEM_USE_CUDA)
+      set_property(SOURCE ${CEED_TESTS_SRCS} PROPERTY LANGUAGE CUDA)
+   endif(MFEM_USE_CUDA)
+   mfem_add_executable(ceed_tests ${CEED_TESTS_SRCS})
+   target_link_libraries(ceed_tests mfem)
+   add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} ceed_tests)
+   # Add CEED tests
+   add_test(NAME ceed_tests COMMAND ceed_tests)
+   if (MFEM_USE_CUDA)
+      add_test(NAME ceed_tests_cuda_ref
                COMMAND ceed_tests --device ceed-cuda:/gpu/cuda/ref)
-        add_test(NAME ceed_tests_cuda_shared
+      add_test(NAME ceed_tests_cuda_shared
                COMMAND ceed_tests --device ceed-cuda:/gpu/cuda/shared)
-        add_test(NAME ceed_tests_cuda_gen
+      add_test(NAME ceed_tests_cuda_gen
                COMMAND ceed_tests --device ceed-cuda:/gpu/cuda/gen)
-    endif()
+   endif()
 endif()
 
 #-----------------------------------------------------------
@@ -304,54 +303,54 @@ endif()
 #-----------------------------------------------------------
 # Define executables and tests
 if (MFEM_USE_MPI)
-    # punit_tests
-    if (MFEM_USE_CUDA)
-        set_property(SOURCE punit_test_main.cpp PROPERTY LANGUAGE CUDA)
-    endif()
-    mfem_add_executable(punit_tests punit_test_main.cpp ${UNIT_TESTS_SRCS})
-    target_link_libraries(punit_tests mfem)
-    add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} punit_tests)
-    foreach(np 1 ${MFEM_MPI_NP})
-        if (MFEM_USE_DOUBLE) # otherwise returns MFEM_SKIP_RETURN_VALUE
-            add_test(NAME punit_tests_np=${np}
+   # punit_tests
+   if (MFEM_USE_CUDA)
+      set_property(SOURCE punit_test_main.cpp PROPERTY LANGUAGE CUDA)
+   endif()
+   mfem_add_executable(punit_tests punit_test_main.cpp ${UNIT_TESTS_SRCS})
+   target_link_libraries(punit_tests mfem)
+   add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} punit_tests)
+   foreach(np 1 ${MFEM_MPI_NP})
+      if (MFEM_USE_DOUBLE) # otherwise returns MFEM_SKIP_RETURN_VALUE
+         add_test(NAME punit_tests_np=${np}
                   COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${np}
                   ${MPIEXEC_PREFLAGS} $<TARGET_FILE:punit_tests>
                   ${MPIEXEC_POSTFLAGS})
-        endif()
-    endforeach()
-    if (MFEM_USE_CUDA)
-        # pgpu_unit_tests
-        set(PGPU_UNIT_TESTS_SRCS pgpu_unit_test_main.cpp)
-        set_property(SOURCE ${PGPU_UNIT_TESTS_SRCS} PROPERTY LANGUAGE CUDA)
-        mfem_add_executable(pgpu_unit_tests ${PGPU_UNIT_TESTS_SRCS} ${UNIT_TESTS_SRCS})
-        add_dependencies(pgpu_unit_tests copy_data)
-        target_link_libraries(pgpu_unit_tests mfem)
-        add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} pgpu_unit_tests)
-        foreach(np 1 ${MFEM_MPI_NP})
-            if (MFEM_USE_DOUBLE) # otherwise returns MFEM_SKIP_RETURN_VALUE
-                add_test(NAME pgpu_unit_tests_np=${np}
+      endif()
+   endforeach()
+   if (MFEM_USE_CUDA)
+       # pgpu_unit_tests
+       set(PGPU_UNIT_TESTS_SRCS pgpu_unit_test_main.cpp)
+       set_property(SOURCE ${PGPU_UNIT_TESTS_SRCS} PROPERTY LANGUAGE CUDA)
+       mfem_add_executable(pgpu_unit_tests ${PGPU_UNIT_TESTS_SRCS} ${UNIT_TESTS_SRCS})
+       add_dependencies(pgpu_unit_tests copy_data)
+       target_link_libraries(pgpu_unit_tests mfem)
+       add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} pgpu_unit_tests)
+       foreach(np 1 ${MFEM_MPI_NP})
+          if (MFEM_USE_DOUBLE) # otherwise returns MFEM_SKIP_RETURN_VALUE
+             add_test(NAME pgpu_unit_tests_np=${np}
                       COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${np}
                       ${MPIEXEC_PREFLAGS} $<TARGET_FILE:pgpu_unit_tests>
                       ${MPIEXEC_POSTFLAGS})
-            endif()
-        endforeach()
+          endif()
+       endforeach()
     endif()
     if (MFEM_USE_HIP)
-        # pgpu_unit_tests
-        set(PGPU_UNIT_TESTS_SRCS pgpu_unit_test_main.cpp)
-        mfem_add_executable(pgpu_unit_tests ${PGPU_UNIT_TESTS_SRCS} ${UNIT_TESTS_SRCS})
-        add_dependencies(pgpu_unit_tests copy_data)
-        target_link_libraries(pgpu_unit_tests mfem)
-        add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} pgpu_unit_tests)
-        foreach(np 1 ${MFEM_MPI_NP})
-            if (MFEM_USE_DOUBLE) # otherwise returns MFEM_SKIP_RETURN_VALUE
-                add_test(NAME pgpu_unit_tests_np=${np}
+       # pgpu_unit_tests
+       set(PGPU_UNIT_TESTS_SRCS pgpu_unit_test_main.cpp)
+       mfem_add_executable(pgpu_unit_tests ${PGPU_UNIT_TESTS_SRCS} ${UNIT_TESTS_SRCS})
+       add_dependencies(pgpu_unit_tests copy_data)
+       target_link_libraries(pgpu_unit_tests mfem)
+       add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} pgpu_unit_tests)
+       foreach(np 1 ${MFEM_MPI_NP})
+          if (MFEM_USE_DOUBLE) # otherwise returns MFEM_SKIP_RETURN_VALUE
+             add_test(NAME pgpu_unit_tests_np=${np}
                       COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${np}
                       ${MPIEXEC_PREFLAGS} $<TARGET_FILE:pgpu_unit_tests>
                       ${MPIEXEC_POSTFLAGS})
-            endif()
-        endforeach()
-    endif()
+          endif()
+       endforeach()
+   endif()
 endif(MFEM_USE_MPI)
 
 #-----------------------------------------------------------
@@ -425,8 +424,8 @@ endfunction(add_parallel_miniapp_test)
 
 # Additional MPI unit tests
 if (MFEM_USE_MPI)
-    add_parallel_miniapp_test(sedov TRUE)
-    add_parallel_miniapp_test(tmop_pa FALSE)
+   add_parallel_miniapp_test(sedov TRUE)
+   add_parallel_miniapp_test(tmop_pa FALSE)
 endif(MFEM_USE_MPI)
 
 #-----------------------------------------------------------
@@ -435,10 +434,10 @@ endif(MFEM_USE_MPI)
 #-----------------------------------------------------------
 set(DEBUG_DEVICE_SRCS miniapps/test_debug_device.cpp)
 if (MFEM_USE_CUDA)
-    set_property(SOURCE ${DEBUG_DEVICE_SRCS} PROPERTY LANGUAGE CUDA)
+   set_property(SOURCE ${DEBUG_DEVICE_SRCS} PROPERTY LANGUAGE CUDA)
 endif()
 if (MFEM_USE_HIP)
-    set_property(SOURCE ${DEBUG_DEVICE_SRCS}
+   set_property(SOURCE ${DEBUG_DEVICE_SRCS}
                 PROPERTY HIP_SOURCE_PROPERTY_FORMAT TRUE)
 endif()
 mfem_add_executable(debug_device_tests ${DEBUG_DEVICE_SRCS})

tests/unit/dfem/test_diffusion.cpp

@@ -255,8 +255,8 @@ void DFemDiffusion(const char *filename, int p, const int r)
          DOperator dop_mf(vsol, {{Coords, mfes}}, pmesh);
          const auto mf_vector_diffusion_qf =
             [] MFEM_HOST_DEVICE (const tensor<dscalar_t, DIM, DIM> &dudxi,
-                                 const tensor<real_t, DIM, DIM> &J,
-                                 const real_t &w)
+               const tensor<real_t, DIM, DIM> &J,
+               const real_t &w)
          {
             const auto invJ = inv(J), TinJ = transpose(invJ);
             return tuple{ (dudxi * invJ) * TinJ * det(J) * w };

tests/unit/dfem/test_lvector_interface.cpp

@@ -1,98 +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 "../unit_tests.hpp"
-#include "mfem.hpp"
-#include <utility>
-
-#ifdef MFEM_USE_MPI
-
-using namespace mfem;
-using namespace mfem::future;
-using mfem::future::tensor;
-
-constexpr int DIM = 3;
-
-namespace kernels
-{
-struct MFApply
-{
-   MFEM_HOST_DEVICE inline auto operator()(const tensor<real_t, DIM> &dudxi,
-                                           const tensor<real_t, DIM, DIM> &J,
-                                           const real_t &w) const
-   {
-      const auto invJ = inv(J);
-      return tuple{ (dudxi * invJ) * transpose(invJ) * det(J) * w };
-   }
-};
-}
-
-TEST_CASE("DFEM L-Vector interface", "[Parallel][DFEM]")
-{
-   constexpr int p = 2; // Polynomial order
-   constexpr int r = 1;
-   constexpr int q = 2 * p + r;
-
-   const auto filename = GENERATE("../../data/fichera.mesh");
-   Mesh smesh(filename);
-   ParMesh pmesh(MPI_COMM_WORLD, smesh);
-   MFEM_VERIFY(pmesh.Dimension() == DIM, "Mesh dimension mismatch");
-
-   pmesh.EnsureNodes();
-   auto *nodes = static_cast<ParGridFunction *>(pmesh.GetNodes());
-   smesh.Clear();
-
-   Array<int> all_domain_attr;
-   if (pmesh.attributes.Size() > 0)
-   {
-      all_domain_attr.SetSize(pmesh.attributes.Max());
-      all_domain_attr = 1;
-   }
-
-   H1_FECollection fec(p, DIM);
-   ParFiniteElementSpace pfes(&pmesh, &fec);
-   ParFiniteElementSpace *mfes = nodes->ParFESpace();
-
-   const auto *ir = &IntRules.Get(pmesh.GetTypicalElementGeometry(), q);
-
-   ParGridFunction x(&pfes), y(&pfes), z(&pfes);
-   Vector X(pfes.GetTrueVSize()), Y(pfes.GetTrueVSize()), Z(pfes.GetTrueVSize());
-
-   X.Randomize(1);
-   x.SetFromTrueDofs(X);
-
-   ParBilinearForm blf_fa(&pfes);
-   blf_fa.AddDomainIntegrator(new DiffusionIntegrator(ir));
-   blf_fa.Assemble();
-   blf_fa.Finalize();
-
-   static constexpr int U = 0, Coords = 1;
-
-   const auto solution = std::vector{FieldDescriptor{U, &pfes}};
-   DifferentiableOperator dop(solution, {{Coords, mfes}}, pmesh);
-
-   kernels::MFApply mf_apply_qf;
-   dop.AddDomainIntegrator(mf_apply_qf,
-                           tuple{Gradient<U>{}, Gradient<Coords>{}, Weight{}},
-                           tuple{Gradient<U>{}}, *ir, all_domain_attr);
-
-   // Use the L-vector interface to multiply
-   dop.SetMultLevel(DifferentiableOperator::MultLevel::LVECTOR);
-   dop.SetParameters({nodes});
-   dop.Mult(x, z);
-
-   blf_fa.Mult(x, y);
-
-   z -= y;
-   REQUIRE(z.Normlinf() == MFEM_Approx(0.0));
-}
-
-#endif