Add unit test for operator ownership in sum and product operators

SumOperator, ProductOperator, and TripleProductOperator use Handle<T> instead of
raw pointers. Explicit ownership flags and destructors can be removed. The
classes now have proper copy and move semantics (rule of zero).

Retain the constructors with explicit ownership flags for backwards
compatibility.

fem/multigrid.hpp

@@ -16,7 +16,7 @@
 #include "bilinearform.hpp"
 
 #include "../linalg/operator.hpp"
-#include "../linalg/handle.hpp"
+#include "../linalg/op_handle.hpp"
 
 namespace mfem
 {

general/CMakeLists.txt

@@ -45,6 +45,7 @@ list(APPEND HDRS
   gecko.hpp
   globals.hpp
   zstr.hpp
+  handle.hpp
   hash.hpp
   isockstream.hpp
   kdtree.hpp

general/handle.hpp

@@ -0,0 +1,200 @@
+// Copyright (c) 2010-2024, 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_HANDLE_HPP
+#define MFEM_HANDLE_HPP
+
+#include "../config/config.hpp"
+#include <memory>
+
+namespace mfem
+{
+
+/// @brief A smart pointer class that may represent either shared ownership, or
+/// a non-owning borrow.
+///
+/// A Handle may either be owning or non-owning. Non-owning Handle%s point to
+/// externally owned data; it is the responsibility of the user both to ensure
+/// that the data remains valid as long as the Handle is alive and to delete the
+/// pointer when it is no longer needed. Owning Handle%s use <a
+/// href="https://en.cppreference.com/w/cpp/memory/shared_ptr">
+/// std::shared_ptr</a> to implement reference counting. The underlying data
+/// will be valid as long as there is at least one live copy. When the last
+/// Handle is destroyed, the pointer is deleted.
+///
+/// Both types of Handle%s can be copied, moved, stored in standard containers,
+/// etc.
+///
+/// A non-owning Handle may assume ownership over its data, but an owning Handle
+/// cannot release ownership over its data.
+///
+/// It is an invariant of this class that **at most** one of the data members
+/// @a not_owned and @a owned will be non-null.
+template <typename T>
+class Handle
+{
+   /// If this is a non-owning handle, @a not_owned will point to the data.
+   T *not_owned = nullptr;
+
+   /// If this is an owning handle, @a owned will point to the data.
+   std::shared_ptr<T> owned = nullptr;
+
+   /// @brief Types @a Handle<T> and @a %Handle\<U\> are friends to allow
+   /// construction of one from another when @a T and @a U are convertible
+   /// types.
+   template <typename U> friend class Handle;
+
+public:
+   /// Create an empty (null) Handle.
+   Handle() = default;
+
+   /// @brief Create a Handle pointing to @a t.
+   ///
+   /// If @a take_ownership is true, then the Handle assumes ownership over the
+   /// pointer, and it should not be deleted externally. Otherwise, the Handle
+   /// will be non-owning, and it is the user's responsibility to ensure the
+   /// correct lifetime of @a t.
+   Handle(T *t, bool take_ownership)
+   {
+      if (take_ownership) { owned.reset(t); }
+      else { not_owned = t; }
+   }
+
+   /// Create a Handle from a std::shared_ptr (sharing ownership with @a t).
+   Handle(const std::shared_ptr<T> &t) : owned(t) { }
+
+   /// @brief Copy constructor.
+   ///
+   /// Copying an owning Handle results in another owning handle. Copying a
+   /// non-owning handle results in a non-owning handle.
+   Handle(const Handle &other) = default;
+
+   /// Move constructor (see Handle(const Handle&)).
+   Handle(Handle &&other) = default;
+
+   /// @brief Constructs a copy of @a u, where type @a U is convertible to @a T.
+   ///
+   /// This allows the construction of Handle<Base> from Handle<Derived>.
+   template <typename U>
+   Handle(const Handle<U> &u) : not_owned(u.not_owned), owned(u.owned) { }
+
+   /// @brief Move-constructs from @a u, where type @a U is convertible to @a T.
+   ///
+   /// See @ref Handle(const Handle<U>&).
+   template <typename U>
+   Handle(Handle<U> &&u) : not_owned(u.not_owned), owned(u.owned) { }
+
+   /// Destructor. If the Handle is owning, decrement the reference count.
+   ~Handle() = default;
+
+   /// Copy assignment (see Handle(const Handle&)).
+   Handle &operator=(const Handle &other) = default;
+
+   /// Move assignment (see Handle(const Handle&)).
+   Handle &operator=(Handle &&other) = default;
+
+   /// Returns the contained pointer (may be null).
+   T *Get() const
+   {
+      if (not_owned) { return not_owned; }
+      else { return owned.get(); }
+   }
+
+   /// @brief If the Handle is owning, return a copy of the underlying shared
+   /// pointer.
+   ///
+   /// @warning If the Handle is non-owning (even if non-null), this will return
+   /// and empty (null) shared pointer.
+   std::shared_ptr<T> GetSharedPtr() const { return owned; }
+
+   /// Dereference operator. The Handle must be non-null.
+   T &operator*() const { return *Get(); }
+
+   /// Member access (arrow) operator. The Handle must be non-null.
+   T *operator->() const { return Get(); }
+
+   /// @brief Returns true if the Handle is owning, false if it is non-owning.
+   ///
+   /// Returns false if the Handle is null (empty).
+   bool IsOwner() const { return owned; }
+
+   /// Returns true if the Handle is non-null.
+   explicit operator bool() const { return not_owned || owned; }
+
+   /// @brief Assume owernship of the data.
+   ///
+   /// If the Handle is already owning, this does nothing.
+   void MakeOwner()
+   {
+      if (owned) { return; }
+      owned.reset(not_owned);
+      not_owned = nullptr;
+   }
+
+   /// @brief Reset the Handle to be empty.
+   ///
+   /// If the Handle is owning, this will decrement the reference count.
+   void Reset()
+   {
+      owned.reset();
+      not_owned = nullptr;
+   }
+
+   /// @brief Reset the Handle to point to @a t.
+   ///
+   /// The Handle may assume ownership of the pointer according to @a
+   /// take_ownership (see @ref Handle(T*, bool)).
+   void Reset(T *t, bool take_ownership)
+   {
+      if (take_ownership)
+      {
+         owned.reset(t);
+         not_owned = nullptr;
+      }
+      else
+      {
+         owned.reset();
+         not_owned = t;
+      }
+   }
+
+   /// Reset the Handle to share ownership with @a t.
+   void Reset(const std::shared_ptr<T> &t)
+   {
+      owned = t;
+      not_owned = nullptr;
+   }
+};
+
+/// @brief Return a new owning Handle, where the pointed-to object is a new
+/// object constructed using the given arguments.
+///
+/// This is analogous to <a
+/// href="https://en.cppreference.com/w/cpp/memory/shared_ptr/make_shared">
+/// std::make_shared</a>.
+template <typename T, typename... Args>
+Handle<T> MakeOwning(Args&&... args)
+{
+   T *t = new T(std::forward<Args>(args)...);
+   return Handle<T>(t, true);
+}
+
+/// Return a new owning Handle pointing to @a t.
+template <typename T>
+Handle<T> Owning(T *t) { return Handle<T>(t, true); }
+
+/// Return a new non-owning Handle pointing to @a t.
+template <typename T>
+Handle<T> NonOwning(T *t) { return Handle<T>(t, false); }
+
+} // namespace mfem
+
+#endif

linalg/CMakeLists.txt

@@ -24,7 +24,7 @@ list(APPEND SRCS
   constraints.cpp
   densemat.cpp
   symmat.cpp
-  handle.cpp
+  op_handle.cpp
   matrix.cpp
   ode.cpp
   operator.cpp
@@ -51,7 +51,7 @@ list(APPEND HDRS
   dinvariants.hpp
   symmat.hpp
   dtensor.hpp
-  handle.hpp
+  op_handle.hpp
   invariants.hpp
   kernels.hpp
   lapack.hpp

linalg/linalg.hpp

@@ -28,7 +28,7 @@
 #include "symmat.hpp"
 #include "ode.hpp"
 #include "solvers.hpp"
-#include "handle.hpp"
+#include "op_handle.hpp"
 #include "invariants.hpp"
 #include "constraints.hpp"
 #include "auxiliary.hpp"

linalg/op_handle.cpp

@@ -9,7 +9,7 @@
 // terms of the BSD-3 license. We welcome feedback and contributions, see file
 // CONTRIBUTING.md for details.
 
-#include "handle.hpp"
+#include "op_handle.hpp"
 #include "sparsemat.hpp"
 #ifdef MFEM_USE_MPI
 #include "petsc.hpp"

linalg/op_handle.hpp

@@ -9,8 +9,8 @@
 // terms of the BSD-3 license. We welcome feedback and contributions, see file
 // CONTRIBUTING.md for details.
 
-#ifndef MFEM_HANDLE_HPP
-#define MFEM_HANDLE_HPP
+#ifndef MFEM_OP_HANDLE_HPP
+#define MFEM_OP_HANDLE_HPP
 
 #include "../config/config.hpp"
 #include "operator.hpp"

linalg/operator.cpp

@@ -183,8 +183,8 @@ Operator * Operator::SetupRAP(const Operator *Pi, const Operator *Po)
    {
       if (!IsIdentityProlongation(Po))
       {
-         TransposeOperator * PoT = new TransposeOperator(Po);
-         rap = new ProductOperator(PoT, this, true,false);
+         rap = new ProductOperator(Owning(new TransposeOperator(Po)),
+                                   NonOwning(this));
       }
       else
       {
@@ -365,11 +365,10 @@ void SecondOrderTimeDependentOperator::ImplicitSolve(const real_t dt0,
    mfem_error("SecondOrderTimeDependentOperator::ImplicitSolve() is not overridden!");
 }
 
-SumOperator::SumOperator(const Operator *A, const real_t alpha,
-                         const Operator *B, const real_t beta,
-                         bool ownA, bool ownB)
-   : Operator(A->Height(), A->Width()),
-     A(A), B(B), alpha(alpha), beta(beta), ownA(ownA), ownB(ownB),
+SumOperator::SumOperator(Handle<const Operator> A_, const real_t alpha,
+                         Handle<const Operator> B_, const real_t beta)
+   : Operator(A_->Height(), A_->Width()),
+     A(A_), B(B_), alpha(alpha), beta(beta),
      z(A->Height())
 {
    MFEM_VERIFY(A->Width() == B->Width(),
@@ -381,53 +380,43 @@ SumOperator::SumOperator(const Operator *A, const real_t alpha,
                << "A->Height() = " << A->Height()
                << ", B->Height() = " << B->Height() );
 
+   if (auto SolverA = dynamic_cast<const Solver*>(A.Get()))
    {
-      const Solver* SolverA = dynamic_cast<const Solver*>(A);
-      const Solver* SolverB = dynamic_cast<const Solver*>(B);
-      if (SolverA)
-      {
-         MFEM_VERIFY(!(SolverA->iterative_mode),
-                     "Operator A of a SumOperator should not be in iterative mode");
-      }
-      if (SolverB)
-      {
-         MFEM_VERIFY(!(SolverB->iterative_mode),
-                     "Operator B of a SumOperator should not be in iterative mode");
-      }
+      MFEM_VERIFY(!(SolverA->iterative_mode),
+                  "Operator A of a SumOperator should not be in iterative mode");
+   }
+   if (auto SolverB = dynamic_cast<const Solver*>(B.Get()))
+   {
+      MFEM_VERIFY(!(SolverB->iterative_mode),
+                  "Operator B of a SumOperator should not be in iterative mode");
    }
-
 }
 
-SumOperator::~SumOperator()
-{
-   if (ownA) { delete A; }
-   if (ownB) { delete B; }
-}
+SumOperator::SumOperator(const Operator *A_, const real_t alpha,
+                         const Operator *B_, const real_t beta,
+                         bool own_A, bool own_B)
+   : SumOperator({A_, own_A}, alpha, {B_, own_B}, beta) { }
 
-ProductOperator::ProductOperator(const Operator *A, const Operator *B,
-                                 bool ownA, bool ownB)
-   : Operator(A->Height(), B->Width()),
-     A(A), B(B), ownA(ownA), ownB(ownB), z(A->Width())
+
+ProductOperator::ProductOperator(Handle<const Operator> A_,
+                                 Handle<const Operator> B_)
+   : Operator(A_->Height(), B_->Width()),
+     A(A_), B(B_), z(A->Width())
 {
    MFEM_VERIFY(A->Width() == B->Height(),
                "incompatible Operators: A->Width() = " << A->Width()
                << ", B->Height() = " << B->Height());
 
+   if (auto SolverB = dynamic_cast<const Solver*>(B.Get()))
    {
-      const Solver* SolverB = dynamic_cast<const Solver*>(B);
-      if (SolverB)
-      {
-         MFEM_VERIFY(!(SolverB->iterative_mode),
-                     "Operator B of a ProductOperator should not be in iterative mode");
-      }
+      MFEM_VERIFY(!(SolverB->iterative_mode),
+                  "Operator B of a ProductOperator should not be in iterative mode");
    }
 }
 
-ProductOperator::~ProductOperator()
-{
-   if (ownA) { delete A; }
-   if (ownB) { delete B; }
-}
+ProductOperator::ProductOperator(const Operator *A_, const Operator *B_,
+                                 bool own_A, bool own_B)
+   : ProductOperator({A_, own_A}, {B_, own_B}) { }
 
 
 RAPOperator::RAPOperator(const Operator &Rt_, const Operator &A_,
@@ -465,11 +454,9 @@ RAPOperator::RAPOperator(const Operator &Rt_, const Operator &A_,
 
 
 TripleProductOperator::TripleProductOperator(
-   const Operator *A, const Operator *B, const Operator *C,
-   bool ownA, bool ownB, bool ownC)
-   : Operator(A->Height(), C->Width())
-   , A(A), B(B), C(C)
-   , ownA(ownA), ownB(ownB), ownC(ownC)
+   Handle<const Operator> A_, Handle<const Operator> B_, Handle<const Operator> C_)
+   : Operator(A_->Height(), C_->Width()),
+     A(A_), B(B_), C(C_)
 {
    MFEM_VERIFY(A->Width() == B->Height(),
                "incompatible Operators: A->Width() = " << A->Width()
@@ -478,20 +465,16 @@ TripleProductOperator::TripleProductOperator(
                "incompatible Operators: B->Width() = " << B->Width()
                << ", C->Height() = " << C->Height());
 
+   if (auto SolverB = dynamic_cast<const Solver*>(B.Get()))
    {
-      const Solver* SolverB = dynamic_cast<const Solver*>(B);
-      if (SolverB)
-      {
-         MFEM_VERIFY(!(SolverB->iterative_mode),
-                     "Operator B of a TripleProductOperator should not be in iterative mode");
-      }
+      MFEM_VERIFY(!(SolverB->iterative_mode),
+                  "Operator B of a TripleProductOperator should not be in iterative mode");
+   }
 
-      const Solver* SolverC = dynamic_cast<const Solver*>(C);
-      if (SolverC)
-      {
-         MFEM_VERIFY(!(SolverC->iterative_mode),
-                     "Operator C of a TripleProductOperator should not be in iterative mode");
-      }
+   if (auto SolverC = dynamic_cast<const Solver*>(C.Get()))
+   {
+      MFEM_VERIFY(!(SolverC->iterative_mode),
+                  "Operator C of a TripleProductOperator should not be in iterative mode");
    }
 
    mem_class = A->GetMemoryClass()*C->GetMemoryClass();
@@ -500,12 +483,10 @@ TripleProductOperator::TripleProductOperator(
    t2.SetSize(B->Height(), mem_type);
 }
 
-TripleProductOperator::~TripleProductOperator()
-{
-   if (ownA) { delete A; }
-   if (ownB) { delete B; }
-   if (ownC) { delete C; }
-}
+TripleProductOperator::TripleProductOperator(
+   const Operator *A_, const Operator *B_, const Operator *C_,
+   bool own_A, bool own_B, bool own_C)
+   : TripleProductOperator({A_, own_A}, {B_, own_B}, {C_, own_C}) { }
 
 
 ConstrainedOperator::ConstrainedOperator(Operator *A, const Array<int> &list,

linalg/operator.hpp

@@ -13,6 +13,7 @@
 #define MFEM_OPERATOR
 
 #include "vector.hpp"
+#include "../general/handle.hpp"
 
 namespace mfem
 {
@@ -869,43 +870,42 @@ public:
 /// General linear combination operator: x -> a A(x) + b B(x).
 class SumOperator : public Operator
 {
-   const Operator *A, *B;
+   Handle<const Operator> A, B;
    const real_t alpha, beta;
-   bool ownA, ownB;
    mutable Vector z;
 
 public:
-   SumOperator(
-      const Operator *A, const real_t alpha,
-      const Operator *B, const real_t beta,
-      bool ownA, bool ownB);
+   SumOperator(Handle<const Operator> A_, const real_t alpha,
+               Handle<const Operator> B_, const real_t beta);
+
+   SumOperator(const Operator *A_, const real_t alpha,
+               const Operator *B_, const real_t beta,
+               bool own_A, bool own_B);
 
    void Mult(const Vector &x, Vector &y) const override
    { z.SetSize(A->Height()); A->Mult(x, z); B->Mult(x, y); add(alpha, z, beta, y, y); }
 
    void MultTranspose(const Vector &x, Vector &y) const override
    { z.SetSize(A->Width()); A->MultTranspose(x, z); B->MultTranspose(x, y); add(alpha, z, beta, y, y); }
-
-   virtual ~SumOperator();
 };
 
 /// General product operator: x -> (A*B)(x) = A(B(x)).
 class ProductOperator : public Operator
 {
-   const Operator *A, *B;
-   bool ownA, ownB;
+   Handle<const Operator> A, B;
    mutable Vector z;
 
 public:
-   ProductOperator(const Operator *A, const Operator *B, bool ownA, bool ownB);
+   ProductOperator(Handle<const Operator> A_, Handle<const Operator> B_);
+
+   ProductOperator(const Operator *A_, const Operator *B_,
+                   bool own_A, bool own_B);
 
    void Mult(const Vector &x, Vector &y) const override
    { B->Mult(x, z); A->Mult(z, y); }
 
    void MultTranspose(const Vector &x, Vector &y) const override
    { A->MultTranspose(x, z); B->MultTranspose(z, y); }
-
-   virtual ~ProductOperator();
 };
 
 
@@ -956,16 +956,19 @@ public:
 /// General triple product operator x -> A*B*C*x, with ownership of the factors.
 class TripleProductOperator : public Operator
 {
-   const Operator *A;
-   const Operator *B;
-   const Operator *C;
-   bool ownA, ownB, ownC;
+   Handle<const Operator> A;
+   Handle<const Operator> B;
+   Handle<const Operator> C;
    mutable Vector t1, t2;
    MemoryClass mem_class;
 
 public:
-   TripleProductOperator(const Operator *A, const Operator *B,
-                         const Operator *C, bool ownA, bool ownB, bool ownC);
+   TripleProductOperator(Handle<const Operator> A_, Handle<const Operator> B_,
+                         Handle<const Operator> C_);
+
+   TripleProductOperator(
+      const Operator *A_, const Operator *B_, const Operator *C_,
+      bool own_A, bool own_B, bool own_C);
 
    MemoryClass GetMemoryClass() const override { return mem_class; }
 
@@ -974,8 +977,6 @@ public:
 
    void MultTranspose(const Vector &x, Vector &y) const override
    { A->MultTranspose(x, t2); B->MultTranspose(t2, t1); C->MultTranspose(t1, y); }
-
-   virtual ~TripleProductOperator();
 };
 
 

linalg/petsc.hpp

@@ -21,7 +21,7 @@
 
 #include <limits>
 
-#include "handle.hpp"
+#include "op_handle.hpp"
 #include "hypre.hpp"
 #include "ode.hpp"
 #include "../general/mem_manager.hpp"

linalg/solvers.hpp

@@ -14,7 +14,7 @@
 
 #include "../config/config.hpp"
 #include "densemat.hpp"
-#include "handle.hpp"
+#include "op_handle.hpp"
 #include <memory>
 
 #ifdef MFEM_USE_MPI

tests/unit/linalg/test_operator.cpp

@@ -104,3 +104,30 @@ TEST_CASE("ConstrainedOperator", "[ConstrainedOperator][Operator]")
    REQUIRE(constrained_mult_application(A, list, x, y_true_zero_transpose, true,
                                         Operator::DiagonalPolicy::DIAG_ZERO) == MFEM_Approx(0.0));
 }
+
+TEST_CASE("Sum and product operators", "[Operator]")
+{
+   const int n = 1;
+   IdentityOperator op_1(n);
+   auto op_2 = MakeOwning<IdentityOperator>(n);
+
+   // op_1 will not be owned, op_2 will be owned
+   SumOperator sum(NonOwning(&op_1), 1.0, op_2, 1.0);
+   ProductOperator product(NonOwning(&op_1), op_2);
+   // Note: it is not a problem for triple to own op_2 'twice'
+   TripleProductOperator triple(NonOwning(&op_1), op_2, op_2);
+
+   // Even though op_2 is reset here, it remains valid in each of the operators
+   op_2.Reset();
+
+   Vector x(n), y(n);
+   x = 1.0;
+   sum.Mult(x, y);
+   REQUIRE(y[0] == 2.0);
+
+   product.Mult(x, y);
+   REQUIRE(y[0] == 1.0);
+
+   triple.Mult(x, y);
+   REQUIRE(y[0] == 1.0);
+}