Added instantiation of TBlockOperator to compile it in blockoperator.cpp.

CHANGELOG

@@ -50,6 +50,8 @@ Miscellaneous
   IntegrationRules IntRules, RefinedIntRules, GeometryRefiner
   GlobGeometryRefiner, and FiniteElement::dof2quad_array.
 
+- Added TBlockOperator for generic handling of block-structured operators.
+
 
 Version 4.6, released on September 27, 2023
 ===========================================

linalg/CMakeLists.txt

@@ -51,6 +51,7 @@ list(APPEND HDRS
   solvers.hpp
   sparsemat.hpp
   sparsesmoothers.hpp
+  tblockoperator.hpp
   tlayout.hpp
   tmatrix.hpp
   ttensor.hpp

linalg/blockoperator.cpp

@@ -18,129 +18,9 @@
 namespace mfem
 {
 
-BlockOperator::BlockOperator(const Array<int> & offsets)
-   : Operator(offsets.Last()),
-     owns_blocks(0),
-     nRowBlocks(offsets.Size() - 1),
-     nColBlocks(offsets.Size() - 1),
-     row_offsets(offsets),
-     col_offsets(offsets),
-     op(nRowBlocks, nRowBlocks),
-     coef(nRowBlocks, nColBlocks)
-{
-   op = static_cast<Operator *>(NULL);
-}
+template class TBlockOperator<const Operator>;
+template class TBlockOperator<Operator>;
 
-BlockOperator::BlockOperator(const Array<int> & row_offsets_,
-                             const Array<int> & col_offsets_)
-   : Operator(row_offsets_.Last(), col_offsets_.Last()),
-     owns_blocks(0),
-     nRowBlocks(row_offsets_.Size()-1),
-     nColBlocks(col_offsets_.Size()-1),
-     row_offsets(row_offsets_),
-     col_offsets(col_offsets_),
-     op(nRowBlocks, nColBlocks),
-     coef(nRowBlocks, nColBlocks)
-{
-   op = static_cast<Operator *>(NULL);
-}
-
-void BlockOperator::SetDiagonalBlock(int iblock, const Operator *opt, double c)
-{
-   SetBlock(iblock, iblock, opt, c);
-}
-
-void BlockOperator::SetBlock(int iRow, int iCol, const Operator *opt, double c)
-{
-   if (owns_blocks && op(iRow, iCol))
-   {
-      delete op(iRow, iCol);
-   }
-   op(iRow, iCol) = opt;
-   coef(iRow, iCol) = c;
-
-   MFEM_VERIFY(row_offsets[iRow+1] - row_offsets[iRow] == opt->NumRows() &&
-               col_offsets[iCol+1] - col_offsets[iCol] == opt->NumCols(),
-               "incompatible Operator dimensions");
-}
-
-// Operator application
-void BlockOperator::Mult (const Vector & x, Vector & y) const
-{
-   MFEM_ASSERT(x.Size() == width, "incorrect input Vector size");
-   MFEM_ASSERT(y.Size() == height, "incorrect output Vector size");
-
-   x.Read();
-   y.Write(); y = 0.0;
-
-   xblock.Update(const_cast<Vector&>(x),col_offsets);
-   yblock.Update(y,row_offsets);
-
-   for (int iRow=0; iRow < nRowBlocks; ++iRow)
-   {
-      tmp.SetSize(row_offsets[iRow+1] - row_offsets[iRow]);
-      for (int jCol=0; jCol < nColBlocks; ++jCol)
-      {
-         if (op(iRow,jCol))
-         {
-            op(iRow,jCol)->Mult(xblock.GetBlock(jCol), tmp);
-            yblock.GetBlock(iRow).Add(coef(iRow,jCol), tmp);
-         }
-      }
-   }
-
-   for (int iRow=0; iRow < nRowBlocks; ++iRow)
-   {
-      yblock.GetBlock(iRow).SyncAliasMemory(y);
-   }
-}
-
-// Action of the transpose operator
-void BlockOperator::MultTranspose (const Vector & x, Vector & y) const
-{
-   MFEM_ASSERT(x.Size() == height, "incorrect input Vector size");
-   MFEM_ASSERT(y.Size() == width, "incorrect output Vector size");
-
-   x.Read();
-   y.Write(); y = 0.0;
-
-   xblock.Update(const_cast<Vector&>(x),row_offsets);
-   yblock.Update(y,col_offsets);
-
-   for (int iRow=0; iRow < nColBlocks; ++iRow)
-   {
-      tmp.SetSize(col_offsets[iRow+1] - col_offsets[iRow]);
-      for (int jCol=0; jCol < nRowBlocks; ++jCol)
-      {
-         if (op(jCol,iRow))
-         {
-            op(jCol,iRow)->MultTranspose(xblock.GetBlock(jCol), tmp);
-            yblock.GetBlock(iRow).Add(coef(jCol,iRow), tmp);
-         }
-      }
-   }
-
-   for (int iRow=0; iRow < nColBlocks; ++iRow)
-   {
-      yblock.GetBlock(iRow).SyncAliasMemory(y);
-   }
-}
-
-BlockOperator::~BlockOperator()
-{
-   if (owns_blocks)
-   {
-      for (int iRow=0; iRow < nRowBlocks; ++iRow)
-      {
-         for (int jCol=0; jCol < nColBlocks; ++jCol)
-         {
-            delete op(iRow,jCol);
-         }
-      }
-   }
-}
-
-//-----------------------------------------------------------------------
 BlockDiagonalPreconditioner::BlockDiagonalPreconditioner(
    const Array<int> & offsets_):
    Solver(offsets_.Last()),

linalg/blockoperator.hpp

@@ -14,13 +14,29 @@
 
 #include "../config/config.hpp"
 #include "../general/array.hpp"
+#include "tblockoperator.hpp"
 #include "operator.hpp"
 #include "blockvector.hpp"
 
 namespace mfem
 {
 
-//! @class BlockOperator
+//! @typedef BlockConstOperator
+/**
+ * \brief A class to handle constant Block systems in a matrix-free implementation.
+ *
+ * Usage:
+ * - Use one of the constructors to define the block structure.
+ * - Use SetDiagonalBlock or SetBlock to fill the BlockConstOperator
+ * - Use the method Mult and MultTranspose to apply the operator to a vector.
+ *
+ * If a block is not set, it is assumed to be a zero block.
+ */
+typedef TBlockOperator<const Operator> BlockConstOperator;
+
+extern template class TBlockOperator<const Operator>;
+
+//! @typedef BlockOperator
 /**
  * \brief A class to handle Block systems in a matrix-free implementation.
  *
@@ -31,104 +47,9 @@ namespace mfem
  *
  * If a block is not set, it is assumed to be a zero block.
  */
-class BlockOperator : public Operator
-{
-public:
-   //! Constructor for BlockOperators with the same block-structure for rows and
-   //! columns.
-   /**
-    *  offsets: offsets that mark the start of each row/column block (size
-    *  nRowBlocks+1).
-    */
-   BlockOperator(const Array<int> & offsets);
-   //! Constructor for general BlockOperators.
-   /**
-    *  row_offsets: offsets that mark the start of each row block (size
-    *  nRowBlocks+1).  col_offsets: offsets that mark the start of each column
-    *  block (size nColBlocks+1).
-    */
-   BlockOperator(const Array<int> & row_offsets, const Array<int> & col_offsets);
+typedef TBlockOperator<Operator> BlockOperator;
 
-   /// Copy assignment is not supported
-   BlockOperator &operator=(const BlockOperator &) = delete;
-
-   /// Move assignment is not supported
-   BlockOperator &operator=(BlockOperator &&) = delete;
-
-   //! Add block op in the block-entry (iblock, iblock).
-   /**
-    * iblock: The block will be inserted in location (iblock, iblock).
-    * op: the Operator to be inserted.
-    * c: optional scalar multiple for this block.
-    */
-   void SetDiagonalBlock(int iblock, const Operator *op, double c = 1.0);
-   //! Add a block op in the block-entry (iblock, jblock).
-   /**
-    * irow, icol: The block will be inserted in location (irow, icol).
-    * op: the Operator to be inserted.
-    * c: optional scalar multiple for this block.
-    */
-   void SetBlock(int iRow, int iCol, const Operator *op, double c = 1.0);
-
-   //! Return the number of row blocks
-   int NumRowBlocks() const { return nRowBlocks; }
-   //! Return the number of column blocks
-   int NumColBlocks() const { return nColBlocks; }
-
-   //! Check if block (i,j) is a zero block
-   int IsZeroBlock(int i, int j) const { return (op(i,j)==NULL) ? 1 : 0; }
-   //! Return a reference to block i,j
-   const Operator & GetBlock(int i, int j) const
-   { MFEM_VERIFY(op(i,j), ""); return *op(i,j); }
-   //! Return the coefficient for block i,j
-   double GetBlockCoef(int i, int j) const
-   { MFEM_VERIFY(op(i,j), ""); return coef(i,j); }
-   //! Set the coefficient for block i,j
-   void SetBlockCoef(int i, int j, double c)
-   { MFEM_VERIFY(op(i,j), ""); coef(i,j) = c; }
-
-   //! Return the row offsets for block starts
-   Array<int> & RowOffsets() { return row_offsets; }
-   //! Read only access to the row offsets for block starts
-   const Array<int> & RowOffsets() const { return row_offsets; }
-   //! Return the columns offsets for block starts
-   Array<int> & ColOffsets() { return col_offsets; }
-   //! Read only access to the columns offsets for block starts
-   const Array<int> & ColOffsets() const { return col_offsets; }
-
-   /// Operator application
-   virtual void Mult (const Vector & x, Vector & y) const;
-
-   /// Action of the transpose operator
-   virtual void MultTranspose (const Vector & x, Vector & y) const;
-
-   ~BlockOperator();
-
-   //! Controls the ownership of the blocks: if nonzero, BlockOperator will
-   //! delete all blocks that are set (non-NULL); the default value is zero.
-   int owns_blocks;
-
-   virtual Type GetType() const { return MFEM_Block_Operator; }
-
-private:
-   //! Number of block rows
-   int nRowBlocks;
-   //! Number of block columns
-   int nColBlocks;
-   //! Row offsets for the starting position of each block
-   Array<int> row_offsets;
-   //! Column offsets for the starting position of each block
-   Array<int> col_offsets;
-   //! 2D array that stores each block of the operator.
-   Array2D<const Operator *> op;
-   //! 2D array that stores a coefficient for each block of the operator.
-   Array2D<double> coef;
-
-   //! Temporary Vectors used to efficiently apply the Mult and MultTranspose methods.
-   mutable BlockVector xblock;
-   mutable BlockVector yblock;
-   mutable Vector tmp;
-};
+extern template class TBlockOperator<Operator>;
 
 //! @class BlockDiagonalPreconditioner
 /**
@@ -136,7 +57,7 @@ private:
  *
  * Usage:
  * - Use the constructors to define the block structure
- * - Use SetDiagonalBlock to fill the BlockOperator
+ * - Use SetDiagonalBlock to fill the BlockDiagonalPreconditioner
  * - Use the method Mult and MultTranspose to apply the operator to a vector.
  *
  * If a block is not set, it is assumed to be an identity block.
@@ -206,7 +127,7 @@ private:
  *
  * Usage:
  * - Use the constructors to define the block structure
- * - Use SetBlock() to fill the BlockOperator
+ * - Use SetBlock() to fill the BlockLowerTriangularOperator
  * - Diagonal blocks of the preconditioner should approximate the inverses of
  *   the diagonal block of the matrix
  * - Off-diagonal blocks of the preconditioner should match/approximate those of
@@ -220,7 +141,7 @@ private:
 class BlockLowerTriangularPreconditioner : public Solver
 {
 public:
-   //! Constructor for BlockLowerTriangularPreconditioners with the same
+   //! Constructor for BlockLowerTriangularPreconditioner%s with the same
    //! block-structure for rows and columns.
    /**
     *  @param offsets  Offsets that mark the start of each row/column block

linalg/tblockoperator.hpp

@@ -0,0 +1,272 @@
+// 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_TBLOCKOPERATOR
+#define MFEM_TBLOCKOPERATOR
+
+#include "../config/config.hpp"
+#include "../general/array.hpp"
+#include "operator.hpp"
+#include "blockvector.hpp"
+
+namespace mfem
+{
+
+//! @class TBlockOperator
+/**
+ * \brief A template to handle Block systems in a matrix-free implementation.
+ *
+ * Usage:
+ * - Use one of the constructors to define the block structure.
+ * - Use SetDiagonalBlock or SetBlock to fill the TBlockOperator
+ * - Use the method Mult and MultTranspose to apply the operator to a vector.
+ *
+ * If a block is not set, it is assumed to be a zero block.
+ */
+template<class T> class TBlockOperator : public Operator
+{
+public:
+   //! Constructor for TBlockOperator%s with the same block-structure for rows and
+   //! columns.
+   /**
+    *  offsets: offsets that mark the start of each row/column block (size
+    *  nRowBlocks+1).
+    */
+   TBlockOperator(const Array<int> & offsets);
+   //! Constructor for general TBlockOperator%s.
+   /**
+    *  row_offsets: offsets that mark the start of each row block (size
+    *  nRowBlocks+1).  col_offsets: offsets that mark the start of each column
+    *  block (size nColBlocks+1).
+    */
+   TBlockOperator(const Array<int> & row_offsets,
+                  const Array<int> & col_offsets);
+
+   /// Copy assignment is not supported
+   TBlockOperator &operator=(const TBlockOperator &) = delete;
+
+   /// Move assignment is not supported
+   TBlockOperator &operator=(TBlockOperator &&) = delete;
+
+   //! Add block op in the block-entry (iblock, iblock).
+   /**
+    * iblock: The block will be inserted in location (iblock, iblock).
+    * op: the Operator to be inserted.
+    * c: optional scalar multiple for this block.
+    */
+   void SetDiagonalBlock(int iblock, T *op, double c = 1.0);
+   //! Add a block op in the block-entry (iblock, jblock).
+   /**
+    * irow, icol: The block will be inserted in location (irow, icol).
+    * op: the Operator to be inserted.
+    * c: optional scalar multiple for this block.
+    */
+   void SetBlock(int iRow, int iCol, T *op, double c = 1.0);
+
+   //! Return the number of row blocks
+   int NumRowBlocks() const { return nRowBlocks; }
+   //! Return the number of column blocks
+   int NumColBlocks() const { return nColBlocks; }
+
+   //! Check if block (i,j) is a zero block
+   int IsZeroBlock(int i, int j) const { return (op(i,j)==NULL) ? 1 : 0; }
+   //! Return a reference to block i,j
+   T& GetBlock(int i, int j)
+   { MFEM_VERIFY(op(i,j), ""); return *op(i,j); }
+   //! Return a reference to block i,j
+   const T& GetBlock(int i, int j) const
+   { MFEM_VERIFY(op(i,j), ""); return *op(i,j); }
+   //! Return the coefficient for block i,j
+   double GetBlockCoef(int i, int j) const
+   { MFEM_VERIFY(op(i,j), ""); return coef(i,j); }
+   //! Set the coefficient for block i,j
+   void SetBlockCoef(int i, int j, double c)
+   { MFEM_VERIFY(op(i,j), ""); coef(i,j) = c; }
+
+   //! Return the row offsets for block starts
+   Array<int> & RowOffsets() { return row_offsets; }
+   //! Read only access to the row offsets for block starts
+   const Array<int> & RowOffsets() const { return row_offsets; }
+   //! Return the columns offsets for block starts
+   Array<int> & ColOffsets() { return col_offsets; }
+   //! Read only access to the columns offsets for block starts
+   const Array<int> & ColOffsets() const { return col_offsets; }
+
+   /// Operator application
+   virtual void Mult (const Vector & x, Vector & y) const;
+
+   /// Action of the transpose operator
+   virtual void MultTranspose (const Vector & x, Vector & y) const;
+
+   virtual ~TBlockOperator();
+
+   //! Controls the ownership of the blocks: if nonzero, TBlockOperator will
+   //! delete all blocks that are set (non-NULL); the default value is zero.
+   int owns_blocks;
+
+   virtual Type GetType() const { return MFEM_Block_Operator; }
+
+protected:
+   //! Number of block rows
+   int nRowBlocks;
+   //! Number of block columns
+   int nColBlocks;
+   //! Row offsets for the starting position of each block
+   Array<int> row_offsets;
+   //! Column offsets for the starting position of each block
+   Array<int> col_offsets;
+   //! 2D array that stores each block of the operator.
+   Array2D<T *> op;
+   //! 2D array that stores a coefficient for each block of the operator.
+   Array2D<double> coef;
+
+   //! Temporary Vectors used to efficiently apply the Mult and MultTranspose methods.
+   mutable BlockVector xblock;
+   mutable BlockVector yblock;
+   mutable Vector tmp;
+};
+
+//-----------------------------------------------------------------------
+
+template<class T>
+TBlockOperator<T>::TBlockOperator(const Array<int> & offsets)
+   : Operator(offsets.Last()),
+     owns_blocks(0),
+     nRowBlocks(offsets.Size() - 1),
+     nColBlocks(offsets.Size() - 1),
+     row_offsets(offsets),
+     col_offsets(offsets),
+     op(nRowBlocks, nRowBlocks),
+     coef(nRowBlocks, nColBlocks)
+{
+   op = static_cast<T*>(NULL);
+}
+
+template<class T>
+TBlockOperator<T>::TBlockOperator(const Array<int> & row_offsets_,
+                                  const Array<int> & col_offsets_)
+   : Operator(row_offsets_.Last(), col_offsets_.Last()),
+     owns_blocks(0),
+     nRowBlocks(row_offsets_.Size()-1),
+     nColBlocks(col_offsets_.Size()-1),
+     row_offsets(row_offsets_),
+     col_offsets(col_offsets_),
+     op(nRowBlocks, nColBlocks),
+     coef(nRowBlocks, nColBlocks)
+{
+   op = static_cast<T*>(NULL);
+}
+
+template<class T>
+void TBlockOperator<T>::SetDiagonalBlock(int iblock, T *opt,
+                                         double c)
+{
+   SetBlock(iblock, iblock, opt, c);
+}
+
+template<class T>
+void TBlockOperator<T>::SetBlock(int iRow, int iCol, T *opt,
+                                 double c)
+{
+   if (owns_blocks && op(iRow, iCol))
+   {
+      delete op(iRow, iCol);
+   }
+   op(iRow, iCol) = opt;
+   coef(iRow, iCol) = c;
+
+   MFEM_VERIFY(row_offsets[iRow+1] - row_offsets[iRow] == opt->NumRows() &&
+               col_offsets[iCol+1] - col_offsets[iCol] == opt->NumCols(),
+               "incompatible Operator dimensions");
+}
+
+// Operator application
+template<class T>
+void TBlockOperator<T>::Mult (const Vector & x, Vector & y) const
+{
+   MFEM_ASSERT(x.Size() == width, "incorrect input Vector size");
+   MFEM_ASSERT(y.Size() == height, "incorrect output Vector size");
+
+   x.Read();
+   y.Write(); y = 0.0;
+
+   xblock.Update(const_cast<Vector&>(x),col_offsets);
+   yblock.Update(y,row_offsets);
+
+   for (int iRow=0; iRow < nRowBlocks; ++iRow)
+   {
+      tmp.SetSize(row_offsets[iRow+1] - row_offsets[iRow]);
+      for (int jCol=0; jCol < nColBlocks; ++jCol)
+      {
+         if (op(iRow,jCol))
+         {
+            op(iRow,jCol)->Mult(xblock.GetBlock(jCol), tmp);
+            yblock.GetBlock(iRow).Add(coef(iRow,jCol), tmp);
+         }
+      }
+   }
+
+   for (int iRow=0; iRow < nRowBlocks; ++iRow)
+   {
+      yblock.GetBlock(iRow).SyncAliasMemory(y);
+   }
+}
+
+// Action of the transpose operator
+template<class T>
+void TBlockOperator<T>::MultTranspose (const Vector & x, Vector & y) const
+{
+   MFEM_ASSERT(x.Size() == height, "incorrect input Vector size");
+   MFEM_ASSERT(y.Size() == width, "incorrect output Vector size");
+
+   x.Read();
+   y.Write(); y = 0.0;
+
+   xblock.Update(const_cast<Vector&>(x),row_offsets);
+   yblock.Update(y,col_offsets);
+
+   for (int iRow=0; iRow < nColBlocks; ++iRow)
+   {
+      tmp.SetSize(col_offsets[iRow+1] - col_offsets[iRow]);
+      for (int jCol=0; jCol < nRowBlocks; ++jCol)
+      {
+         if (op(jCol,iRow))
+         {
+            op(jCol,iRow)->MultTranspose(xblock.GetBlock(jCol), tmp);
+            yblock.GetBlock(iRow).Add(coef(jCol,iRow), tmp);
+         }
+      }
+   }
+
+   for (int iRow=0; iRow < nColBlocks; ++iRow)
+   {
+      yblock.GetBlock(iRow).SyncAliasMemory(y);
+   }
+}
+
+template<class T>
+TBlockOperator<T>::~TBlockOperator()
+{
+   if (owns_blocks)
+   {
+      for (int iRow=0; iRow < nRowBlocks; ++iRow)
+      {
+         for (int jCol=0; jCol < nColBlocks; ++jCol)
+         {
+            delete op(iRow,jCol);
+         }
+      }
+   }
+}
+
+}
+
+#endif /* MFEM_TBLOCKOPERATOR */

miniapps/dpg/pacoustics.cpp

@@ -616,8 +616,8 @@ int main(int argc, char *argv[])
       a->FormLinearSystem(ess_tdof_list,x,Ah, X,B);
 
       ComplexOperator * Ahc = Ah.As<ComplexOperator>();
-      BlockOperator * BlockA_r = dynamic_cast<BlockOperator *>(&Ahc->real());
-      BlockOperator * BlockA_i = dynamic_cast<BlockOperator *>(&Ahc->imag());
+      auto * BlockA_r = dynamic_cast<TBlockOperator<HypreParMatrix> *>(&Ahc->real());
+      auto * BlockA_i = dynamic_cast<TBlockOperator<HypreParMatrix> *>(&Ahc->imag());
 
       int num_blocks = BlockA_r->NumRowBlocks();
       Array<int> tdof_offsets(2*num_blocks+1);
@@ -650,12 +650,10 @@ int main(int argc, char *argv[])
 
       if (!static_cond)
       {
-         HypreBoomerAMG * solver_p = new HypreBoomerAMG((HypreParMatrix &)
-                                                        BlockA_r->GetBlock(0,0));
+         HypreBoomerAMG * solver_p = new HypreBoomerAMG(BlockA_r->GetBlock(0,0));
          solver_p->SetPrintLevel(0);
          solver_p->SetSystemsOptions(dim);
-         HypreBoomerAMG * solver_u = new HypreBoomerAMG((HypreParMatrix &)
-                                                        BlockA_r->GetBlock(1,1));
+         HypreBoomerAMG * solver_u = new HypreBoomerAMG(BlockA_r->GetBlock(1,1));
          solver_u->SetPrintLevel(0);
          solver_u->SetSystemsOptions(dim);
          M.SetDiagonalBlock(0,solver_p);
@@ -664,22 +662,22 @@ int main(int argc, char *argv[])
          M.SetDiagonalBlock(num_blocks+1,solver_u);
       }
 
-      HypreBoomerAMG * solver_hatp = new HypreBoomerAMG((HypreParMatrix &)
-                                                        BlockA_r->GetBlock(skip,skip));
+      HypreBoomerAMG * solver_hatp = new HypreBoomerAMG(BlockA_r->GetBlock(skip,
+                                                                           skip));
       solver_hatp->SetPrintLevel(0);
 
       HypreSolver * solver_hatu = nullptr;
       if (dim == 2)
       {
          // AMS preconditioner for 2D H(div) (trace) space
-         solver_hatu = new HypreAMS((HypreParMatrix &)BlockA_r->GetBlock(skip+1,skip+1),
+         solver_hatu = new HypreAMS(BlockA_r->GetBlock(skip+1,skip+1),
                                     hatu_fes);
          dynamic_cast<HypreAMS*>(solver_hatu)->SetPrintLevel(0);
       }
       else
       {
          // ADS preconditioner for 3D H(div) (trace) space
-         solver_hatu = new HypreADS((HypreParMatrix &)BlockA_r->GetBlock(skip+1,skip+1),
+         solver_hatu = new HypreADS(BlockA_r->GetBlock(skip+1,skip+1),
                                     hatu_fes);
          dynamic_cast<HypreADS*>(solver_hatu)->SetPrintLevel(0);
       }

miniapps/dpg/pconvection-diffusion.cpp

@@ -449,23 +449,22 @@ int main(int argc, char *argv[])
       Vector X,B;
       a->FormLinearSystem(ess_tdof_list,x,Ah,X,B);
 
-      BlockOperator * A = Ah.As<BlockOperator>();
+      auto * A = Ah.As<TBlockOperator<HypreParMatrix>>();
 
       BlockDiagonalPreconditioner M(A->RowOffsets());
       M.owns_blocks = 1;
       int skip = 0;
       if (!static_cond)
       {
-         HypreBoomerAMG * amg0 = new HypreBoomerAMG((HypreParMatrix &)A->GetBlock(0,0));
-         HypreBoomerAMG * amg1 = new HypreBoomerAMG((HypreParMatrix &)A->GetBlock(1,1));
+         HypreBoomerAMG * amg0 = new HypreBoomerAMG(A->GetBlock(0,0));
+         HypreBoomerAMG * amg1 = new HypreBoomerAMG(A->GetBlock(1,1));
          amg0->SetPrintLevel(0);
          amg1->SetPrintLevel(0);
          M.SetDiagonalBlock(0,amg0);
          M.SetDiagonalBlock(1,amg1);
          skip = 2;
       }
-      HypreBoomerAMG * amg2 = new HypreBoomerAMG((HypreParMatrix &)A->GetBlock(skip,
-                                                                               skip));
+      HypreBoomerAMG * amg2 = new HypreBoomerAMG(A->GetBlock(skip, skip));
       amg2->SetPrintLevel(0);
       M.SetDiagonalBlock(skip,amg2);
 
@@ -473,7 +472,7 @@ int main(int argc, char *argv[])
       if (dim == 2)
       {
          // AMS preconditioner for 2D H(div) (trace) space
-         prec = new HypreAMS((HypreParMatrix &)A->GetBlock(skip+1,skip+1), hatf_fes);
+         prec = new HypreAMS(A->GetBlock(skip+1,skip+1), hatf_fes);
       }
       else
       {

miniapps/dpg/pdiffusion.cpp

@@ -374,35 +374,34 @@ int main(int argc, char *argv[])
       OperatorPtr Ah;
       a->FormLinearSystem(ess_tdof_list,x,Ah,X,B);
 
-      BlockOperator * A = Ah.As<BlockOperator>();
+      auto * A = Ah.As<TBlockOperator<HypreParMatrix>>();
 
       BlockDiagonalPreconditioner M(A->RowOffsets());
       M.owns_blocks = 1;
       int skip = 0;
       if (!static_cond)
       {
-         HypreBoomerAMG * amg0 = new HypreBoomerAMG((HypreParMatrix &)A->GetBlock(0,0));
-         HypreBoomerAMG * amg1 = new HypreBoomerAMG((HypreParMatrix &)A->GetBlock(1,1));
+         HypreBoomerAMG * amg0 = new HypreBoomerAMG(A->GetBlock(0,0));
+         HypreBoomerAMG * amg1 = new HypreBoomerAMG(A->GetBlock(1,1));
          amg0->SetPrintLevel(0);
          amg1->SetPrintLevel(0);
          M.SetDiagonalBlock(0,amg0);
          M.SetDiagonalBlock(1,amg1);
          skip=2;
       }
-      HypreBoomerAMG * amg2 = new HypreBoomerAMG((HypreParMatrix &)A->GetBlock(skip,
-                                                                               skip));
+      HypreBoomerAMG * amg2 = new HypreBoomerAMG(A->GetBlock(skip, skip));
       amg2->SetPrintLevel(0);
       M.SetDiagonalBlock(skip,amg2);
       HypreSolver * prec;
       if (dim == 2)
       {
          // AMS preconditioner for 2D H(div) (trace) space
-         prec = new HypreAMS((HypreParMatrix &)A->GetBlock(skip+1,skip+1), hatsigma_fes);
+         prec = new HypreAMS(A->GetBlock(skip+1,skip+1), hatsigma_fes);
       }
       else
       {
          // ADS preconditioner for 3D H(div) (trace) space
-         prec = new HypreADS((HypreParMatrix &)A->GetBlock(skip+1,skip+1), hatsigma_fes);
+         prec = new HypreADS(A->GetBlock(skip+1,skip+1), hatsigma_fes);
       }
       M.SetDiagonalBlock(skip+1,prec);
 

miniapps/dpg/pmaxwell.cpp

@@ -820,8 +820,8 @@ int main(int argc, char *argv[])
 
       ComplexOperator * Ahc = Ah.As<ComplexOperator>();
 
-      BlockOperator * BlockA_r = dynamic_cast<BlockOperator *>(&Ahc->real());
-      BlockOperator * BlockA_i = dynamic_cast<BlockOperator *>(&Ahc->imag());
+      auto * BlockA_r = dynamic_cast<TBlockOperator<HypreParMatrix>*>(&Ahc->real());
+      auto * BlockA_i = dynamic_cast<TBlockOperator<HypreParMatrix>*>(&Ahc->imag());
 
       int num_blocks = BlockA_r->NumRowBlocks();
       Array<int> tdof_offsets(2*num_blocks+1);
@@ -853,12 +853,10 @@ int main(int argc, char *argv[])
 
       if (!static_cond)
       {
-         HypreBoomerAMG * solver_E = new HypreBoomerAMG((HypreParMatrix &)
-                                                        BlockA_r->GetBlock(0,0));
+         HypreBoomerAMG * solver_E = new HypreBoomerAMG(BlockA_r->GetBlock(0,0));
          solver_E->SetPrintLevel(0);
          solver_E->SetSystemsOptions(dim);
-         HypreBoomerAMG * solver_H = new HypreBoomerAMG((HypreParMatrix &)
-                                                        BlockA_r->GetBlock(1,1));
+         HypreBoomerAMG * solver_H = new HypreBoomerAMG(BlockA_r->GetBlock(1,1));
          solver_H->SetPrintLevel(0);
          solver_H->SetSystemsOptions(dim);
          M.SetDiagonalBlock(0,solver_E);
@@ -868,20 +866,17 @@ int main(int argc, char *argv[])
       }
 
       HypreSolver * solver_hatH = nullptr;
-      HypreAMS * solver_hatE = new HypreAMS((HypreParMatrix &)BlockA_r->GetBlock(skip,
-                                                                                 skip),
+      HypreAMS * solver_hatE = new HypreAMS(BlockA_r->GetBlock(skip, skip),
                                             hatE_fes);
       solver_hatE->SetPrintLevel(0);
       if (dim == 2)
       {
-         solver_hatH = new HypreBoomerAMG((HypreParMatrix &)BlockA_r->GetBlock(skip+1,
-                                                                               skip+1));
+         solver_hatH = new HypreBoomerAMG(BlockA_r->GetBlock(skip+1, skip+1));
          dynamic_cast<HypreBoomerAMG*>(solver_hatH)->SetPrintLevel(0);
       }
       else
       {
-         solver_hatH = new HypreAMS((HypreParMatrix &)BlockA_r->GetBlock(skip+1,skip+1),
-                                    hatH_fes);
+         solver_hatH = new HypreAMS(BlockA_r->GetBlock(skip+1,skip+1), hatH_fes);
          dynamic_cast<HypreAMS*>(solver_hatH)->SetPrintLevel(0);
       }
 

miniapps/dpg/util/blockstaticcond.cpp

@@ -547,7 +547,7 @@ void BlockStaticCondensation::BuildProlongation()
 void BlockStaticCondensation::BuildParallelProlongation()
 {
    MFEM_VERIFY(parallel, "BuildParallelProlongation: wrong code path");
-   pP = new BlockOperator(rdof_offsets, rtdof_offsets);
+   pP = new TBlockOperator<HypreParMatrix>(rdof_offsets, rtdof_offsets);
    R = new BlockMatrix(rtdof_offsets, rdof_offsets);
    pP->owns_blocks = 0;
    R->owns_blocks = 0;
@@ -571,8 +571,8 @@ void BlockStaticCondensation::ParallelAssemble(BlockMatrix *m)
 {
    if (!pP) { BuildParallelProlongation(); }
 
-   pS = new BlockOperator(rtdof_offsets);
-   pS_e = new BlockOperator(rtdof_offsets);
+   pS = new TBlockOperator<HypreParMatrix>(rtdof_offsets);
+   pS_e = new TBlockOperator<HypreParMatrix>(rtdof_offsets);
    pS->owns_blocks = 1;
    pS_e->owns_blocks = 1;
    HypreParMatrix * A = nullptr;
@@ -584,7 +584,7 @@ void BlockStaticCondensation::ParallelAssemble(BlockMatrix *m)
    {
       if (!tr_fes[i]) { continue; }
       pfes_i = dynamic_cast<ParFiniteElementSpace*>(tr_fes[i]);
-      HypreParMatrix * Pi = (HypreParMatrix*)(&pP->GetBlock(skip_i,skip_i));
+      HypreParMatrix * Pi = &pP->GetBlock(skip_i,skip_i);
       int skip_j=0;
       for (int j = 0; j<nblocks; j++)
       {
@@ -602,7 +602,7 @@ void BlockStaticCondensation::ParallelAssemble(BlockMatrix *m)
          else
          {
             pfes_j = dynamic_cast<ParFiniteElementSpace*>(tr_fes[j]);
-            HypreParMatrix * Pj = (HypreParMatrix*)(&pP->GetBlock(skip_j,skip_j));
+            HypreParMatrix * Pj = &pP->GetBlock(skip_j,skip_j);
             A = new HypreParMatrix(pfes_i->GetComm(), pfes_i->GlobalVSize(),
                                    pfes_j->GlobalVSize(), pfes_i->GetDofOffsets(),
                                    pfes_j->GetDofOffsets(), &m->GetBlock(skip_i,skip_j));
@@ -865,7 +865,7 @@ void BlockStaticCondensation::ReduceSystem(Vector &x, Vector &X,
       for (int j = 0; j<rblocks; j++)
       {
          if (!ess_tdofs[j]->Size()) { continue; }
-         HypreParMatrix *Ah = (HypreParMatrix *)(&pS->GetBlock(j,j));
+         HypreParMatrix *Ah = &pS->GetBlock(j,j);
          Vector diag;
          Ah->GetDiag(diag);
          for (int i = 0; i < ess_tdofs[j]->Size(); i++)

miniapps/dpg/util/blockstaticcond.hpp

@@ -65,9 +65,9 @@ class BlockStaticCondensation
    BlockMatrix * R = nullptr; // Block Restriction
 
 #ifdef MFEM_USE_MPI
-   BlockOperator * pS = nullptr;
-   BlockOperator * pS_e = nullptr;
-   BlockOperator * pP = nullptr;
+   TBlockOperator<HypreParMatrix> * pS = nullptr;
+   TBlockOperator<HypreParMatrix> * pS_e = nullptr;
+   TBlockOperator<HypreParMatrix> * pP = nullptr;
 #endif
 
    bool Parallel() const { return parallel; }
@@ -152,10 +152,10 @@ public:
 
 #ifdef MFEM_USE_MPI
    /// Return the parallel Schur complement matrix.
-   BlockOperator &GetParallelSchurMatrix() { return *pS; }
+   TBlockOperator<HypreParMatrix> &GetParallelSchurMatrix() { return *pS; }
 
    /// Return the eliminated part of the parallel Schur complement matrix.
-   BlockOperator &GetParallelSchurMatrixElim() { return *pS_e; }
+   TBlockOperator<HypreParMatrix> &GetParallelSchurMatrixElim() { return *pS_e; }
 
    void ParallelAssemble(BlockMatrix *m);
 #endif

miniapps/dpg/util/complexstaticcond.cpp

@@ -592,7 +592,7 @@ void ComplexBlockStaticCondensation::BuildProlongation()
 void ComplexBlockStaticCondensation::BuildParallelProlongation()
 {
    MFEM_VERIFY(parallel, "BuildParallelProlongation: wrong code path");
-   pP = new BlockOperator(rdof_offsets, rtdof_offsets);
+   pP = new TBlockOperator<HypreParMatrix>(rdof_offsets, rtdof_offsets);
    R = new BlockMatrix(rtdof_offsets, rdof_offsets);
    pP->owns_blocks = 0;
    R->owns_blocks = 0;
@@ -618,10 +618,10 @@ void ComplexBlockStaticCondensation::ParallelAssemble(BlockMatrix *m_r,
 
    if (!pP) { BuildParallelProlongation(); }
 
-   pS_r = new BlockOperator(rtdof_offsets);
-   pS_e_r = new BlockOperator(rtdof_offsets);
-   pS_i = new BlockOperator(rtdof_offsets);
-   pS_e_i = new BlockOperator(rtdof_offsets);
+   pS_r = new TBlockOperator<HypreParMatrix>(rtdof_offsets);
+   pS_e_r = new TBlockOperator<HypreParMatrix>(rtdof_offsets);
+   pS_i = new TBlockOperator<HypreParMatrix>(rtdof_offsets);
+   pS_e_i = new TBlockOperator<HypreParMatrix>(rtdof_offsets);
    pS_r->owns_blocks = 1;
    pS_i->owns_blocks = 1;
    pS_e_r->owns_blocks = 1;
@@ -637,7 +637,7 @@ void ComplexBlockStaticCondensation::ParallelAssemble(BlockMatrix *m_r,
    {
       if (!tr_fes[i]) { continue; }
       pfes_i = dynamic_cast<ParFiniteElementSpace*>(tr_fes[i]);
-      HypreParMatrix * Pi = (HypreParMatrix*)(&pP->GetBlock(skip_i,skip_i));
+      HypreParMatrix * Pi = &pP->GetBlock(skip_i,skip_i);
       int skip_j=0;
       for (int j = 0; j<nblocks; j++)
       {
@@ -663,7 +663,7 @@ void ComplexBlockStaticCondensation::ParallelAssemble(BlockMatrix *m_r,
          else
          {
             pfes_j = dynamic_cast<ParFiniteElementSpace*>(tr_fes[j]);
-            HypreParMatrix * Pj = (HypreParMatrix*)(&pP->GetBlock(skip_j,skip_j));
+            HypreParMatrix * Pj = &pP->GetBlock(skip_j,skip_j);
             A_r = new HypreParMatrix(pfes_i->GetComm(), pfes_i->GlobalVSize(),
                                      pfes_j->GlobalVSize(), pfes_i->GetDofOffsets(),
                                      pfes_j->GetDofOffsets(), &m_r->GetBlock(skip_i,skip_j));

miniapps/dpg/util/complexstaticcond.hpp

@@ -65,12 +65,12 @@ class ComplexBlockStaticCondensation
    BlockMatrix * R = nullptr; // Block Restriction
 
 #ifdef MFEM_USE_MPI
-   BlockOperator * pS_r = nullptr;
-   BlockOperator * pS_e_r = nullptr;
-   BlockOperator * pS_i = nullptr;
-   BlockOperator * pS_e_i = nullptr;
+   TBlockOperator<HypreParMatrix> * pS_r = nullptr;
+   TBlockOperator<HypreParMatrix> * pS_e_r = nullptr;
+   TBlockOperator<HypreParMatrix> * pS_i = nullptr;
+   TBlockOperator<HypreParMatrix> * pS_e_i = nullptr;
    // Block HypreParMatrix for Prolongation
-   BlockOperator * pP = nullptr;
+   TBlockOperator<HypreParMatrix> * pP = nullptr;
 #endif
 
    bool Parallel() const { return parallel; }
@@ -179,12 +179,12 @@ public:
 
 #ifdef MFEM_USE_MPI
    /// Return the parallel Schur complement matrix.
-   BlockOperator &GetParallelSchurMatrix_r() { return *pS_r; }
-   BlockOperator &GetParallelSchurMatrix_i() { return *pS_i; }
+   TBlockOperator<HypreParMatrix> &GetParallelSchurMatrix_r() { return *pS_r; }
+   TBlockOperator<HypreParMatrix> &GetParallelSchurMatrix_i() { return *pS_i; }
 
    /// Return the eliminated part of the parallel Schur complement matrix.
-   BlockOperator &GetParallelSchurMatrixElim_r() { return *pS_e_r; }
-   BlockOperator &GetParallelSchurMatrixElim_i() { return *pS_e_i; }
+   TBlockOperator<HypreParMatrix> &GetParallelSchurMatrixElim_r() { return *pS_e_r; }
+   TBlockOperator<HypreParMatrix> &GetParallelSchurMatrixElim_i() { return *pS_e_i; }
 
    void ParallelAssemble(BlockMatrix *m_r, BlockMatrix*m_i);
 #endif

miniapps/dpg/util/pcomplexweakform.cpp

@@ -41,10 +41,10 @@ void ParComplexDPGWeakForm::ParallelAssemble(BlockMatrix *m_r,
 {
    if (!P) { BuildProlongation(); }
 
-   p_mat_r = new BlockOperator(tdof_offsets);
-   p_mat_i = new BlockOperator(tdof_offsets);
-   p_mat_e_r = new BlockOperator(tdof_offsets);
-   p_mat_e_i = new BlockOperator(tdof_offsets);
+   p_mat_r = new TBlockOperator<HypreParMatrix>(tdof_offsets);
+   p_mat_i = new TBlockOperator<HypreParMatrix>(tdof_offsets);
+   p_mat_e_r = new TBlockOperator<HypreParMatrix>(tdof_offsets);
+   p_mat_e_i = new TBlockOperator<HypreParMatrix>(tdof_offsets);
    p_mat_r->owns_blocks = 1;
    p_mat_i->owns_blocks = 1;
    p_mat_e_r->owns_blocks = 1;
@@ -55,7 +55,7 @@ void ParComplexDPGWeakForm::ParallelAssemble(BlockMatrix *m_r,
    HypreParMatrix * PtAP_i = nullptr;
    for (int i = 0; i < nblocks; i++)
    {
-      HypreParMatrix * Pi = (HypreParMatrix*)(&P->GetBlock(i,i));
+      HypreParMatrix * Pi = &P->GetBlock(i,i);
       for (int j = 0; j<nblocks; j++)
       {
          if (m_r->IsZeroBlock(i,j)) { continue; }
@@ -78,7 +78,7 @@ void ParComplexDPGWeakForm::ParallelAssemble(BlockMatrix *m_r,
          }
          else
          {
-            HypreParMatrix * Pj = (HypreParMatrix*)(&P->GetBlock(j,j));
+            HypreParMatrix * Pj = &P->GetBlock(j,j);
             A_r = new HypreParMatrix(trial_pfes[i]->GetComm(), trial_pfes[i]->GlobalVSize(),
                                      trial_pfes[j]->GlobalVSize(), trial_pfes[i]->GetDofOffsets(),
                                      trial_pfes[j]->GetDofOffsets(), &m_r->GetBlock(i,j));
@@ -104,7 +104,7 @@ void ParComplexDPGWeakForm::ParallelAssemble(BlockMatrix *m_r,
 
 void ParComplexDPGWeakForm::BuildProlongation()
 {
-   P = new BlockOperator(dof_offsets, tdof_offsets);
+   P = new TBlockOperator<HypreParMatrix>(dof_offsets, tdof_offsets);
    R = new BlockMatrix(tdof_offsets, dof_offsets);
    P->owns_blocks = 0;
    R->owns_blocks = 0;

miniapps/dpg/util/pcomplexweakform.hpp

@@ -37,15 +37,15 @@ protected:
    void FillEssTdofLists(const Array<int> & ess_tdof_list);
 
    // Block Prolongation
-   BlockOperator * P = nullptr;
+   TBlockOperator<HypreParMatrix> * P = nullptr;
    // Block Restriction
    BlockMatrix * R = nullptr;
 
    ComplexOperator * p_mat = nullptr;
-   BlockOperator * p_mat_r = nullptr;
-   BlockOperator * p_mat_i = nullptr;
-   BlockOperator * p_mat_e_r = nullptr;
-   BlockOperator * p_mat_e_i = nullptr;
+   TBlockOperator<HypreParMatrix> * p_mat_r = nullptr;
+   TBlockOperator<HypreParMatrix> * p_mat_i = nullptr;
+   TBlockOperator<HypreParMatrix> * p_mat_e_r = nullptr;
+   TBlockOperator<HypreParMatrix> * p_mat_e_i = nullptr;
 
    void BuildProlongation();
 

miniapps/dpg/util/pweakform.cpp

@@ -39,15 +39,15 @@ void ParDPGWeakForm::ParallelAssemble(BlockMatrix *m)
 {
    if (!P) { BuildProlongation(); }
 
-   p_mat = new BlockOperator(tdof_offsets);
-   p_mat_e = new BlockOperator(tdof_offsets);
+   p_mat = new TBlockOperator<HypreParMatrix>(tdof_offsets);
+   p_mat_e = new TBlockOperator<HypreParMatrix>(tdof_offsets);
    p_mat->owns_blocks = 1;
    p_mat_e->owns_blocks = 1;
    HypreParMatrix * A = nullptr;
    HypreParMatrix * PtAP = nullptr;
    for (int i = 0; i<nblocks; i++)
    {
-      HypreParMatrix * Pi = (HypreParMatrix*)(&P->GetBlock(i,i));
+      HypreParMatrix * Pi = &P->GetBlock(i,i);
       HypreParMatrix * Pit = Pi->Transpose();
       for (int j = 0; j<nblocks; j++)
       {
@@ -63,7 +63,7 @@ void ParDPGWeakForm::ParallelAssemble(BlockMatrix *m)
          }
          else
          {
-            HypreParMatrix * Pj = (HypreParMatrix*)(&P->GetBlock(j,j));
+            HypreParMatrix * Pj = &P->GetBlock(j,j);
             A = new HypreParMatrix(trial_pfes[i]->GetComm(), trial_pfes[i]->GlobalVSize(),
                                    trial_pfes[j]->GlobalVSize(), trial_pfes[i]->GetDofOffsets(),
                                    trial_pfes[j]->GetDofOffsets(), &m->GetBlock(i,j));
@@ -83,7 +83,7 @@ void ParDPGWeakForm::ParallelAssemble(BlockMatrix *m)
 
 void ParDPGWeakForm::BuildProlongation()
 {
-   P = new BlockOperator(dof_offsets, tdof_offsets);
+   P = new TBlockOperator<HypreParMatrix>(dof_offsets, tdof_offsets);
    R = new BlockMatrix(tdof_offsets, dof_offsets);
    P->owns_blocks = 0;
    R->owns_blocks = 0;

miniapps/dpg/util/pweakform.hpp

@@ -37,12 +37,12 @@ protected:
    void FillEssTdofLists(const Array<int> & ess_tdof_list);
 
    // Block operator of HypreParMatrix
-   BlockOperator * P = nullptr; // Block Prolongation
+   TBlockOperator<HypreParMatrix> * P = nullptr; // Block Prolongation
    BlockMatrix * R = nullptr; // Block Restriction
 
    // Block operator of HypreParMatrix
-   BlockOperator * p_mat = nullptr;
-   BlockOperator * p_mat_e = nullptr;
+   TBlockOperator<HypreParMatrix> * p_mat = nullptr;
+   TBlockOperator<HypreParMatrix> * p_mat_e = nullptr;
 
    void BuildProlongation();
 

miniapps/solvers/darcy_solver.hpp

@@ -44,7 +44,7 @@ public:
 /// Wrapper for the block-diagonal-preconditioned MINRES employed in ex5p.cpp
 class BDPMinresSolver : public DarcySolver
 {
-   BlockOperator op_;
+   BlockConstOperator op_;
    BlockDiagonalPreconditioner prec_;
    OperatorPtr BT_;
    OperatorPtr S_;   // S_ = B diag(M)^{-1} B^T

miniapps/solvers/div_free_solver.cpp

@@ -319,15 +319,15 @@ DivFreeSolver::DivFreeSolver(const HypreParMatrix &M, const HypreParMatrix& B,
      blk_Ps_(ops_.Size()-1), smoothers_(ops_.Size())
 {
    ops_offsets_.back().MakeRef(DarcySolver::offsets_);
-   ops_.Last() = new BlockOperator(ops_offsets_.back());
-   ops_.Last()->SetBlock(0, 0, const_cast<HypreParMatrix*>(&M));
-   ops_.Last()->SetBlock(1, 0, const_cast<HypreParMatrix*>(&B));
-   ops_.Last()->SetBlock(0, 1, BT_.Ptr());
+   ops_.Last() = new TBlockOperator<const HypreParMatrix>(ops_offsets_.back());
+   ops_.Last()->SetBlock(0, 0, &M);
+   ops_.Last()->SetBlock(1, 0, &B);
+   ops_.Last()->SetBlock(0, 1, BT_.As<HypreParMatrix>());
 
    for (int l = data.P_l2.size(); l >= 0; --l)
    {
-      auto& M_f = static_cast<const HypreParMatrix&>(ops_[l]->GetBlock(0, 0));
-      auto& B_f = static_cast<const HypreParMatrix&>(ops_[l]->GetBlock(1, 0));
+      auto& M_f = ops_[l]->GetBlock(0, 0);
+      auto& B_f = ops_[l]->GetBlock(1, 0);
 
       if (l == 0)
       {
@@ -378,11 +378,12 @@ DivFreeSolver::DivFreeSolver(const HypreParMatrix &M, const HypreParMatrix& B,
       ops_offsets_[l-1][1] = M_c->NumRows();
       ops_offsets_[l-1][2] = M_c->NumRows() + B_c->NumRows();
 
-      blk_Ps_[l-1] = new BlockOperator(ops_offsets_[l], ops_offsets_[l-1]);
+      blk_Ps_[l-1] = new TBlockOperator<const HypreParMatrix>(ops_offsets_[l],
+                                                              ops_offsets_[l-1]);
       blk_Ps_[l-1]->SetBlock(0, 0, &P_hdiv_l);
       blk_Ps_[l-1]->SetBlock(1, 1, &P_l2_l);
 
-      ops_[l-1] = new BlockOperator(ops_offsets_[l-1]);
+      ops_[l-1] = new TBlockOperator<const HypreParMatrix>(ops_offsets_[l-1]);
       ops_[l-1]->SetBlock(0, 0, M_c);
       ops_[l-1]->SetBlock(1, 0, B_c);
       ops_[l-1]->SetBlock(0, 1, B_c->Transpose());

miniapps/solvers/div_free_solver.hpp

@@ -178,8 +178,8 @@ class DivFreeSolver : public DarcySolver
    OperatorPtr BT_;
    BBTSolver BBT_solver_;
    std::vector<Array<int>> ops_offsets_;
-   Array<BlockOperator*> ops_;
-   Array<BlockOperator*> blk_Ps_;
+   Array<TBlockOperator<const HypreParMatrix>*> ops_;
+   Array<TBlockOperator<const HypreParMatrix>*> blk_Ps_;
    Array<Solver*> smoothers_;
    OperatorPtr prec_;
    OperatorPtr solver_;