Fix single precision build.

examples/ex10.cpp

@@ -412,7 +412,7 @@ void ReducedSystemOperator::Mult(const Vector &k, Vector &y) const
 Operator &ReducedSystemOperator::GetGradient(const Vector &k) const
 {
    delete Jacobian;
-   Jacobian = Add(1.0, M->SpMat(), dt, S->SpMat());
+   Jacobian = Add((real_t)1.0, M->SpMat(), dt, S->SpMat());
    add(*v, dt, k, w);
    add(*x, dt, w, z);
    SparseMatrix *grad_H = dynamic_cast<SparseMatrix *>(&H->GetGradient(z));

examples/ex10p.cpp

@@ -476,7 +476,7 @@ void ReducedSystemOperator::Mult(const Vector &k, Vector &y) const
 Operator &ReducedSystemOperator::GetGradient(const Vector &k) const
 {
    delete Jacobian;
-   SparseMatrix *localJ = Add(1.0, M->SpMat(), dt, S->SpMat());
+   SparseMatrix *localJ = Add((real_t)1.0, M->SpMat(), dt, S->SpMat());
    add(*v, dt, k, w);
    add(*x, dt, w, z);
    localJ->Add(dt*dt, H->GetLocalGradient(z));

examples/ex16.cpp

@@ -323,7 +323,7 @@ void ConductionOperator::ImplicitSolve(const real_t dt,
    // for du_dt, where K is linearized by using u from the previous timestep
    if (!T)
    {
-      T = Add(1.0, Mmat, dt, Kmat);
+      T = Add((real_t)1.0, Mmat, dt, Kmat);
       current_dt = dt;
       T_solver.SetOperator(*T);
    }

examples/ex16p.cpp

@@ -414,7 +414,7 @@ void ConductionOperator::ImplicitSolve(const real_t dt,
    // for du_dt, where K is linearized by using u from the previous timestep
    if (!T)
    {
-      T = Add(1.0, Mmat, dt, Kmat);
+      T = Add((real_t)1.0, Mmat, dt, Kmat);
       current_dt = dt;
       T_solver.SetOperator(*T);
    }

examples/ex23.cpp

@@ -139,7 +139,7 @@ void WaveOperator::ImplicitSolve(const real_t fac0, const real_t fac1,
    // for d2udt2
    if (!T)
    {
-      T = Add(1.0, Mmat, fac0, Kmat);
+      T = Add((real_t)1.0, Mmat, fac0, Kmat);
       T_solver.SetOperator(*T);
    }
    K->FullMult(u, z);

examples/ex3.cpp

@@ -56,6 +56,51 @@ void f_exact(const Vector &, Vector &);
 real_t freq = 1.0, kappa;
 int dim;
 
+void SolveSingle(SparseMatrix &A, const Vector &B, Vector &X)
+{
+   VectorMP<float> Bs, Xs;
+   const real_t *data = A.GetData();
+
+   const int n = A.GetI()[A.NumRows()];
+   int *Icopy = new int[A.NumRows() + 1];
+   int *Jcopy = new int[n];
+
+   float *sdata = new float[n];
+   for (int i=0; i<n; ++i)
+   {
+      sdata[i] = data[i];
+      Jcopy[i] = A.GetJ()[i];
+   }
+
+   for (int i=0; i<A.NumRows() + 1; ++i)
+   {
+      Icopy[i] = A.GetI()[i];
+   }
+
+   SparseMatrixMP<float> As(Icopy, Jcopy, sdata, A.NumRows(), A.NumCols());
+
+   Bs.SetSize(B.Size());
+   Xs.SetSize(X.Size());
+
+   for (int i=0; i<B.Size(); ++i)
+   {
+      Bs[i] = B[i];
+   }
+
+   for (int i=0; i<X.Size(); ++i)
+   {
+      Xs[i] = X[i];
+   }
+
+   GSSmootherMP<float> Ms(As);
+   PCG<float>(As, Ms, Bs, Xs, 1, 500, 1e-12, 0.0);
+
+   for (int i=0; i<X.Size(); ++i)
+   {
+      X[i] = Xs[i];
+   }
+}
+
 int main(int argc, char *argv[])
 {
    // 1. Parse command-line options.
@@ -185,6 +230,7 @@ int main(int argc, char *argv[])
 
    cout << "Size of linear system: " << A->Height() << endl;
 
+   /*
    // 11. Solve the linear system A X = B.
    if (pa) // Jacobi preconditioning in partial assembly mode
    {
@@ -193,20 +239,23 @@ int main(int argc, char *argv[])
    }
    else
    {
-#ifndef MFEM_USE_SUITESPARSE
+   #ifndef MFEM_USE_SUITESPARSE
       // 11. Define a simple symmetric Gauss-Seidel preconditioner and use it to
       //     solve the system Ax=b with PCG.
       GSSmoother M((SparseMatrix&)(*A));
       PCG(*A, M, B, X, 1, 500, 1e-12, 0.0);
-#else
+   #else
       // 11. If MFEM was compiled with SuiteSparse, use UMFPACK to solve the
       //     system.
       UMFPackSolver umf_solver;
       umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
       umf_solver.SetOperator(*A);
       umf_solver.Mult(B, X);
-#endif
+   #endif
    }
+   */
+
+   SolveSingle((SparseMatrix&)(*A), B, X);
 
    // 12. Recover the solution as a finite element grid function.
    a->RecoverFEMSolution(X, *b, x);

examples/ex3p.cpp

@@ -297,6 +297,18 @@ int main(int argc, char *argv[])
       sol_sock << "solution\n" << *pmesh << x << flush;
    }
 
+   VectorMP<float> xf(x.Size());
+   for (int i=0; i<x.Size(); ++i)
+   {
+      xf[i] = x[i];
+   }
+
+   if (myid == 0)
+   {
+      cout << "Norm of x " << x.Norml2() << endl;
+      cout << "Norm of xf " << xf.Norml2() << endl;
+   }
+
    // 18. Free the used memory.
    delete a;
    delete sigma;

fem/ceed/interface/util.hpp

@@ -33,7 +33,6 @@ class FiniteElement;
 class FiniteElementSpace;
 class ElementTransformation;
 class IntegrationRule;
-class Vector;
 
 /** @brief Function that determines if a CEED kernel should be used, based on
     the current mfem::Device configuration. */

fem/fe/fe_nurbs.cpp

@@ -81,7 +81,7 @@ void NURBS1DFiniteElement::CalcHessian (const IntegrationPoint &ip,
 
    sum = 1.0/sum;
    add(sum, hess, -2*dsum*sum*sum, grad, hess);
-   add(1.0, hess, (-d2sum + 2*dsum*dsum*sum)*sum*sum, shape_x, hess);
+   add((real_t)1.0, hess, (-d2sum + 2*dsum*dsum*sum)*sum*sum, shape_x, hess);
 }
 
 

fem/fe/fe_pyramid.cpp

@@ -1574,7 +1574,7 @@ void FuentesPyramid::V_R(int p, Vector s, const DenseMatrix &grad_s,
    {
       // dphi_E_i.GetRow(i, dphi);
       for (int l=0; l<3; l++) { dphi[l] = dphi_E_i(i, l); }
-      add(t * t, dphi, 2.0 * t * phi_E_i(i), dt3, dphit2);
+      add(t * t, dphi, 2 * t * phi_E_i(i), dt3, dphit2);
       dphit2.cross3D(dmu3, dphixdmu);
       // u.SetRow(i, dphixdmu);
       for (int l=0; l<3; l++) { u(i, l) = dphixdmu(l); }

fem/linearform_ext.hpp

@@ -18,7 +18,6 @@
 namespace mfem
 {
 
-class Operator;
 class LinearForm;
 
 /// Class extending the LinearForm class to support assembly on devices.

general/optparser.hpp

@@ -14,12 +14,11 @@
 
 #include "../config/config.hpp"
 #include "array.hpp"
+#include "../linalg/vector.hpp"
 
 namespace mfem
 {
 
-class Vector;
-
 /** Class for parsing command-line options.
 
     The class is initialized with argc and argv, and new options are added with

linalg/matrix.cpp

@@ -19,7 +19,8 @@
 namespace mfem
 {
 
-void Matrix::Print (std::ostream & os, int width_) const
+template <class T>
+void MatrixMP<T>::Print(std::ostream & os, int width_) const
 {
    using namespace std;
    // output flags = scientific + show sign
@@ -40,4 +41,10 @@ void Matrix::Print (std::ostream & os, int width_) const
    os << '\n';
 }
 
+template class MatrixMP<float>;
+template class MatrixMP<double>;
+
+template class AbstractSparseMatrixMP<float>;
+template class AbstractSparseMatrixMP<double>;
+
 }

linalg/matrix.hpp

@@ -21,31 +21,39 @@ namespace mfem
 
 // Abstract data types matrix, inverse matrix
 
-class MatrixInverse;
+template <class T>
+class MatrixInverseMP;
 
 /// Abstract data type matrix
-class Matrix : public Operator
+
+template <class T>
+class MatrixMP : public OperatorMP<T>
 {
-   friend class MatrixInverse;
+   friend class MatrixInverseMP<T>;
+
+protected:
+   using OperatorBase::height;
+   using OperatorBase::width;
+
 public:
 
    /// Creates a square matrix of size s.
-   explicit Matrix(int s) : Operator(s) { }
+   explicit MatrixMP(int s) : OperatorMP<T>(s) { }
 
    /// Creates a matrix of the given height and width.
-   explicit Matrix(int h, int w) : Operator(h, w) { }
+   explicit MatrixMP(int h, int w) : OperatorMP<T>(h, w) { }
 
    /// Returns whether the matrix is a square matrix.
    bool IsSquare() const { return (height == width); }
 
    /// Returns reference to a_{ij}.
-   virtual real_t &Elem(int i, int j) = 0;
+   virtual T &Elem(int i, int j) = 0;
 
    /// Returns constant reference to a_{ij}.
-   virtual const real_t &Elem(int i, int j) const = 0;
+   virtual const T &Elem(int i, int j) const = 0;
 
    /// Returns a pointer to (an approximation) of the matrix inverse.
-   virtual MatrixInverse *Inverse() const = 0;
+   virtual MatrixInverseMP<T> *Inverse() const = 0;
 
    /// Finalizes the matrix initialization.
    virtual void Finalize(int) { }
@@ -54,30 +62,35 @@ public:
    virtual void Print(std::ostream & os = mfem::out, int width_ = 4) const;
 
    /// Destroys matrix.
-   virtual ~Matrix() { }
+   virtual ~MatrixMP() { }
 };
 
+using Matrix = MatrixMP<real_t>;
 
 /// Abstract data type for matrix inverse
-class MatrixInverse : public Solver
+template <class T>
+class MatrixInverseMP : public SolverMP<T>
 {
 public:
-   MatrixInverse() { }
+   MatrixInverseMP() { }
 
    /// Creates approximation of the inverse of square matrix
-   MatrixInverse(const Matrix &mat)
-      : Solver(mat.height, mat.width) { }
+   MatrixInverseMP(const MatrixMP<T> &mat)
+      : SolverMP<T>(mat.height, mat.width) { }
 };
 
+using MatrixInverse = MatrixInverseMP<real_t>;
+
 /// Abstract data type for sparse matrices
-class AbstractSparseMatrix : public Matrix
+template <class T>
+class AbstractSparseMatrixMP : public MatrixMP<T>
 {
 public:
    /// Creates a square matrix of the given size.
-   explicit AbstractSparseMatrix(int s = 0) : Matrix(s) { }
+   explicit AbstractSparseMatrixMP(int s = 0) : MatrixMP<T>(s) { }
 
    /// Creates a matrix of the given height and width.
-   explicit AbstractSparseMatrix(int h, int w) : Matrix(h, w) { }
+   explicit AbstractSparseMatrixMP(int h, int w) : MatrixMP<T>(h, w) { }
 
    /// Returns the number of non-zeros in a matrix
    virtual int NumNonZeroElems() const = 0;
@@ -86,30 +99,33 @@ public:
    /** Returns:
        - 0 if @a cols and @a srow are copies of the values in the matrix.
        - 1 if @a cols and @a srow are views of the values in the matrix. */
-   virtual int GetRow(const int row, Array<int> &cols, Vector &srow) const = 0;
+   virtual int GetRow(const int row, Array<int> &cols,
+                      VectorMP<T> &srow) const = 0;
 
    /** @brief If the matrix is square, this method will place 1 on the diagonal
        (i,i) if row i has "almost" zero l1-norm.
 
        If entry (i,i) does not belong to the sparsity pattern of A, then an
        error will occur. */
-   virtual void EliminateZeroRows(const real_t threshold = 1e-12) = 0;
+   virtual void EliminateZeroRows(const T threshold = 1e-12) = 0;
 
    /// Matrix-Vector Multiplication y = A*x
-   void Mult(const Vector &x, Vector &y) const override = 0;
+   void Mult(const VectorMP<T> &x, VectorMP<T> &y) const override = 0;
    /// Matrix-Vector Multiplication y = y + val*A*x
-   void AddMult(const Vector &x, Vector &y,
-                const real_t val = 1.) const override = 0;
+   void AddMult(const VectorMP<T> &x, VectorMP<T> &y,
+                const T val = 1.) const override = 0;
    /// MatrixTranspose-Vector Multiplication y = A'*x
-   void MultTranspose(const Vector &x, Vector &y) const override = 0;
+   void MultTranspose(const VectorMP<T> &x, VectorMP<T> &y) const override = 0;
    /// MatrixTranspose-Vector Multiplication y = y + val*A'*x
-   void AddMultTranspose(const Vector &x, Vector &y,
-                         const real_t val = 1.) const override = 0;
+   void AddMultTranspose(const VectorMP<T> &x, VectorMP<T> &y,
+                         const T val = 1.) const override = 0;
 
    /// Destroys AbstractSparseMatrix.
-   virtual ~AbstractSparseMatrix() { }
+   virtual ~AbstractSparseMatrixMP() { }
 };
 
+using AbstractSparseMatrix = AbstractSparseMatrixMP<real_t>;
+
 }
 
 #endif

linalg/mma.hpp

@@ -23,7 +23,11 @@
 namespace mfem
 {
 // forward declaration
-class Vector;
+template <class T>
+class VectorMP;
+
+using Vector = VectorMP<real_t>;
+
 
 /** \brief MMA (Method of Moving Asymptotes) solves a nonlinear optimization
  *         problem involving an objective function, inequality constraints,

linalg/operator.cpp

@@ -19,10 +19,12 @@
 namespace mfem
 {
 
-void Operator::InitTVectors(const Operator *Po, const Operator *Ri,
-                            const Operator *Pi,
-                            Vector &x, Vector &b,
-                            Vector &X, Vector &B) const
+template <class T>
+void OperatorMP<T>::InitTVectors(const OperatorMP<T> *Po,
+                                 const OperatorMP<T> *Ri,
+                                 const OperatorMP<T> *Pi,
+                                 VectorMP<T> &x, VectorMP<T> &b,
+                                 VectorMP<T> &X, VectorMP<T> &B) const
 {
    if (!IsIdentityProlongation(Po))
    {
@@ -48,23 +50,27 @@ void Operator::InitTVectors(const Operator *Po, const Operator *Ri,
    }
 }
 
-void Operator::AddMult(const Vector &x, Vector &y, const real_t a) const
+template <class T>
+void OperatorMP<T>::AddMult(const VectorMP<T> &x, VectorMP<T> &y,
+                            const T a) const
 {
-   mfem::Vector z(y.Size());
+   mfem::VectorMP<T> z(y.Size());
    Mult(x, z);
    y.Add(a, z);
 }
 
-void Operator::AddMultTranspose(const Vector &x, Vector &y,
-                                const real_t a) const
+template <class T>
+void OperatorMP<T>::AddMultTranspose(const VectorMP<T> &x, VectorMP<T> &y,
+                                     const T a) const
 {
-   mfem::Vector z(y.Size());
+   mfem::VectorMP<T> z(y.Size());
    MultTranspose(x, z);
    y.Add(a, z);
 }
 
-void Operator::ArrayMult(const Array<const Vector *> &X,
-                         Array<Vector *> &Y) const
+template <class T>
+void OperatorMP<T>::ArrayMult(const Array<const VectorMP<T> *> &X,
+                              Array<VectorMP<T> *> &Y) const
 {
    MFEM_ASSERT(X.Size() == Y.Size(),
                "Number of columns mismatch in Operator::Mult!");
@@ -75,8 +81,9 @@ void Operator::ArrayMult(const Array<const Vector *> &X,
    }
 }
 
-void Operator::ArrayMultTranspose(const Array<const Vector *> &X,
-                                  Array<Vector *> &Y) const
+template <class T>
+void OperatorMP<T>::ArrayMultTranspose(const Array<const VectorMP<T> *> &X,
+                                       Array<VectorMP<T> *> &Y) const
 {
    MFEM_ASSERT(X.Size() == Y.Size(),
                "Number of columns mismatch in Operator::MultTranspose!");
@@ -87,8 +94,10 @@ void Operator::ArrayMultTranspose(const Array<const Vector *> &X,
    }
 }
 
-void Operator::ArrayAddMult(const Array<const Vector *> &X, Array<Vector *> &Y,
-                            const real_t a) const
+template <class T>
+void OperatorMP<T>::ArrayAddMult(const Array<const VectorMP<T> *> &X,
+                                 Array<VectorMP<T> *> &Y,
+                                 const T a) const
 {
    MFEM_ASSERT(X.Size() == Y.Size(),
                "Number of columns mismatch in Operator::AddMult!");
@@ -99,8 +108,9 @@ void Operator::ArrayAddMult(const Array<const Vector *> &X, Array<Vector *> &Y,
    }
 }
 
-void Operator::ArrayAddMultTranspose(const Array<const Vector *> &X,
-                                     Array<Vector *> &Y, const real_t a) const
+template <class T>
+void OperatorMP<T>::ArrayAddMultTranspose(const Array<const VectorMP<T> *> &X,
+                                          Array<VectorMP<T> *> &Y, const T a) const
 {
    MFEM_ASSERT(X.Size() == Y.Size(),
                "Number of columns mismatch in Operator::AddMultTranspose!");
@@ -111,44 +121,48 @@ void Operator::ArrayAddMultTranspose(const Array<const Vector *> &X,
    }
 }
 
-void Operator::FormLinearSystem(const Array<int> &ess_tdof_list,
-                                Vector &x, Vector &b,
-                                Operator* &Aout, Vector &X, Vector &B,
-                                int copy_interior)
+template <class T>
+void OperatorMP<T>::FormLinearSystem(const Array<int> &ess_tdof_list,
+                                     VectorMP<T> &x, VectorMP<T> &b,
+                                     OperatorMP<T>* &Aout, VectorMP<T> &X, VectorMP<T> &B,
+                                     int copy_interior)
 {
-   const Operator *P = this->GetProlongation();
-   const Operator *R = this->GetRestriction();
+   const OperatorMP<T> *P = this->GetProlongation();
+   const OperatorMP<T> *R = this->GetRestriction();
    InitTVectors(P, R, P, x, b, X, B);
 
    if (!copy_interior) { X.SetSubVectorComplement(ess_tdof_list, 0.0); }
 
-   ConstrainedOperator *constrainedA;
+   ConstrainedOperatorMP<T> *constrainedA;
    FormConstrainedSystemOperator(ess_tdof_list, constrainedA);
    constrainedA->EliminateRHS(X, B);
    Aout = constrainedA;
 }
 
-void Operator::FormRectangularLinearSystem(
+template <class T>
+void OperatorMP<T>::FormRectangularLinearSystem(
    const Array<int> &trial_tdof_list,
-   const Array<int> &test_tdof_list, Vector &x, Vector &b,
-   Operator* &Aout, Vector &X, Vector &B)
+   const Array<int> &test_tdof_list, VectorMP<T> &x, VectorMP<T> &b,
+   OperatorMP<T>* &Aout, VectorMP<T> &X, VectorMP<T> &B)
 {
-   const Operator *Pi = this->GetProlongation();
-   const Operator *Po = this->GetOutputProlongation();
-   const Operator *Ri = this->GetRestriction();
+   const OperatorMP<T> *Pi = this->GetProlongation();
+   const OperatorMP<T> *Po = this->GetOutputProlongation();
+   const OperatorMP<T> *Ri = this->GetRestriction();
    InitTVectors(Po, Ri, Pi, x, b, X, B);
 
-   RectangularConstrainedOperator *constrainedA;
+   RectangularConstrainedOperatorMP<T> *constrainedA;
    FormRectangularConstrainedSystemOperator(trial_tdof_list, test_tdof_list,
                                             constrainedA);
    constrainedA->EliminateRHS(X, B);
    Aout = constrainedA;
 }
 
-void Operator::RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
+template <class T>
+void OperatorMP<T>::RecoverFEMSolution(const VectorMP<T> &X,
+                                       const VectorMP<T> &b, VectorMP<T> &x)
 {
    // Same for Rectangular and Square operators
-   const Operator *P = this->GetProlongation();
+   const OperatorMP<T> *P = this->GetProlongation();
    if (!IsIdentityProlongation(P))
    {
       // Apply conforming prolongation
@@ -165,26 +179,28 @@ void Operator::RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
    }
 }
 
-Operator * Operator::SetupRAP(const Operator *Pi, const Operator *Po)
+template <class T>
+OperatorMP<T> * OperatorMP<T>::SetupRAP(const OperatorMP<T> *Pi,
+                                        const OperatorMP<T> *Po)
 {
-   Operator *rap;
+   OperatorMP<T> *rap;
    if (!IsIdentityProlongation(Pi))
    {
       if (!IsIdentityProlongation(Po))
       {
-         rap = new RAPOperator(*Po, *this, *Pi);
+         rap = new RAPOperatorMP<T>(*Po, *this, *Pi);
       }
       else
       {
-         rap = new ProductOperator(this, Pi, false,false);
+         rap = new ProductOperatorMP<T>(this, Pi, false, false);
       }
    }
    else
    {
       if (!IsIdentityProlongation(Po))
       {
-         TransposeOperator * PoT = new TransposeOperator(Po);
-         rap = new ProductOperator(PoT, this, true,false);
+         TransposeOperatorMP<T> * PoT = new TransposeOperatorMP<T>(Po);
+         rap = new ProductOperatorMP<T>(PoT, this, true, false);
       }
       else
       {
@@ -194,67 +210,74 @@ Operator * Operator::SetupRAP(const Operator *Pi, const Operator *Po)
    return rap;
 }
 
-void Operator::FormConstrainedSystemOperator(
-   const Array<int> &ess_tdof_list, ConstrainedOperator* &Aout)
+template <class T>
+void OperatorMP<T>::FormConstrainedSystemOperator(
+   const Array<int> &ess_tdof_list, ConstrainedOperatorMP<T>* &Aout)
 {
-   const Operator *P = this->GetProlongation();
-   Operator *rap = SetupRAP(P, P);
+   const OperatorMP<T> *P = this->GetProlongation();
+   OperatorMP<T> *rap = SetupRAP(P, P);
 
    // Impose the boundary conditions through a ConstrainedOperator, which owns
    // the rap operator when P and R are non-trivial
-   ConstrainedOperator *A = new ConstrainedOperator(rap, ess_tdof_list,
-                                                    rap != this);
+   ConstrainedOperatorMP<T> *A = new ConstrainedOperatorMP<T>(rap, ess_tdof_list,
+                                                              rap != this);
    Aout = A;
 }
 
-void Operator::FormRectangularConstrainedSystemOperator(
+template <class T>
+void OperatorMP<T>::FormRectangularConstrainedSystemOperator(
    const Array<int> &trial_tdof_list, const Array<int> &test_tdof_list,
-   RectangularConstrainedOperator* &Aout)
+   RectangularConstrainedOperatorMP<T>* &Aout)
 {
-   const Operator *Pi = this->GetProlongation();
-   const Operator *Po = this->GetOutputProlongation();
-   Operator *rap = SetupRAP(Pi, Po);
+   const OperatorMP<T> *Pi = this->GetProlongation();
+   const OperatorMP<T> *Po = this->GetOutputProlongation();
+   OperatorMP<T> *rap = SetupRAP(Pi, Po);
 
    // Impose the boundary conditions through a RectangularConstrainedOperator,
    // which owns the rap operator when P and R are non-trivial
-   RectangularConstrainedOperator *A
-      = new RectangularConstrainedOperator(rap,
-                                           trial_tdof_list, test_tdof_list,
-                                           rap != this);
+   RectangularConstrainedOperatorMP<T> *A
+      = new RectangularConstrainedOperatorMP<T>(rap,
+                                                trial_tdof_list, test_tdof_list,
+                                                rap != this);
    Aout = A;
 }
 
-void Operator::FormSystemOperator(const Array<int> &ess_tdof_list,
-                                  Operator* &Aout)
+template <class T>
+void OperatorMP<T>::FormSystemOperator(const Array<int> &ess_tdof_list,
+                                       OperatorMP<T>* &Aout)
 {
-   ConstrainedOperator *A;
+   ConstrainedOperatorMP<T> *A;
    FormConstrainedSystemOperator(ess_tdof_list, A);
    Aout = A;
 }
 
-void Operator::FormRectangularSystemOperator(const Array<int> &trial_tdof_list,
-                                             const Array<int> &test_tdof_list,
-                                             Operator* &Aout)
+template <class T>
+void OperatorMP<T>::FormRectangularSystemOperator(const Array<int>
+                                                  &trial_tdof_list,
+                                                  const Array<int> &test_tdof_list,
+                                                  OperatorMP<T>* &Aout)
 {
-   RectangularConstrainedOperator *A;
+   RectangularConstrainedOperatorMP<T> *A;
    FormRectangularConstrainedSystemOperator(trial_tdof_list, test_tdof_list, A);
    Aout = A;
 }
 
-void Operator::FormDiscreteOperator(Operator* &Aout)
+template <class T>
+void OperatorMP<T>::FormDiscreteOperator(OperatorMP<T>* &Aout)
 {
-   const Operator *Pin  = this->GetProlongation();
-   const Operator *Rout = this->GetOutputRestriction();
-   Aout = new TripleProductOperator(Rout, this, Pin,false, false, false);
+   const OperatorMP<T> *Pin  = this->GetProlongation();
+   const OperatorMP<T> *Rout = this->GetOutputRestriction();
+   Aout = new TripleProductOperatorMP<T>(Rout, this, Pin, false, false, false);
 }
 
-void Operator::PrintMatlab(std::ostream & os, int n, int m) const
+template <class T>
+void OperatorMP<T>::PrintMatlab(std::ostream & os, int n, int m) const
 {
    using namespace std;
    if (n == 0) { n = width; }
    if (m == 0) { m = height; }
 
-   Vector x(n), y(m);
+   VectorMP<T> x(n), y(m);
    x = 0.0;
 
    os << setiosflags(ios::scientific | ios::showpos);
@@ -273,7 +296,8 @@ void Operator::PrintMatlab(std::ostream & os, int n, int m) const
    }
 }
 
-void Operator::PrintMatlab(std::ostream &os) const
+template <class T>
+void OperatorMP<T>::PrintMatlab(std::ostream &os) const
 {
    PrintMatlab(os, width, height);
 }
@@ -404,9 +428,11 @@ SumOperator::~SumOperator()
    if (ownB) { delete B; }
 }
 
-ProductOperator::ProductOperator(const Operator *A, const Operator *B,
-                                 bool ownA, bool ownB)
-   : Operator(A->Height(), B->Width()),
+template <class T>
+ProductOperatorMP<T>::ProductOperatorMP(const OperatorMP<T> *A,
+                                        const OperatorMP<T> *B,
+                                        bool ownA, bool ownB)
+   : OperatorMP<T>(A->Height(), B->Width()),
      A(A), B(B), ownA(ownA), ownB(ownB), z(A->Width())
 {
    MFEM_VERIFY(A->Width() == B->Height(),
@@ -423,16 +449,18 @@ ProductOperator::ProductOperator(const Operator *A, const Operator *B,
    }
 }
 
-ProductOperator::~ProductOperator()
+template <class T>
+ProductOperatorMP<T>::~ProductOperatorMP()
 {
    if (ownA) { delete A; }
    if (ownB) { delete B; }
 }
 
-
-RAPOperator::RAPOperator(const Operator &Rt_, const Operator &A_,
-                         const Operator &P_)
-   : Operator(Rt_.Width(), P_.Width()), Rt(Rt_), A(A_), P(P_)
+template <class T>
+RAPOperatorMP<T>::RAPOperatorMP(const OperatorMP<T> &Rt_,
+                                const OperatorMP<T> &A_,
+                                const OperatorMP<T> &P_)
+   : OperatorMP<T>(Rt_.Width(), P_.Width()), Rt(Rt_), A(A_), P(P_)
 {
    MFEM_VERIFY(Rt.Height() == A.Height(),
                "incompatible Operators: Rt.Height() = " << Rt.Height()
@@ -463,11 +491,11 @@ RAPOperator::RAPOperator(const Operator &Rt_, const Operator &A_,
    APx.SetSize(A.Height(), mem_type);
 }
 
-
-TripleProductOperator::TripleProductOperator(
-   const Operator *A, const Operator *B, const Operator *C,
+template <class T>
+TripleProductOperatorMP<T>::TripleProductOperatorMP(
+   const OperatorMP<T> *A, const OperatorMP<T> *B, const OperatorMP<T> *C,
    bool ownA, bool ownB, bool ownC)
-   : Operator(A->Height(), C->Width())
+   : OperatorMP<T>(A->Height(), C->Width())
    , A(A), B(B), C(C)
    , ownA(ownA), ownB(ownB), ownC(ownC)
 {
@@ -500,18 +528,20 @@ TripleProductOperator::TripleProductOperator(
    t2.SetSize(B->Height(), mem_type);
 }
 
-TripleProductOperator::~TripleProductOperator()
+template <class T>
+TripleProductOperatorMP<T>::~TripleProductOperatorMP()
 {
    if (ownA) { delete A; }
    if (ownB) { delete B; }
    if (ownC) { delete C; }
 }
 
-
-ConstrainedOperator::ConstrainedOperator(Operator *A, const Array<int> &list,
-                                         bool own_A_,
-                                         DiagonalPolicy diag_policy_)
-   : Operator(A->Height(), A->Width()), A(A), own_A(own_A_),
+template <class T>
+ConstrainedOperatorMP<T>::ConstrainedOperatorMP(OperatorMP<T> *A,
+                                                const Array<int> &list,
+                                                bool own_A_,
+                                                DiagonalPolicy diag_policy_)
+   : OperatorMP<T>(A->Height(), A->Width()), A(A), own_A(own_A_),
      diag_policy(diag_policy_)
 {
    // 'mem_class' should work with A->Mult() and mfem::forall():
@@ -521,11 +551,12 @@ ConstrainedOperator::ConstrainedOperator(Operator *A, const Array<int> &list,
    constraint_list.MakeRef(list);
    // typically z and w are large vectors, so use the device (GPU) to perform
    // operations on them
-   z.SetSize(height, mem_type); z.UseDevice(true);
-   w.SetSize(height, mem_type); w.UseDevice(true);
+   z.SetSize(this->height, mem_type); z.UseDevice(true);
+   w.SetSize(this->height, mem_type); w.UseDevice(true);
 }
 
-void ConstrainedOperator::AssembleDiagonal(Vector &diag) const
+template <class T>
+void ConstrainedOperatorMP<T>::AssembleDiagonal(VectorMP<T> &diag) const
 {
    A->AssembleDiagonal(diag);
 
@@ -556,7 +587,9 @@ void ConstrainedOperator::AssembleDiagonal(Vector &diag) const
    }
 }
 
-void ConstrainedOperator::EliminateRHS(const Vector &x, Vector &b) const
+template <class T>
+void ConstrainedOperatorMP<T>::EliminateRHS(const VectorMP<T> &x,
+                                            VectorMP<T> &b) const
 {
    w = 0.0;
    const int csz = constraint_list.Size();
@@ -583,8 +616,10 @@ void ConstrainedOperator::EliminateRHS(const Vector &x, Vector &b) const
    });
 }
 
-void ConstrainedOperator::ConstrainedMult(const Vector &x, Vector &y,
-                                          const bool transpose) const
+template <class T>
+void ConstrainedOperatorMP<T>::ConstrainedMult(const VectorMP<T> &x,
+                                               VectorMP<T> &y,
+                                               const bool transpose) const
 {
    const int csz = constraint_list.Size();
    if (csz == 0)
@@ -645,8 +680,10 @@ void ConstrainedOperator::ConstrainedMult(const Vector &x, Vector &y,
    }
 }
 
-void ConstrainedOperator::ConstrainedAbsMult(const Vector &x, Vector &y,
-                                             const bool transpose) const
+template <class T>
+void ConstrainedOperatorMP<T>::ConstrainedAbsMult(const VectorMP<T> &x,
+                                                  VectorMP<T> &y,
+                                                  const bool transpose) const
 {
    const int csz = constraint_list.Size();
    if (csz == 0)
@@ -707,43 +744,52 @@ void ConstrainedOperator::ConstrainedAbsMult(const Vector &x, Vector &y,
    }
 }
 
-void ConstrainedOperator::Mult(const Vector &x, Vector &y) const
+template <class T>
+void ConstrainedOperatorMP<T>::Mult(const VectorMP<T> &x, VectorMP<T> &y) const
 {
    constexpr bool transpose = false;
    ConstrainedMult(x, y, transpose);
 }
 
-void ConstrainedOperator::AbsMult(const Vector &x, Vector &y) const
+template <class T>
+void ConstrainedOperatorMP<T>::AbsMult(const VectorMP<T> &x,
+                                       VectorMP<T> &y) const
 {
    constexpr bool transpose = false;
    ConstrainedAbsMult(x, y, transpose);
 }
 
-void ConstrainedOperator::MultTranspose(const Vector &x, Vector &y) const
+template <class T>
+void ConstrainedOperatorMP<T>::MultTranspose(const VectorMP<T> &x,
+                                             VectorMP<T> &y) const
 {
    constexpr bool transpose = true;
    ConstrainedMult(x, y, transpose);
 }
 
-void ConstrainedOperator::AbsMultTranspose(const Vector &x, Vector &y) const
+template <class T>
+void ConstrainedOperatorMP<T>::AbsMultTranspose(const VectorMP<T> &x,
+                                                VectorMP<T> &y) const
 {
    constexpr bool transpose = true;
    ConstrainedAbsMult(x, y, transpose);
 }
 
-void ConstrainedOperator::AddMult(const Vector &x, Vector &y,
-                                  const real_t a) const
+template <class T>
+void ConstrainedOperatorMP<T>::AddMult(const VectorMP<T> &x, VectorMP<T> &y,
+                                       const T a) const
 {
    Mult(x, w);
    y.Add(a, w);
 }
 
-RectangularConstrainedOperator::RectangularConstrainedOperator(
-   Operator *A,
+template <class T>
+RectangularConstrainedOperatorMP<T>::RectangularConstrainedOperatorMP(
+   OperatorMP<T> *A,
    const Array<int> &trial_list,
    const Array<int> &test_list,
    bool own_A_)
-   : Operator(A->Height(), A->Width()), A(A), own_A(own_A_)
+   : OperatorMP<T>(A->Height(), A->Width()), A(A), own_A(own_A_)
 {
    // 'mem_class' should work with A->Mult() and mfem::forall():
    mem_class = A->GetMemoryClass()*Device::GetMemoryClass();
@@ -753,12 +799,13 @@ RectangularConstrainedOperator::RectangularConstrainedOperator(
    trial_constraints.MakeRef(trial_list);
    test_constraints.MakeRef(test_list);
    // typically z and w are large vectors, so store them on the device
-   z.SetSize(height, mem_type); z.UseDevice(true);
-   w.SetSize(width, mem_type); w.UseDevice(true);
+   z.SetSize(this->height, mem_type); z.UseDevice(true);
+   w.SetSize(this->width, mem_type); w.UseDevice(true);
 }
 
-void RectangularConstrainedOperator::EliminateRHS(const Vector &x,
-                                                  Vector &b) const
+template <class T>
+void RectangularConstrainedOperatorMP<T>::EliminateRHS(const VectorMP<T> &x,
+                                                       VectorMP<T> &b) const
 {
    w = 0.0;
    const int trial_csz = trial_constraints.Size();
@@ -783,7 +830,9 @@ void RectangularConstrainedOperator::EliminateRHS(const Vector &x,
    });
 }
 
-void RectangularConstrainedOperator::Mult(const Vector &x, Vector &y) const
+template <class T>
+void RectangularConstrainedOperatorMP<T>::Mult(const VectorMP<T> &x,
+                                               VectorMP<T> &y) const
 {
    const int trial_csz = trial_constraints.Size();
    const int test_csz = test_constraints.Size();
@@ -817,8 +866,9 @@ void RectangularConstrainedOperator::Mult(const Vector &x, Vector &y) const
    }
 }
 
-void RectangularConstrainedOperator::MultTranspose(const Vector &x,
-                                                   Vector &y) const
+template <class T>
+void RectangularConstrainedOperatorMP<T>::MultTranspose(const VectorMP<T> &x,
+                                                        VectorMP<T> &y) const
 {
    const int trial_csz = trial_constraints.Size();
    const int test_csz = test_constraints.Size();
@@ -852,7 +902,9 @@ void RectangularConstrainedOperator::MultTranspose(const Vector &x,
    }
 }
 
-real_t InnerProductOperator::Dot(const Vector &x, const Vector &y) const
+template <class T>
+T InnerProductOperatorMP<T>::Dot(const VectorMP<T> &x,
+                                 const VectorMP<T> &y) const
 {
 #ifndef MFEM_USE_MPI
    return (x * y);
@@ -927,4 +979,24 @@ real_t PowerMethod::EstimateLargestEigenvalue(Operator& opr, Vector& v0,
    return eigenvalue;
 }
 
+template class OperatorMP<float>;
+template class OperatorMP<double>;
+
+template class ConstrainedOperatorMP<float>;
+template class ConstrainedOperatorMP<double>;
+
+template class RectangularConstrainedOperatorMP<float>;
+template class RectangularConstrainedOperatorMP<double>;
+
+template class RAPOperatorMP<float>;
+template class RAPOperatorMP<double>;
+
+template class ProductOperatorMP<float>;
+template class ProductOperatorMP<double>;
+
+template class TripleProductOperatorMP<float>;
+template class TripleProductOperatorMP<double>;
+
+template class InnerProductOperatorMP<float>;
+template class InnerProductOperatorMP<double>;
 }

linalg/operator.hpp

@@ -17,32 +17,44 @@
 namespace mfem
 {
 
-class ConstrainedOperator;
-class RectangularConstrainedOperator;
+template <class T>
+class ConstrainedOperatorMP;
 
-/// Abstract operator
-class Operator
+template <class T>
+class RectangularConstrainedOperatorMP;
+
+class OperatorBase
 {
 protected:
    int height; ///< Dimension of the output / number of rows in the matrix.
    int width;  ///< Dimension of the input / number of columns in the matrix.
 
-   /// see FormSystemOperator()
-   /** @note Uses DiagonalPolicy::DIAG_ONE. */
-   void FormConstrainedSystemOperator(
-      const Array<int> &ess_tdof_list, ConstrainedOperator* &Aout);
-
-   /// see FormRectangularSystemOperator()
-   void FormRectangularConstrainedSystemOperator(
-      const Array<int> &trial_tdof_list,
-      const Array<int> &test_tdof_list,
-      RectangularConstrainedOperator* &Aout);
-
-   /** @brief Returns RAP Operator of this, using input/output Prolongation matrices
-       @a Pi corresponds to "P", @a Po corresponds to "Rt" */
-   Operator *SetupRAP(const Operator *Pi, const Operator *Po);
-
 public:
+   /// Get the height (size of output) of the Operator. Synonym with NumRows().
+   inline int Height() const { return height; }
+
+   /// Get the width (size of input) of the Operator. Synonym with NumCols().
+   inline int Width() const { return width; }
+
+   enum Type
+   {
+      ANY_TYPE,         ///< ID for the base class Operator, i.e. any type.
+      MFEM_SPARSEMAT,   ///< ID for class SparseMatrix.
+      Hypre_ParCSR,     ///< ID for class HypreParMatrix.
+      PETSC_MATAIJ,     ///< ID for class PetscParMatrix, MATAIJ format.
+      PETSC_MATIS,      ///< ID for class PetscParMatrix, MATIS format.
+      PETSC_MATSHELL,   ///< ID for class PetscParMatrix, MATSHELL format.
+      PETSC_MATNEST,    ///< ID for class PetscParMatrix, MATNEST format.
+      PETSC_MATHYPRE,   ///< ID for class PetscParMatrix, MATHYPRE format.
+      PETSC_MATGENERIC, ///< ID for class PetscParMatrix, unspecified format.
+      Complex_Operator, ///< ID for class ComplexOperator.
+      MFEM_ComplexSparseMat, ///< ID for class ComplexSparseMatrix.
+      Complex_Hypre_ParCSR,   ///< ID for class ComplexHypreParMatrix.
+      Complex_DenseMat,  ///< ID for class ComplexDenseMatrix
+      MFEM_Block_Matrix,     ///< ID for class BlockMatrix.
+      MFEM_Block_Operator   ///< ID for the base class BlockOperator.
+   };
+
    /// Defines operator diagonal policy upon elimination of rows and/or columns.
    enum DiagonalPolicy
    {
@@ -50,26 +62,45 @@ public:
       DIAG_ONE,  ///< Set the diagonal value to one
       DIAG_KEEP  ///< Keep the diagonal value
    };
+};
 
+/// Abstract operator
+template <class T>
+class OperatorMP : public OperatorBase
+{
+protected:
+   /// see FormSystemOperator()
+   /** @note Uses DiagonalPolicy::DIAG_ONE. */
+   void FormConstrainedSystemOperator(
+      const Array<int> &ess_tdof_list, ConstrainedOperatorMP<T>* &Aout);
+
+   /// see FormRectangularSystemOperator()
+   void FormRectangularConstrainedSystemOperator(
+      const Array<int> &trial_tdof_list,
+      const Array<int> &test_tdof_list,
+      RectangularConstrainedOperatorMP<T>* &Aout);
+
+   /** @brief Returns RAP Operator of this, using input/output Prolongation matrices
+       @a Pi corresponds to "P", @a Po corresponds to "Rt" */
+   OperatorMP *SetupRAP(const OperatorMP<T> *Pi, const OperatorMP<T> *Po);
+
+public:
    /// Initializes memory for true vectors of linear system
-   void InitTVectors(const Operator *Po, const Operator *Ri, const Operator *Pi,
-                     Vector &x, Vector &b, Vector &X, Vector &B) const;
+   void InitTVectors(const OperatorMP<T> *Po, const OperatorMP<T> *Ri,
+                     const OperatorMP<T> *Pi,
+                     VectorMP<T> &x, VectorMP<T> &b, VectorMP<T> &X, VectorMP<T> &B) const;
 
    /// Construct a square Operator with given size s (default 0).
-   explicit Operator(int s = 0) { height = width = s; }
+   explicit OperatorMP(int s = 0) { height = width = s; }
 
    /** @brief Construct an Operator with the given height (output size) and
        width (input size). */
-   Operator(int h, int w) { height = h; width = w; }
+   OperatorMP(int h, int w) { height = h; width = w; }
 
-   /// Get the height (size of output) of the Operator. Synonym with NumRows().
-   inline int Height() const { return height; }
    /** @brief Get the number of rows (size of output) of the Operator. Synonym
        with Height(). */
    inline int NumRows() const { return height; }
 
-   /// Get the width (size of input) of the Operator. Synonym with NumCols().
-   inline int Width() const { return width; }
    /** @brief Get the number of columns (size of input) of the Operator. Synonym
        with Width(). */
    inline int NumCols() const { return width; }
@@ -86,61 +117,63 @@ public:
    virtual MemoryClass GetMemoryClass() const { return MemoryClass::HOST; }
 
    /// Operator application: `y=A(x)`.
-   virtual void Mult(const Vector &x, Vector &y) const = 0;
+   virtual void Mult(const VectorMP<T> &x, VectorMP<T> &y) const = 0;
 
    /** @brief Action of the absolute-value operator: `y=|A|(x)`. The default
      behavior in class Operator is to generate an error. If the Operator is a
      composition of several operators, the composition unfold into a product
      of absolute-value operators too. */
-   virtual void AbsMult(const Vector &x, Vector &y) const
+   virtual void AbsMult(const VectorMP<T> &x, VectorMP<T> &y) const
    { MFEM_ABORT("Operator::AbsMult() is not overridden!"); }
 
    /** @brief Action of the transpose operator: `y=A^t(x)`. The default behavior
        in class Operator is to generate an error. */
-   virtual void MultTranspose(const Vector &x, Vector &y) const
+   virtual void MultTranspose(const VectorMP<T> &x, VectorMP<T> &y) const
    { MFEM_ABORT("Operator::MultTranspose() is not overridden!"); }
 
    /** @brief Action of the transpose absolute-value operator: `y=|A|^t(x)`.
       The default behavior in class Operator is to generate an error. */
-   virtual void AbsMultTranspose(const Vector &x, Vector &y) const
+   virtual void AbsMultTranspose(const VectorMP<T> &x, VectorMP<T> &y) const
    { MFEM_ABORT("Operator::AbsMultTranspose() is not overridden!"); }
 
    /// Operator application: `y+=A(x)` (default) or `y+=a*A(x)`.
-   virtual void AddMult(const Vector &x, Vector &y, const real_t a = 1.0) const;
+   virtual void AddMult(const VectorMP<T> &x, VectorMP<T> &y,
+                        const T a = 1.0) const;
 
    /// Operator transpose application: `y+=A^t(x)` (default) or `y+=a*A^t(x)`.
-   virtual void AddMultTranspose(const Vector &x, Vector &y,
-                                 const real_t a = 1.0) const;
+   virtual void AddMultTranspose(const VectorMP<T> &x, VectorMP<T> &y,
+                                 const T a = 1.0) const;
 
    /// Operator application on a matrix: `Y=A(X)`.
-   virtual void ArrayMult(const Array<const Vector *> &X,
-                          Array<Vector *> &Y) const;
+   virtual void ArrayMult(const Array<const VectorMP<T> *> &X,
+                          Array<VectorMP<T> *> &Y) const;
 
    /// Action of the transpose operator on a matrix: `Y=A^t(X)`.
-   virtual void ArrayMultTranspose(const Array<const Vector *> &X,
-                                   Array<Vector *> &Y) const;
+   virtual void ArrayMultTranspose(const Array<const VectorMP<T> *> &X,
+                                   Array<VectorMP<T> *> &Y) const;
 
    /// Operator application on a matrix: `Y+=A(X)` (default) or `Y+=a*A(X)`.
-   virtual void ArrayAddMult(const Array<const Vector *> &X, Array<Vector *> &Y,
-                             const real_t a = 1.0) const;
+   virtual void ArrayAddMult(const Array<const VectorMP<T> *> &X,
+                             Array<VectorMP<T> *> &Y,
+                             const T a = 1.0) const;
 
    /** @brief Operator transpose application on a matrix: `Y+=A^t(X)` (default)
        or `Y+=a*A^t(X)`. */
-   virtual void ArrayAddMultTranspose(const Array<const Vector *> &X,
-                                      Array<Vector *> &Y, const real_t a = 1.0) const;
+   virtual void ArrayAddMultTranspose(const Array<const VectorMP<T> *> &X,
+                                      Array<VectorMP<T> *> &Y, const T a = 1.0) const;
 
    /** @brief Evaluate the gradient operator at the point @a x. The default
        behavior in class Operator is to generate an error. */
-   virtual Operator &GetGradient(const Vector &x) const
+   virtual OperatorMP<T> &GetGradient(const VectorMP<T> &x) const
    {
       MFEM_ABORT("Operator::GetGradient() is not overridden!");
-      return const_cast<Operator &>(*this);
+      return const_cast<OperatorMP<T> &>(*this);
    }
 
    /** @brief Computes the diagonal entries into @a diag. Typically, this
        operation only makes sense for linear Operator%s. In some cases, only an
        approximation of the diagonal is computed. */
-   virtual void AssembleDiagonal(Vector &diag) const
+   virtual void AssembleDiagonal(VectorMP<T> &diag) const
    {
       MFEM_CONTRACT_VAR(diag);
       MFEM_ABORT("Not relevant or not implemented for this Operator.");
@@ -148,15 +181,15 @@ public:
 
    /** @brief Prolongation operator from linear algebra (linear system) vectors,
        to input vectors for the operator. `NULL` means identity. */
-   virtual const Operator *GetProlongation() const { return NULL; }
+   virtual const OperatorMP<T> *GetProlongation() const { return NULL; }
 
    /** @brief Restriction operator from input vectors for the operator to linear
        algebra (linear system) vectors. `NULL` means identity. */
-   virtual const Operator *GetRestriction() const  { return NULL; }
+   virtual const OperatorMP<T> *GetRestriction() const  { return NULL; }
 
    /** @brief Prolongation operator from linear algebra (linear system) vectors,
        to output vectors for the operator. `NULL` means identity. */
-   virtual const Operator *GetOutputProlongation() const
+   virtual const OperatorMP<T> *GetOutputProlongation() const
    {
       return GetProlongation(); // Assume square unless specialized
    }
@@ -165,11 +198,11 @@ public:
        form to facilitate matrix-free RAP-type operators.
 
        `NULL` means identity. */
-   virtual const Operator *GetOutputRestrictionTranspose() const { return NULL; }
+   virtual const OperatorMP<T> *GetOutputRestrictionTranspose() const { return NULL; }
 
    /** @brief Restriction operator from output vectors for the operator to linear
        algebra (linear system) vectors. `NULL` means identity. */
-   virtual const Operator *GetOutputRestriction() const
+   virtual const OperatorMP<T> *GetOutputRestriction() const
    {
       return GetRestriction(); // Assume square unless specialized
    }
@@ -205,8 +238,8 @@ public:
        @note If there are no transformations, @a X simply reuses the data of @a
        x. */
    void FormLinearSystem(const Array<int> &ess_tdof_list,
-                         Vector &x, Vector &b,
-                         Operator* &A, Vector &X, Vector &B,
+                         VectorMP<T> &x, VectorMP<T> &b,
+                         OperatorMP<T>* &A, VectorMP<T> &X, VectorMP<T> &B,
                          int copy_interior = 0);
 
    /** @brief Form a column-constrained linear system using a matrix-free approach.
@@ -237,8 +270,8 @@ public:
        x. */
    void FormRectangularLinearSystem(const Array<int> &trial_tdof_list,
                                     const Array<int> &test_tdof_list,
-                                    Vector &x, Vector &b,
-                                    Operator* &A, Vector &X, Vector &B);
+                                    VectorMP<T> &x, VectorMP<T> &b,
+                                    OperatorMP<T>* &A, VectorMP<T> &X, VectorMP<T> &B);
 
    /** @brief Reconstruct a solution vector @a x (e.g. a GridFunction) from the
        solution @a X of a constrained linear system obtained from
@@ -249,7 +282,8 @@ public:
        @a x, for this Operator (presumably a finite element grid function). This
        method has identical signature to the analogous method for bilinear
        forms, though currently @a b is not used in the implementation. */
-   virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x);
+   virtual void RecoverFEMSolution(const VectorMP<T> &X, const VectorMP<T> &b,
+                                   VectorMP<T> &x);
 
    /** @brief Return in @a A a parallel (on truedofs) version of this square
        operator.
@@ -257,7 +291,7 @@ public:
        This returns the same operator as FormLinearSystem(), but does without
        the transformations of the right-hand side and initial guess. */
    void FormSystemOperator(const Array<int> &ess_tdof_list,
-                           Operator* &A);
+                           OperatorMP<T>* &A);
 
    /** @brief Return in @a A a parallel (on truedofs) version of this
        rectangular operator (including constraints).
@@ -266,7 +300,7 @@ public:
        without the transformations of the right-hand side. */
    void FormRectangularSystemOperator(const Array<int> &trial_tdof_list,
                                       const Array<int> &test_tdof_list,
-                                      Operator* &A);
+                                      OperatorMP<T>* &A);
 
    /** @brief Return in @a A a parallel (on truedofs) version of this
        rectangular operator.
@@ -279,7 +313,7 @@ public:
        Operator maps between. These are e.g. available through the (parallel)
        finite element space of any (parallel) bilinear form operator. We have:
        `A(X)=[Rout (*this) Pin](X)`. */
-   void FormDiscreteOperator(Operator* &A);
+   void FormDiscreteOperator(OperatorMP<T>* &A);
 
    /// Prints operator with input size n and output size m in Matlab format.
    void PrintMatlab(std::ostream & out, int n, int m = 0) const;
@@ -288,28 +322,7 @@ public:
    virtual void PrintMatlab(std::ostream & out) const;
 
    /// Virtual destructor.
-   virtual ~Operator() { }
-
-   /// Enumeration defining IDs for some classes derived from Operator.
-   /** This enumeration is primarily used with class OperatorHandle. */
-   enum Type
-   {
-      ANY_TYPE,         ///< ID for the base class Operator, i.e. any type.
-      MFEM_SPARSEMAT,   ///< ID for class SparseMatrix.
-      Hypre_ParCSR,     ///< ID for class HypreParMatrix.
-      PETSC_MATAIJ,     ///< ID for class PetscParMatrix, MATAIJ format.
-      PETSC_MATIS,      ///< ID for class PetscParMatrix, MATIS format.
-      PETSC_MATSHELL,   ///< ID for class PetscParMatrix, MATSHELL format.
-      PETSC_MATNEST,    ///< ID for class PetscParMatrix, MATNEST format.
-      PETSC_MATHYPRE,   ///< ID for class PetscParMatrix, MATHYPRE format.
-      PETSC_MATGENERIC, ///< ID for class PetscParMatrix, unspecified format.
-      Complex_Operator, ///< ID for class ComplexOperator.
-      MFEM_ComplexSparseMat, ///< ID for class ComplexSparseMatrix.
-      Complex_Hypre_ParCSR,   ///< ID for class ComplexHypreParMatrix.
-      Complex_DenseMat,  ///< ID for class ComplexDenseMatrix
-      MFEM_Block_Matrix,     ///< ID for class BlockMatrix.
-      MFEM_Block_Operator   ///< ID for the base class BlockOperator.
-   };
+   virtual ~OperatorMP() { }
 
    /// Return the type ID of the Operator class.
    /** This method is intentionally non-virtual, so that it returns the ID of
@@ -319,6 +332,7 @@ public:
    Type GetType() const { return ANY_TYPE; }
 };
 
+using Operator = OperatorMP<real_t>;
 
 /// Base abstract class for first order time dependent operators.
 /** Operator of the form: (u,t) -> k(u,t), where k generally solves the
@@ -788,7 +802,8 @@ public:
 
 
 /// Base class for solvers
-class Solver : public Operator
+template <class T>
+class SolverMP : public OperatorMP<T>
 {
 public:
    /// If true, use the second argument of Mult() as an initial guess.
@@ -798,37 +813,42 @@ public:
 
        @warning Use a Boolean expression for the second parameter (not an int)
        to distinguish this call from the general rectangular constructor. */
-   explicit Solver(int s = 0, bool iter_mode = false)
-      : Operator(s) { iterative_mode = iter_mode; }
+   explicit SolverMP(int s = 0, bool iter_mode = false)
+      : OperatorMP<T>(s) { iterative_mode = iter_mode; }
 
    /// Initialize a Solver with height @a h and width @a w.
-   Solver(int h, int w, bool iter_mode = false)
-      : Operator(h, w) { iterative_mode = iter_mode; }
+   SolverMP(int h, int w, bool iter_mode = false)
+      : OperatorMP<T>(h, w) { iterative_mode = iter_mode; }
 
    /// Set/update the solver for the given operator.
-   virtual void SetOperator(const Operator &op) = 0;
+   virtual void SetOperator(const OperatorMP<T> &op) = 0;
 };
 
+using Solver = SolverMP<real_t>;
 
 /// Identity Operator I: x -> x.
-class IdentityOperator : public Operator
+template <class T>
+class IdentityOperatorMP : public OperatorMP<T>
 {
 public:
    /// Create an identity operator of size @a n.
-   explicit IdentityOperator(int n) : Operator(n) { }
+   explicit IdentityOperatorMP(int n) : OperatorMP<T>(n) { }
 
    /// Operator application
-   void Mult(const Vector &x, Vector &y) const override { y = x; }
+   void Mult(const VectorMP<T> &x, VectorMP<T> &y) const override { y = x; }
 
    /// Application of the transpose
-   void MultTranspose(const Vector &x, Vector &y) const override { y = x; }
+   void MultTranspose(const VectorMP<T> &x, VectorMP<T> &y) const override { y = x; }
 };
 
+using IdentityOperator = IdentityOperatorMP<real_t>;
+
 /// Returns true if P is the identity prolongation, i.e. if it is either NULL or
 /// an IdentityOperator.
-inline bool IsIdentityProlongation(const Operator *P)
+template <class T>
+inline bool IsIdentityProlongation(const OperatorMP<T> *P)
 {
-   return !P || dynamic_cast<const IdentityOperator*>(P);
+   return !P || dynamic_cast<const IdentityOperatorMP<T>*>(P);
 }
 
 /// Scaled Operator B: x -> a A(x).
@@ -855,29 +875,32 @@ public:
 
 /** @brief The transpose of a given operator. Switches the roles of the methods
     Mult() and MultTranspose(). */
-class TransposeOperator : public Operator
+template <class T>
+class TransposeOperatorMP : public OperatorMP<T>
 {
 private:
-   const Operator &A;
+   const OperatorMP<T> &A;
 
 public:
    /// Construct the transpose of a given operator @a *a.
-   TransposeOperator(const Operator *a)
-      : Operator(a->Width(), a->Height()), A(*a) { }
+   TransposeOperatorMP(const OperatorMP<T> *a)
+      : OperatorMP<T>(a->Width(), a->Height()), A(*a) { }
 
    /// Construct the transpose of a given operator @a a.
-   TransposeOperator(const Operator &a)
-      : Operator(a.Width(), a.Height()), A(a) { }
+   TransposeOperatorMP(const OperatorMP<T> &a)
+      : OperatorMP<T>(a.Width(), a.Height()), A(a) { }
 
    /// Operator application. Apply the transpose of the original Operator.
-   void Mult(const Vector &x, Vector &y) const override
+   void Mult(const VectorMP<T> &x, VectorMP<T> &y) const override
    { A.MultTranspose(x, y); }
 
    /// Application of the transpose. Apply the original Operator.
-   void MultTranspose(const Vector &x, Vector &y) const override
+   void MultTranspose(const VectorMP<T> &x, VectorMP<T> &y) const override
    { A.Mult(x, y); }
 };
 
+using TransposeOperator = TransposeOperatorMP<real_t>;
+
 /// General linear combination operator: x -> a A(x) + b B(x).
 class SumOperator : public Operator
 {
@@ -902,48 +925,53 @@ public:
 };
 
 /// General product operator: x -> (A*B)(x) = A(B(x)).
-class ProductOperator : public Operator
+template <class T>
+class ProductOperatorMP : public OperatorMP<T>
 {
-   const Operator *A, *B;
+   const OperatorMP<T> *A, *B;
    bool ownA, ownB;
-   mutable Vector z;
+   mutable VectorMP<T> z;
 
 public:
-   ProductOperator(const Operator *A, const Operator *B, bool ownA, bool ownB);
+   ProductOperatorMP(const OperatorMP<T> *A, const OperatorMP<T> *B, bool ownA,
+                     bool ownB);
 
-   void Mult(const Vector &x, Vector &y) const override
+   void Mult(const VectorMP<T> &x, VectorMP<T> &y) const override
    { B->Mult(x, z); A->Mult(z, y); }
 
-   void MultTranspose(const Vector &x, Vector &y) const override
+   void MultTranspose(const VectorMP<T> &x, VectorMP<T> &y) const override
    { A->MultTranspose(x, z); B->MultTranspose(z, y); }
 
-   virtual ~ProductOperator();
+   virtual ~ProductOperatorMP<T>();
 };
 
+using ProductOperator = ProductOperatorMP<real_t>;
 
 /// The operator x -> R*A*P*x constructed through the actions of R^T, A and P
-class RAPOperator : public Operator
+template <class T>
+class RAPOperatorMP : public OperatorMP<T>
 {
 private:
-   const Operator & Rt;
-   const Operator & A;
-   const Operator & P;
-   mutable Vector Px;
-   mutable Vector APx;
+   const OperatorMP<T> & Rt;
+   const OperatorMP<T> & A;
+   const OperatorMP<T> & P;
+   mutable VectorMP<T> Px;
+   mutable VectorMP<T> APx;
    MemoryClass mem_class;
 
 public:
    /// Construct the RAP operator given R^T, A and P.
-   RAPOperator(const Operator &Rt_, const Operator &A_, const Operator &P_);
+   RAPOperatorMP<T>(const OperatorMP<T> &Rt_, const OperatorMP<T> &A_,
+                    const OperatorMP<T> &P_);
 
    MemoryClass GetMemoryClass() const override { return mem_class; }
 
    /// Operator application.
-   void Mult(const Vector & x, Vector & y) const override
+   void Mult(const VectorMP<T> & x, VectorMP<T> & y) const override
    { P.Mult(x, Px); A.Mult(Px, APx); Rt.MultTranspose(APx, y); }
 
    /// Operator-wise absolute-value application.
-   void AbsMult(const Vector & x, Vector & y) const override
+   void AbsMult(const VectorMP<T> & x, VectorMP<T> & y) const override
    { P.AbsMult(x, Px); A.AbsMult(Px, APx); Rt.AbsMultTranspose(APx, y); }
 
    /// Approximate diagonal of the RAP Operator.
@@ -953,7 +981,7 @@ public:
        When P is the FE space prolongation operator on a mesh without hanging
        nodes and Rt = P, the returned diagonal is exact, as long as the diagonal
        of A is also exact. */
-   void AssembleDiagonal(Vector &diag) const override
+   void AssembleDiagonal(VectorMP<T> &diag) const override
    {
       A.AssembleDiagonal(APx);
       P.MultTranspose(APx, diag);
@@ -964,11 +992,11 @@ public:
    }
 
    /// Application of the transpose.
-   void MultTranspose(const Vector & x, Vector & y) const override
+   void MultTranspose(const VectorMP<T> & x, VectorMP<T> & y) const override
    { Rt.Mult(x, APx); A.MultTranspose(APx, Px); P.MultTranspose(Px, y); }
 
    /// Operator-wise absolute-value application of the transpose
-   void AbsMultTranspose(const Vector & x, Vector & y) const override
+   void AbsMultTranspose(const VectorMP<T> & x, VectorMP<T> & y) const override
    {
       Rt.AbsMult(x, APx);
       A.AbsMultTranspose(APx, Px);
@@ -976,32 +1004,35 @@ public:
    }
 };
 
+using RAPOperator = RAPOperatorMP<real_t>;
 
 /// General triple product operator x -> A*B*C*x, with ownership of the factors.
-class TripleProductOperator : public Operator
+template <class T>
+class TripleProductOperatorMP : public OperatorMP<T>
 {
-   const Operator *A;
-   const Operator *B;
-   const Operator *C;
+   const OperatorMP<T> *A;
+   const OperatorMP<T> *B;
+   const OperatorMP<T> *C;
    bool ownA, ownB, ownC;
-   mutable Vector t1, t2;
+   mutable VectorMP<T> t1, t2;
    MemoryClass mem_class;
 
 public:
-   TripleProductOperator(const Operator *A, const Operator *B,
-                         const Operator *C, bool ownA, bool ownB, bool ownC);
+   TripleProductOperatorMP(const OperatorMP<T> *A, const OperatorMP<T> *B,
+                           const OperatorMP<T> *C, bool ownA, bool ownB, bool ownC);
 
    MemoryClass GetMemoryClass() const override { return mem_class; }
 
-   void Mult(const Vector &x, Vector &y) const override
+   void Mult(const VectorMP<T> &x, VectorMP<T> &y) const override
    { C->Mult(x, t1); B->Mult(t1, t2); A->Mult(t2, y); }
 
-   void MultTranspose(const Vector &x, Vector &y) const override
+   void MultTranspose(const VectorMP<T> &x, VectorMP<T> &y) const override
    { A->MultTranspose(x, t2); B->MultTranspose(t2, t1); C->MultTranspose(t1, y); }
 
-   virtual ~TripleProductOperator();
+   virtual ~TripleProductOperatorMP<T>();
 };
 
+using TripleProductOperator = TripleProductOperatorMP<real_t>;
 
 /** @brief Square Operator for imposing essential boundary conditions using only
     the action, Mult(), of a given unconstrained Operator.
@@ -1012,13 +1043,19 @@ public:
 
     Do not confuse with ConstrainedSolver, which despite the name has very
     different functionality. */
-class ConstrainedOperator : public Operator
+template <class T>
+class ConstrainedOperatorMP : public OperatorMP<T>
 {
+   using DiagonalPolicy = OperatorBase::DiagonalPolicy;
+   using OperatorBase::DIAG_ONE;
+   using OperatorBase::DIAG_KEEP;
+   using OperatorBase::DIAG_ZERO;
+
 protected:
    Array<int> constraint_list;  ///< List of constrained indices/dofs.
-   Operator *A;                 ///< The unconstrained Operator.
+   OperatorMP<T> *A;                 ///< The unconstrained Operator.
    bool own_A;                  ///< Ownership flag for A.
-   mutable Vector z, w;         ///< Auxiliary vectors.
+   mutable VectorMP<T> z, w;         ///< Auxiliary vectors.
    MemoryClass mem_class;
    DiagonalPolicy diag_policy;  ///< Diagonal policy for constrained dofs
 
@@ -1031,8 +1068,9 @@ public:
        ownership flag @a own_A is true, the operator @a *A will be destroyed
        when this object is destroyed. The @a diag_policy determines how the
        operator sets entries corresponding to essential dofs. */
-   ConstrainedOperator(Operator *A, const Array<int> &list, bool own_A = false,
-                       DiagonalPolicy diag_policy = DIAG_ONE);
+   ConstrainedOperatorMP(OperatorMP<T> *A, const Array<int> &list,
+                         bool own_A = false,
+                         DiagonalPolicy diag_policy = DIAG_ONE);
 
    /// Returns the type of memory in which the solution and temporaries are stored.
    MemoryClass GetMemoryClass() const override { return mem_class; }
@@ -1042,7 +1080,7 @@ public:
    { diag_policy = diag_policy_; }
 
    /// Diagonal of A, modified according to the used DiagonalPolicy.
-   void AssembleDiagonal(Vector &diag) const override;
+   void AssembleDiagonal(VectorMP<T> &diag) const override;
 
    /** @brief Eliminate "essential boundary condition" values specified in @a x
        from the given right-hand side @a b.
@@ -1055,7 +1093,7 @@ public:
        the vectors, and "_i" -- the rest of the entries.
 
        @note This method is consistent with `DiagonalPolicy::DIAG_ONE`. */
-   void EliminateRHS(const Vector &x, Vector &b) const;
+   void EliminateRHS(const VectorMP<T> &x, VectorMP<T> &b) const;
 
    /** @brief Constrained operator action.
 
@@ -1065,29 +1103,33 @@ public:
 
        where the "_b" subscripts denote the essential (boundary) indices/dofs of
        the vectors, and "_i" -- the rest of the entries. */
-   void Mult(const Vector &x, Vector &y) const override;
+   void Mult(const VectorMP<T> &x, VectorMP<T> &y) const override;
 
-   void AddMult(const Vector &x, Vector &y, const real_t a = 1.0) const override;
+   void AddMult(const VectorMP<T> &x, VectorMP<T> &y,
+                const T a = 1.0) const override;
 
-   void AbsMult(const Vector &x, Vector &y) const override;
+   void AbsMult(const VectorMP<T> &x, VectorMP<T> &y) const override;
 
-   void MultTranspose(const Vector &x, Vector &y) const override;
+   void MultTranspose(const VectorMP<T> &x, VectorMP<T> &y) const override;
 
-   void AbsMultTranspose(const Vector &x, Vector &y) const override;
+   void AbsMultTranspose(const VectorMP<T> &x, VectorMP<T> &y) const override;
 
    /** @brief Implementation of Mult or MultTranspose.
-       TODO - Generalize to allow constraining rows and columns differently. */
-   void ConstrainedMult(const Vector &x, Vector &y, const bool transpose) const;
+    *  TODO - Generalize to allow constraining rows and columns differently. */
+   void ConstrainedMult(const VectorMP<T> &x, VectorMP<T> &y,
+                        const bool transpose) const;
 
    /** @brief Implementation of AbsMult or AbsMultTranspose.
        TODO - Generalize to allow constraining rows and columns differently. */
-   void ConstrainedAbsMult(const Vector &x, Vector &y,
+   void ConstrainedAbsMult(const VectorMP<T> &x, VectorMP<T> &y,
                            const bool transpose) const;
 
    /// Destructor: destroys the unconstrained Operator, if owned.
-   ~ConstrainedOperator() override { if (own_A) { delete A; } }
+   ~ConstrainedOperatorMP<T>() override { if (own_A) { delete A; } }
 };
 
+using ConstrainedOperator = ConstrainedOperatorMP<real_t>;
+
 /** @brief Rectangular Operator for imposing essential boundary conditions on
     the input space using only the action, Mult(), of a given unconstrained
     Operator.
@@ -1095,13 +1137,14 @@ public:
     Rectangular operator constrained by fixing certain entries in the solution
     to given "essential boundary condition" values. This class is used by the
     general matrix-free formulation of Operator::FormRectangularLinearSystem. */
-class RectangularConstrainedOperator : public Operator
+template <class T>
+class RectangularConstrainedOperatorMP : public OperatorMP<T>
 {
 protected:
    Array<int> trial_constraints, test_constraints;
-   Operator *A;
+   OperatorMP<T> *A;
    bool own_A;
-   mutable Vector z, w;
+   mutable VectorMP<T> z, w;
    MemoryClass mem_class;
 
 public:
@@ -1112,8 +1155,8 @@ public:
        constrain, i.e. each entry @a trial_list[i] represents an essential trial
        dof. If the ownership flag @a own_A is true, the operator @a *A will be
        destroyed when this object is destroyed. */
-   RectangularConstrainedOperator(Operator *A, const Array<int> &trial_list,
-                                  const Array<int> &test_list, bool own_A = false);
+   RectangularConstrainedOperatorMP(OperatorMP<T> *A, const Array<int> &trial_list,
+                                    const Array<int> &test_list, bool own_A = false);
    /// Returns the type of memory in which the solution and temporaries are stored.
    MemoryClass GetMemoryClass() const override { return mem_class; }
    /** @brief Eliminate columns corresponding to "essential boundary condition"
@@ -1126,7 +1169,7 @@ public:
 
        where the "_b" subscripts denote the essential (boundary) indices and the
        "_j" subscript denotes the essential test indices */
-   void EliminateRHS(const Vector &x, Vector &b) const;
+   void EliminateRHS(const VectorMP<T> &x, VectorMP<T> &b) const;
    /** @brief Rectangular-constrained operator action.
 
        Performs the following steps:
@@ -1136,16 +1179,19 @@ public:
 
        where the "_i" subscripts denote all the nonessential (boundary) trial
        indices and the "_j" subscript denotes the essential test indices */
-   void Mult(const Vector &x, Vector &y) const override;
-   void MultTranspose(const Vector &x, Vector &y) const override;
-   virtual ~RectangularConstrainedOperator() { if (own_A) { delete A; } }
+   void Mult(const VectorMP<T> &x, VectorMP<T> &y) const override;
+   void MultTranspose(const VectorMP<T> &x, VectorMP<T> &y) const override;
+   virtual ~RectangularConstrainedOperatorMP<T>() { if (own_A) { delete A; } }
 };
 
+using RectangularConstrainedOperator = RectangularConstrainedOperatorMP<real_t>;
+
 /** @brief Abstract class for defining inner products. The method Eval()
     must be implemented in derived classes to compute the inner product
     of two vectors according to a specific inner product definition.
 */
-class InnerProductOperator : public Operator
+template <class T>
+class InnerProductOperatorMP : public OperatorMP<T>
 {
 #ifdef MFEM_USE_MPI
 private:
@@ -1153,16 +1199,16 @@ private:
    int dot_prod_type = 0; // 0: local, 1: global
 
 public:
-   InnerProductOperator(MPI_Comm comm_) : Operator(1)
+   InnerProductOperatorMP(MPI_Comm comm_) : OperatorMP<T>(1)
    { comm = comm_; dot_prod_type = 1; }
 #endif
 protected:
    /// @brief Standard global/local $\ell_2$ inner product.
-   virtual real_t Dot(const Vector &x, const Vector &y) const;
+   virtual T Dot(const VectorMP<T> &x, const VectorMP<T> &y) const;
 
 public:
    /// Create an operator of size 1 (scalar).
-   InnerProductOperator() : Operator(1)
+   InnerProductOperatorMP() : OperatorMP<T>(1)
    {
 #ifdef MFEM_USE_MPI
       dot_prod_type = 0;
@@ -1171,7 +1217,7 @@ public:
 
    /// Operator application - not always needed/used but added
    /// to satisfy the abstract base class interface.
-   virtual void Mult(const Vector &x, Vector &y) const override
+   virtual void Mult(const VectorMP<T> &x, VectorMP<T> &y) const override
    {
       MFEM_ABORT("Mult is not implemented.");
    }
@@ -1179,9 +1225,11 @@ public:
    /** @brief Compute the inner product (x,y) of vectors x and y.
               This is an abstract method that must be
               implemented in derived classes. */
-   virtual real_t Eval(const Vector &x, const Vector &y) = 0;
+   virtual real_t Eval(const VectorMP<T> &x, const VectorMP<T> &y) = 0;
 };
 
+using InnerProductOperator = InnerProductOperatorMP<real_t>;
+
 /** @brief PowerMethod helper class to estimate the largest eigenvalue of an
            operator using the iterative power method. */
 class PowerMethod

linalg/solvers.cpp

@@ -26,8 +26,9 @@ namespace mfem
 
 using namespace std;
 
-IterativeSolver::IterativeSolver()
-   : Solver(0, true)
+template <class T>
+IterativeSolverMP<T>::IterativeSolverMP()
+   : SolverMP<T>(0, true)
 {
    oper = NULL;
    prec = NULL;
@@ -41,8 +42,9 @@ IterativeSolver::IterativeSolver()
 
 #ifdef MFEM_USE_MPI
 
-IterativeSolver::IterativeSolver(MPI_Comm comm_)
-   : Solver(0, true)
+template <class T>
+IterativeSolverMP<T>::IterativeSolverMP(MPI_Comm comm_)
+   : SolverMP<T>(0, true)
 {
    oper = NULL;
    prec = NULL;
@@ -55,7 +57,8 @@ IterativeSolver::IterativeSolver(MPI_Comm comm_)
 
 #endif // MFEM_USE_MPI
 
-real_t IterativeSolver::Dot(const Vector &x, const Vector &y) const
+template <class T>
+T IterativeSolverMP<T>::Dot(const VectorMP<T> &x, const VectorMP<T> &y) const
 {
    if (dot_oper) { return dot_oper->Eval(x,y); } // Use custom inner product (if provided)
 
@@ -73,7 +76,8 @@ real_t IterativeSolver::Dot(const Vector &x, const Vector &y) const
 #endif
 }
 
-void IterativeSolver::SetPrintLevel(int print_lvl)
+template <class T>
+void IterativeSolverMP<T>::SetPrintLevel(int print_lvl)
 {
    print_options = FromLegacyPrintLevel(print_lvl);
    int print_level_ = print_lvl;
@@ -94,7 +98,8 @@ void IterativeSolver::SetPrintLevel(int print_lvl)
    print_level = print_level_;
 }
 
-void IterativeSolver::SetPrintLevel(PrintLevel options)
+template <class T>
+void IterativeSolverMP<T>::SetPrintLevel(PrintLevel options)
 {
    print_options = options;
 
@@ -116,7 +121,9 @@ void IterativeSolver::SetPrintLevel(PrintLevel options)
    print_level = derived_print_level;
 }
 
-IterativeSolver::PrintLevel IterativeSolver::FromLegacyPrintLevel(
+template <class T>
+typename IterativeSolverMP<T>::PrintLevel
+IterativeSolverMP<T>::FromLegacyPrintLevel(
    int print_level_)
 {
 #ifdef MFEM_USE_MPI
@@ -151,7 +158,8 @@ IterativeSolver::PrintLevel IterativeSolver::FromLegacyPrintLevel(
    }
 }
 
-int IterativeSolver::GuessLegacyPrintLevel(PrintLevel print_options_)
+template <class T>
+int IterativeSolverMP<T>::GuessLegacyPrintLevel(PrintLevel print_options_)
 {
    if (print_options_.iterations)
    {
@@ -175,25 +183,28 @@ int IterativeSolver::GuessLegacyPrintLevel(PrintLevel print_options_)
    }
 }
 
-void IterativeSolver::SetPreconditioner(Solver &pr)
+template <class T>
+void IterativeSolverMP<T>::SetPreconditioner(SolverMP<T> &pr)
 {
    prec = &pr;
    prec->iterative_mode = false;
 }
 
-void IterativeSolver::SetOperator(const Operator &op)
+template <class T>
+void IterativeSolverMP<T>::SetOperator(const OperatorMP<T> &op)
 {
    oper = &op;
-   height = op.Height();
-   width = op.Width();
+   this->height = op.Height();
+   this->width = op.Width();
    if (prec)
    {
       prec->SetOperator(*oper);
    }
 }
 
-bool IterativeSolver::Monitor(int it, real_t norm, const Vector& r,
-                              const Vector& x, bool final) const
+template <class T>
+bool IterativeSolverMP<T>::Monitor(int it, T norm, const VectorMP<T>& r,
+                                   const VectorMP<T>& x, bool final) const
 {
    if (controller != nullptr)
    {
@@ -851,30 +862,31 @@ void SLI(const Operator &A, Solver &B, const Vector &b, Vector &x,
    sli.Mult(b, x);
 }
 
-
-void CGSolver::UpdateVectors()
+template <class T>
+void CGSolverMP<T>::UpdateVectors()
 {
-   MemoryType mt = GetMemoryType(oper->GetMemoryClass());
+   MemoryType mt = GetMemoryType(this->oper->GetMemoryClass());
 
-   r.SetSize(width, mt);
+   r.SetSize(this->width, mt);
    r.UseDevice(true);
 
-   d.SetSize(width, mt);
+   d.SetSize(this->width, mt);
    d.UseDevice(true);
 
-   z.SetSize(width, mt);
+   z.SetSize(this->width, mt);
    z.UseDevice(true);
 }
 
-void CGSolver::Mult(const Vector &b, Vector &x) const
+template <class T>
+void CGSolverMP<T>::Mult(const VectorMP<T> &b, VectorMP<T> &x) const
 {
    int i;
-   real_t r0, den, nom, nom0, betanom, alpha, beta;
+   T r0, den, nom, nom0, betanom, alpha, beta;
 
    x.UseDevice(true);
-   if (iterative_mode)
+   if (this->iterative_mode)
    {
-      oper->Mult(x, r);
+      this->oper->Mult(x, r);
       subtract(b, r, r); // r = b - A x
    }
    else
@@ -883,56 +895,56 @@ void CGSolver::Mult(const Vector &b, Vector &x) const
       x = 0.0;
    }
 
-   if (prec)
+   if (this->prec)
    {
-      prec->Mult(r, z); // z = B r
+      this->prec->Mult(r, z); // z = B r
       d = z;
    }
    else
    {
       d = r;
    }
-   nom0 = nom = Dot(d, r);
-   if (nom0 >= 0.0) { initial_norm = sqrt(nom0); }
+   nom0 = nom = this->Dot(d, r);
+   if (nom0 >= 0.0) { this->initial_norm = sqrt(nom0); }
    MFEM_VERIFY(IsFinite(nom), "nom = " << nom);
-   if (print_options.iterations || print_options.first_and_last)
+   if (this->print_options.iterations || this->print_options.first_and_last)
    {
       mfem::out << "   Iteration : " << setw(3) << 0 << "  (B r, r) = "
-                << nom << (print_options.first_and_last ? " ...\n" : "\n");
+                << nom << (this->print_options.first_and_last ? " ...\n" : "\n");
    }
 
    if (nom < 0.0)
    {
-      if (print_options.warnings)
+      if (this->print_options.warnings)
       {
          mfem::out << "PCG: The preconditioner is not positive definite. (Br, r) = "
                    << nom << '\n';
       }
       converged = false;
       final_iter = 0;
-      initial_norm = nom;
+      this->initial_norm = nom;
       final_norm = nom;
 
-      Monitor(0, nom, r, x, true);
+      this->Monitor(0, nom, r, x, true);
       return;
    }
-   r0 = std::max(nom*rel_tol*rel_tol, abs_tol*abs_tol);
-   if (Monitor(0, nom, r, x) || nom <= r0)
+   r0 = std::max(nom*this->rel_tol*this->rel_tol, this->abs_tol*this->abs_tol);
+   if (this->Monitor(0, nom, r, x) || nom <= r0)
    {
       converged = true;
       final_iter = 0;
       final_norm = sqrt(nom);
 
-      Monitor(0, nom, r, x, true);
+      this->Monitor(0, nom, r, x, true);
       return;
    }
 
    oper->Mult(d, z);  // z = A d
-   den = Dot(z, d);
+   den = this->Dot(z, d);
    MFEM_VERIFY(IsFinite(den), "den = " << den);
    if (den <= 0.0)
    {
-      if (Dot(d, d) > 0.0 && print_options.warnings)
+      if (this->Dot(d, d) > 0.0 && print_options.warnings)
       {
          mfem::out << "PCG: The operator is not positive definite. (Ad, d) = "
                    << den << '\n';
@@ -943,7 +955,7 @@ void CGSolver::Mult(const Vector &b, Vector &x) const
          final_iter = 0;
          final_norm = sqrt(nom);
 
-         Monitor(0, nom, r, x, true);
+         this->Monitor(0, nom, r, x, true);
          return;
       }
    }
@@ -960,11 +972,11 @@ void CGSolver::Mult(const Vector &b, Vector &x) const
       if (prec)
       {
          prec->Mult(r, z);      //  z = B r
-         betanom = Dot(r, z);
+         betanom = this->Dot(r, z);
       }
       else
       {
-         betanom = Dot(r, r);
+         betanom = this->Dot(r, r);
       }
       MFEM_VERIFY(IsFinite(betanom), "betanom = " << betanom);
       if (betanom < 0.0)
@@ -985,7 +997,7 @@ void CGSolver::Mult(const Vector &b, Vector &x) const
                    << betanom << std::endl;
       }
 
-      if (Monitor(i, betanom, r, x) || betanom <= r0)
+      if (this->Monitor(i, betanom, r, x) || betanom <= r0)
       {
          converged = true;
          final_iter = i;
@@ -1007,11 +1019,11 @@ void CGSolver::Mult(const Vector &b, Vector &x) const
          add(r, beta, d, d);
       }
       oper->Mult(d, z);       //  z = A d
-      den = Dot(d, z);
+      den = this->Dot(d, z);
       MFEM_VERIFY(IsFinite(den), "den = " << den);
       if (den <= 0.0)
       {
-         if (Dot(d, d) > 0.0 && print_options.warnings)
+         if (this->Dot(d, d) > 0.0 && print_options.warnings)
          {
             mfem::out << "PCG: The operator is not positive definite. (Ad, d) = "
                       << den << '\n';
@@ -1046,7 +1058,7 @@ void CGSolver::Mult(const Vector &b, Vector &x) const
 
    final_norm = sqrt(betanom);
 
-   Monitor(final_iter, final_norm, r, x, true);
+   this->Monitor(final_iter, final_norm, r, x, true);
 }
 
 void CG(const Operator &A, const Vector &b, Vector &x,
@@ -1064,22 +1076,35 @@ void CG(const Operator &A, const Vector &b, Vector &x,
    cg.Mult(b, x);
 }
 
-void PCG(const Operator &A, Solver &B, const Vector &b, Vector &x,
+template <class T>
+void PCG(const OperatorMP<T> &A, SolverMP<T> &B, const VectorMP<T> &b,
+         VectorMP<T> &x,
          int print_iter, int max_num_iter,
-         real_t RTOLERANCE, real_t ATOLERANCE)
+         double RTOLERANCE, double ATOLERANCE)
 {
    MFEM_PERF_FUNCTION;
 
-   CGSolver pcg;
+   CGSolverMP<T> pcg;
    pcg.SetPrintLevel(print_iter);
    pcg.SetMaxIter(max_num_iter);
-   pcg.SetRelTol(sqrt(RTOLERANCE));
-   pcg.SetAbsTol(sqrt(ATOLERANCE));
+   pcg.SetRelTol(sqrt((T)RTOLERANCE));
+   pcg.SetAbsTol(sqrt((T)ATOLERANCE));
    pcg.SetOperator(A);
    pcg.SetPreconditioner(B);
    pcg.Mult(b, x);
 }
 
+template
+void PCG<float>(const OperatorMP<float> &A, SolverMP<float> &B,
+                const VectorMP<float> &b, VectorMP<float> &x,
+                int print_iter, int max_num_iter,
+                double RTOLERANCE, double ATOLERANCE);
+
+template
+void PCG<double>(const OperatorMP<double> &A, SolverMP<double> &B,
+                 const VectorMP<double> &b, VectorMP<double> &x,
+                 int print_iter, int max_num_iter,
+                 double RTOLERANCE, double ATOLERANCE);
 
 inline void GeneratePlaneRotation(real_t &dx, real_t &dy,
                                   real_t &cs, real_t &sn)
@@ -1940,7 +1965,7 @@ void MINRESSolver::Mult(const Vector &b, Vector &x) const
       }
       else if (it == 2)
       {
-         add(1./rho1, *z, -rho2/rho1, w1, w0);   // (w0 == 0)
+         add((real_t) 1./rho1, *z, -rho2/rho1, w1, w0);   // (w0 == 0)
       }
       else
       {
@@ -3057,7 +3082,7 @@ void BlockILU::SetOperator(const Operator &op)
    Factorize();
 }
 
-void BlockILU::CreateBlockPattern(const SparseMatrix &A)
+void BlockILU::CreateBlockPattern(const SparseMatrixMP<real_t> &A)
 {
    MFEM_VERIFY(k_fill == 0, "Only block ILU(0) is currently supported.");
    if (A.Height() % block_size != 0)
@@ -4624,4 +4649,10 @@ void NNLSSolver::Solve(const Vector& rhs_lb, const Vector& rhs_ub,
 }
 #endif // MFEM_USE_LAPACK
 
+template class CGSolverMP<float>;
+template class CGSolverMP<double>;
+
+template class IterativeSolverMP<float>;
+template class IterativeSolverMP<double>;
+
 }

linalg/solvers.hpp

@@ -14,6 +14,7 @@
 
 #include "../config/config.hpp"
 #include "densemat.hpp"
+#include "sparsemat.hpp"
 #include "handle.hpp"
 #include <memory>
 
@@ -32,21 +33,26 @@ namespace mfem
 
 class BilinearForm;
 
+
+template <class T>
+class IterativeSolverMP;
+
 /// Abstract base class for an iterative solver controller
-class IterativeSolverController
+template <class T>
+class IterativeSolverControllerMP
 {
 protected:
    /// The last IterativeSolver to which this controller was attached.
-   const class IterativeSolver *iter_solver;
+   const class IterativeSolverMP<T> *iter_solver;
 
    /// In MonitorResidual or MonitorSolution, this member variable can be set
    /// to true to indicate early convergence.
    bool converged = false;
 
 public:
-   IterativeSolverController() : iter_solver(nullptr) {}
+   IterativeSolverControllerMP() : iter_solver(nullptr) {}
 
-   virtual ~IterativeSolverController() {}
+   virtual ~IterativeSolverControllerMP() {}
 
    /// Has the solver converged?
    ///
@@ -61,13 +67,13 @@ public:
    virtual void Reset() { converged = false; }
 
    /// Monitor the solution vector r
-   virtual void MonitorResidual(int it, real_t norm, const Vector &r,
+   virtual void MonitorResidual(int it, T norm, const VectorMP<T> &r,
                                 bool final)
    {
    }
 
    /// Monitor the solution vector x
-   virtual void MonitorSolution(int it, real_t norm, const Vector &x,
+   virtual void MonitorSolution(int it, T norm, const VectorMP<T> &x,
                                 bool final)
    {
    }
@@ -79,15 +85,16 @@ public:
 
    /** @brief This method is invoked by IterativeSolver::SetController(),
        informing the controller which IterativeSolver is using it. */
-   void SetIterativeSolver(const IterativeSolver &solver)
+   void SetIterativeSolver(const IterativeSolverMP<T> &solver)
    { iter_solver = &solver; }
 };
 
 /// Keeping the alias for backward compatibility
-using IterativeSolverMonitor = IterativeSolverController;
+using IterativeSolverMonitor = IterativeSolverControllerMP<real_t>;
 
 /// Abstract base class for iterative solver
-class IterativeSolver : public Solver
+template <class T>
+class IterativeSolverMP : public SolverMP<T>
 {
 public:
    /** @brief Settings for the output behavior of the IterativeSolver.
@@ -142,10 +149,10 @@ private:
 #endif
 
 protected:
-   const Operator *oper;
-   Solver *prec;
-   IterativeSolverController *controller = nullptr;
-   InnerProductOperator *dot_oper = nullptr;
+   const OperatorMP<T> *oper;
+   SolverMP<T> *prec;
+   IterativeSolverControllerMP<T> *controller = nullptr;
+   InnerProductOperatorMP<T> *dot_oper = nullptr;
 
    /// @name Reporting (protected attributes and member functions)
    ///@{
@@ -177,10 +184,10 @@ protected:
    int max_iter;
 
    /// Relative tolerance.
-   real_t rel_tol;
+   T rel_tol;
 
    /// Absolute tolerance.
-   real_t abs_tol;
+   T abs_tol;
 
    ///@}
 
@@ -190,7 +197,7 @@ protected:
 
    mutable int final_iter = -1;
    mutable bool converged = false;
-   mutable real_t initial_norm = -1.0, final_norm = -1.0;
+   mutable T initial_norm = -1.0, final_norm = -1.0;
 
    ///@}
 
@@ -200,23 +207,23 @@ protected:
        @details Overriding this method in a derived class enables a
        custom inner product.
       */
-   virtual real_t Dot(const Vector &x, const Vector &y) const;
+   virtual T Dot(const VectorMP<T> &x, const VectorMP<T> &y) const;
 
    /// Return the inner product norm of @a x, using the inner product defined by Dot()
-   real_t Norm(const Vector &x) const { return sqrt(Dot(x, x)); }
+   T Norm(const VectorMP<T> &x) const { return sqrt(Dot(x, x)); }
 
    /// Indicated if the controller requires an update of the solution
    bool ControllerRequiresUpdate() const { return controller && controller->RequiresUpdatedSolution(); }
 
    /// Monitor both the residual @a r and the solution @a x
-   bool Monitor(int it, real_t norm, const Vector& r, const Vector& x,
+   bool Monitor(int it, T norm, const VectorMP<T>& r, const VectorMP<T>& x,
                 bool final=false) const;
 
 public:
-   IterativeSolver();
+   IterativeSolverMP();
 
 #ifdef MFEM_USE_MPI
-   IterativeSolver(MPI_Comm comm_);
+   IterativeSolverMP(MPI_Comm comm_);
 #endif
 
    /** @name Convergence
@@ -235,8 +242,8 @@ public:
        X depends on the specific iterative solver.
       */
    ///@{
-   void SetRelTol(real_t rtol) { rel_tol = rtol; }
-   void SetAbsTol(real_t atol) { abs_tol = atol; }
+   void SetRelTol(T rtol) { rel_tol = rtol; }
+   void SetAbsTol(T atol) { abs_tol = atol; }
    void SetMaxIter(int max_it) { max_iter = max_it; }
    ///@}
 
@@ -294,20 +301,20 @@ public:
    /// This function returns the norm of the residual (or preconditioned
    /// residual, depending on the solver), computed before the start of the
    /// iteration.
-   real_t GetInitialNorm() const { return initial_norm; }
+   T GetInitialNorm() const { return initial_norm; }
    /// @brief Returns the final residual norm after termination of the solver
    /// during the last call to Mult().
    ///
    /// This function returns the norm of the residual (or preconditioned
    /// residual, depending on the solver), corresponding to the returned
    /// solution.
-   real_t GetFinalNorm() const { return final_norm; }
+   T GetFinalNorm() const { return final_norm; }
    /// @brief Returns the final residual norm after termination of the solver
    /// during the last call to Mult(), divided by the initial residual norm.
    /// Returns -1 if one of these norms is left undefined by the solver.
    ///
    /// @sa GetFinalNorm(), GetInitialNorm()
-   real_t GetFinalRelNorm() const
+   T GetFinalRelNorm() const
    {
       if (final_norm < 0.0 || initial_norm < 0.0) { return -1.0; }
       return final_norm / initial_norm;
@@ -316,20 +323,20 @@ public:
    ///@}
 
    /// This should be called before SetOperator
-   virtual void SetPreconditioner(Solver &pr);
+   virtual void SetPreconditioner(SolverMP<T> &pr);
 
    /// Also calls SetOperator for the preconditioner if there is one
-   void SetOperator(const Operator &op) override;
+   void SetOperator(const OperatorMP<T> &op) override;
 
    /// Set the iterative solver controller
-   void SetController(IterativeSolverController &c)
+   void SetController(IterativeSolverControllerMP<T> &c)
    { controller = &c; c.SetIterativeSolver(*this); }
 
    /// An alias of SetController() for backward compatibility
-   void SetMonitor(IterativeSolverMonitor &m) { SetController(m); }
+   void SetMonitor(IterativeSolverControllerMP<T> &m) { SetController(m); }
 
    /// Set a user-defined inner product operator (not owned)
-   void SetInnerProduct(InnerProductOperator *ipo) { dot_oper = ipo; }
+   void SetInnerProduct(InnerProductOperatorMP<T> *ipo) { dot_oper = ipo; }
 
 #ifdef MFEM_USE_MPI
    /** @brief Return the associated MPI communicator, or MPI_COMM_NULL if no
@@ -339,6 +346,7 @@ public:
 #endif
 };
 
+using IterativeSolver = IterativeSolverMP<real_t>;
 
 /** @brief Inner product operator constrained to a list of indices/dofs.
     The method Eval() computes the inner product of two vectors
@@ -623,37 +631,50 @@ void SLI(const Operator &A, Solver &B, const Vector &b, Vector &x,
 
 
 /// Conjugate gradient method
-class CGSolver : public IterativeSolver
+template <class T>
+class CGSolverMP : public IterativeSolverMP<T>
 {
 protected:
-   mutable Vector r, d, z;
+   mutable VectorMP<T> r, d, z;
 
    void UpdateVectors();
 
+   using IterativeSolverMP<T>::converged;
+   using IterativeSolverMP<T>::final_iter;
+   using IterativeSolverMP<T>::final_norm;
+   using IterativeSolverMP<T>::oper;
+   using IterativeSolverMP<T>::print_options;
+   using IterativeSolverMP<T>::max_iter;
+   using IterativeSolverMP<T>::prec;
+
 public:
-   CGSolver() { }
+   CGSolverMP() { }
 
 #ifdef MFEM_USE_MPI
-   CGSolver(MPI_Comm comm_) : IterativeSolver(comm_) { }
+   CGSolverMP(MPI_Comm comm_) : IterativeSolverMP<T>(comm_) { }
 #endif
 
-   void SetOperator(const Operator &op) override
-   { IterativeSolver::SetOperator(op); UpdateVectors(); }
+   void SetOperator(const OperatorMP<T> &op) override
+   { IterativeSolverMP<T>::SetOperator(op); UpdateVectors(); }
 
    /** @brief Iterative solution of the linear system using the Conjugate
        Gradient method. */
-   void Mult(const Vector &b, Vector &x) const override;
+   void Mult(const VectorMP<T> &b, VectorMP<T> &x) const override;
 };
 
+using CGSolver = CGSolverMP<real_t>;
+
 /// Conjugate gradient method. (tolerances are squared)
 void CG(const Operator &A, const Vector &b, Vector &x,
         int print_iter = 0, int max_num_iter = 1000,
         real_t RTOLERANCE = 1e-12, real_t ATOLERANCE = 1e-24);
 
 /// Preconditioned conjugate gradient method. (tolerances are squared)
-void PCG(const Operator &A, Solver &B, const Vector &b, Vector &x,
+template <class T>
+void PCG(const OperatorMP<T> &A, SolverMP<T> &B, const VectorMP<T> &b,
+         VectorMP<T> &x,
          int print_iter = 0, int max_num_iter = 1000,
-         real_t RTOLERANCE = 1e-12, real_t ATOLERANCE = 1e-24);
+         double RTOLERANCE = 1e-12, double ATOLERANCE = 1e-24);
 
 
 /// GMRES method
@@ -1152,7 +1173,7 @@ public:
 
 private:
    /// Set up the block CSR structure corresponding to a sparse matrix @a A
-   void CreateBlockPattern(const class SparseMatrix &A);
+   void CreateBlockPattern(const class SparseMatrixMP<real_t> &A);
 
    /// Perform the block ILU factorization
    void Factorize();

linalg/sparsemat.hpp

@@ -34,22 +34,26 @@
 namespace mfem
 {
 
+template <class T>
 class
 #if defined(__alignas_is_defined)
-   alignas(real_t)
+   alignas(T)
 #endif
    RowNode
 {
 public:
-   real_t Value;
+   T Value;
    RowNode *Prev;
    int Column;
 };
 
 /// Data type sparse matrix
-class SparseMatrix : public AbstractSparseMatrix
+template <class T>
+class SparseMatrixMP : public AbstractSparseMatrixMP<T>
 {
 protected:
+   using OperatorBase::height;
+   using OperatorBase::width;
    /// @name Arrays used by the CSR storage format.
    /** */
    ///@{
@@ -65,22 +69,22 @@ protected:
    Memory<int> J;
    /** @brief %Array with size #I[#height], containing the actual entries of the
        sparse matrix, as indexed by the #I array. */
-   Memory<real_t> A;
+   Memory<T> A;
    ///@}
 
    /** @brief %Array of linked lists, one for every row. This array represents
        the linked list (LIL) storage format. */
-   RowNode **Rows;
+   RowNode<T> **Rows;
 
    mutable int current_row;
    mutable int* ColPtrJ;
-   mutable RowNode ** ColPtrNode;
+   mutable RowNode<T> ** ColPtrNode;
 
    /// Transpose of A. Owned. Used to perform MultTranspose() on devices.
-   mutable SparseMatrix *At;
+   mutable SparseMatrixMP<T> *At;
 
 #ifdef MFEM_USE_MEMALLOC
-   typedef MemAlloc <RowNode, 1024> RowNodeAlloc;
+   typedef MemAlloc <RowNode<T>, 1024> RowNodeAlloc;
    RowNodeAlloc * NodesMem;
 #endif
 
@@ -136,7 +140,7 @@ protected:
 
 public:
    /// Create an empty SparseMatrix.
-   SparseMatrix()
+   SparseMatrixMP()
    {
       SetEmpty();
 
@@ -148,25 +152,25 @@ public:
    /** New entries are added as needed by methods like AddSubMatrix(),
        SetSubMatrix(), etc. Calling Finalize() will convert the SparseMatrix to
        the more compact compressed sparse row (CSR) format. */
-   explicit SparseMatrix(int nrows, int ncols = -1);
+   explicit SparseMatrixMP(int nrows, int ncols = -1);
 
    /** @brief Create a sparse matrix in CSR format. Ownership of @a i, @a j, and
        @a data is transferred to the SparseMatrix. */
-   SparseMatrix(int *i, int *j, real_t *data, int m, int n);
+   SparseMatrixMP(int *i, int *j, T *data, int m, int n);
 
    /** @brief Create a sparse matrix in CSR format. Ownership of @a i, @a j, and
        @a data is optionally transferred to the SparseMatrix. */
    /** If the parameter @a data is NULL, then the internal #A array is allocated
        by this constructor (initializing it with zeros and taking ownership,
        regardless of the parameter @a owna). */
-   SparseMatrix(int *i, int *j, real_t *data, int m, int n, bool ownij,
-                bool owna, bool issorted);
+   SparseMatrixMP(int *i, int *j, T *data, int m, int n, bool ownij,
+                  bool owna, bool issorted);
 
    /** @brief Create a sparse matrix in CSR format where each row has space
        allocated for exactly @a rowsize entries. */
    /** SetRow() can then be called or the #I, #J, #A arrays can be used
        directly. */
-   SparseMatrix(int nrows, int ncols, int rowsize);
+   SparseMatrixMP(int nrows, int ncols, int rowsize);
 
    /// Copy constructor (deep copy).
    /** If @a mat is finalized and @a copy_graph is false, the #I and #J arrays
@@ -176,11 +180,16 @@ public:
        SparseMatrix's #I, #J, and #A arrays will be the same as @a mat,
        otherwise the type will be @a mt for those arrays that are deep
        copied. */
-   SparseMatrix(const SparseMatrix &mat, bool copy_graph = true,
-                MemoryType mt = MemoryType::PRESERVE);
+   SparseMatrixMP(const SparseMatrixMP<T> &mat, bool copy_graph = true,
+                  MemoryType mt = MemoryType::PRESERVE);
 
    /// Create a SparseMatrix with diagonal @a v, i.e. A = Diag(v)
-   SparseMatrix(const Vector & v);
+   SparseMatrixMP(const VectorMP<T> & v);
+
+   using DiagonalPolicy = OperatorBase::DiagonalPolicy;
+   using OperatorBase::DIAG_ZERO;
+   using OperatorBase::DIAG_ONE;
+   using OperatorBase::DIAG_KEEP;
 
    /// @brief Sets the height and width of the matrix.
    /** @warning This does not modify in any way the underlying CSR or LIL
@@ -202,16 +211,18 @@ public:
    void UseCuSparse(bool useCuSparse_ = true) { UseGPUSparse(useCuSparse_); }
 
    /// Assignment operator: deep copy
-   SparseMatrix& operator=(const SparseMatrix &rhs);
+   SparseMatrixMP<T>& operator=(const SparseMatrixMP<T> &rhs);
 
    /** @brief Clear the contents of the SparseMatrix and make it a reference to
        @a master */
    /** After this call, the matrix will point to the same data as @a master but
        it will not own its data. The @a master must be finalized. */
-   void MakeRef(const SparseMatrix &master);
+   void MakeRef(const SparseMatrixMP<T> &master);
+
+   //using int OperatorBase::Height();
 
    /// For backward compatibility, define Size() to be synonym of Height().
-   int Size() const { return Height(); }
+   int Size() const { return this->Height(); }
 
    /// Clear the contents of the SparseMatrix.
    void Clear() { Destroy(); SetEmpty(); }
@@ -237,25 +248,25 @@ public:
    inline const int *GetJ() const { return J; }
 
    /// Return the element data, i.e. the array #A.
-   inline real_t *GetData() { return A; }
+   inline T *GetData() { return A; }
    /// Return the element data, i.e. the array #A, const version.
-   inline const real_t *GetData() const { return A; }
+   inline const T *GetData() const { return A; }
 
    // Memory access methods for the #I array.
    Memory<int> &GetMemoryI() { return I; }
    const Memory<int> &GetMemoryI() const { return I; }
    const int *ReadI(bool on_dev = true) const
-   { return mfem::Read(I, Height()+1, on_dev); }
+   { return mfem::Read(I, this->Height()+1, on_dev); }
    int *WriteI(bool on_dev = true)
-   { return mfem::Write(I, Height()+1, on_dev); }
+   { return mfem::Write(I, this->Height()+1, on_dev); }
    int *ReadWriteI(bool on_dev = true)
-   { return mfem::ReadWrite(I, Height()+1, on_dev); }
+   { return mfem::ReadWrite(I, this->Height()+1, on_dev); }
    const int *HostReadI() const
-   { return mfem::Read(I, Height()+1, false); }
+   { return mfem::Read(I, this->Height()+1, false); }
    int *HostWriteI()
-   { return mfem::Write(I, Height()+1, false); }
+   { return mfem::Write(I, this->Height()+1, false); }
    int *HostReadWriteI()
-   { return mfem::ReadWrite(I, Height()+1, false); }
+   { return mfem::ReadWrite(I, this->Height()+1, false); }
 
    // Memory access methods for the #J array.
    Memory<int> &GetMemoryJ() { return J; }
@@ -274,19 +285,19 @@ public:
    { return mfem::ReadWrite(J, J.Capacity(), false); }
 
    // Memory access methods for the #A array.
-   Memory<real_t> &GetMemoryData() { return A; }
-   const Memory<real_t> &GetMemoryData() const { return A; }
-   const real_t *ReadData(bool on_dev = true) const
+   Memory<T> &GetMemoryData() { return A; }
+   const Memory<T> &GetMemoryData() const { return A; }
+   const T *ReadData(bool on_dev = true) const
    { return mfem::Read(A, A.Capacity(), on_dev); }
-   real_t *WriteData(bool on_dev = true)
+   T *WriteData(bool on_dev = true)
    { return mfem::Write(A, A.Capacity(), on_dev); }
-   real_t *ReadWriteData(bool on_dev = true)
+   T *ReadWriteData(bool on_dev = true)
    { return mfem::ReadWrite(A, A.Capacity(), on_dev); }
-   const real_t *HostReadData() const
+   const T *HostReadData() const
    { return mfem::Read(A, A.Capacity(), false); }
-   real_t *HostWriteData()
+   T *HostWriteData()
    { return mfem::Write(A, A.Capacity(), false); }
-   real_t *HostReadWriteData()
+   T *HostReadWriteData()
    { return mfem::ReadWrite(A, A.Capacity(), false); }
 
    /// Returns the number of elements in row @a i.
@@ -301,9 +312,9 @@ public:
    const int *GetRowColumns(const int row) const;
 
    /// Return a pointer to the entries in a row.
-   real_t *GetRowEntries(const int row);
+   T *GetRowEntries(const int row);
    /// Return a pointer to the entries in a row, const version.
-   const real_t *GetRowEntries(const int row) const;
+   const T *GetRowEntries(const int row) const;
 
    /// Change the width of a SparseMatrix.
    /*!
@@ -327,16 +338,16 @@ public:
    void MoveDiagonalFirst();
 
    /// Returns reference to a_{ij}.
-   real_t &Elem(int i, int j) override;
+   T &Elem(int i, int j) override;
 
    /// Returns constant reference to a_{ij}.
-   const real_t &Elem(int i, int j) const override;
+   const T &Elem(int i, int j) const override;
 
    /// Returns reference to A[i][j].
-   real_t &operator()(int i, int j);
+   T &operator()(int i, int j);
 
    /// Returns reference to A[i][j].
-   const real_t &operator()(int i, int j) const;
+   const T &operator()(int i, int j) const;
 
    /// Returns the Diagonal of A
    void GetDiag(Vector & d) const;
@@ -354,24 +365,24 @@ public:
    }
 
    /// Matrix vector multiplication.
-   void Mult(const Vector &x, Vector &y) const override;
+   void Mult(const VectorMP<T> &x, VectorMP<T> &y) const override;
 
    /// y += A * x (default)  or  y += a * A * x
-   void AddMult(const Vector &x, Vector &y,
-                const real_t a = 1.0) const override;
+   void AddMult(const VectorMP<T> &x, VectorMP<T> &y,
+                const T a = 1.0) const override;
 
    /// Multiply a vector with the transposed matrix. y = At * x
    /** If the matrix is modified, call ResetTranspose() and optionally
        EnsureMultTranspose() to make sure this method uses the correct updated
        transpose. */
-   void MultTranspose(const Vector &x, Vector &y) const override;
+   void MultTranspose(const VectorMP<T> &x, VectorMP<T> &y) const override;
 
    /// y += At * x (default)  or  y += a * At * x
    /** If the matrix is modified, call ResetTranspose() and optionally
        EnsureMultTranspose() to make sure this method uses the correct updated
        transpose. */
-   void AddMultTranspose(const Vector &x, Vector &y,
-                         const real_t a = 1.0) const override;
+   void AddMultTranspose(const VectorMP<T> &x, VectorMP<T> &y,
+                         const T a = 1.0) const override;
 
    /** @brief Build and store internally the transpose of this matrix which will
        be used in the methods AddMultTranspose(), MultTranspose(), and
@@ -415,7 +426,7 @@ public:
 
    void PartMult(const Array<int> &rows, const Vector &x, Vector &y) const;
    void PartAddMult(const Array<int> &rows, const Vector &x, Vector &y,
-                    const real_t a=1.0) const;
+                    const T a=1.0) const;
 
    /// y = A * x, treating all entries as booleans (zero=false, nonzero=true).
    /** The actual values stored in the data array, #A, are not used - this means
@@ -430,27 +441,27 @@ public:
    void BooleanMultTranspose(const Array<int> &x, Array<int> &y) const;
 
    /// y = |A| * x, using entry-wise absolute values of matrix A
-   void AbsMult(const Vector &x, Vector &y) const override;
+   void AbsMult(const VectorMP<T> &x, VectorMP<T> &y) const override;
 
    /// y = |At| * x, using entry-wise absolute values of the transpose of matrix A
    /** If the matrix is modified, call ResetTranspose() and optionally
        EnsureMultTranspose() to make sure this method uses the correct updated
        transpose. */
-   void AbsMultTranspose(const Vector &x, Vector &y) const override;
+   void AbsMultTranspose(const VectorMP<T> &x, VectorMP<T> &y) const override;
 
    /// Compute y^t A x
-   real_t InnerProduct(const Vector &x, const Vector &y) const;
+   T InnerProduct(const VectorMP<T> &x, const VectorMP<T> &y) const;
 
    /// For all i compute $ x_i = \sum_j A_{ij} $
-   void GetRowSums(Vector &x) const;
+   void GetRowSums(VectorMP<T> &x) const;
    /// For i = irow compute $ x_i = \sum_j | A_{i, j} | $
-   real_t GetRowNorml1(int irow) const;
+   T GetRowNorml1(int irow) const;
 
    /// This virtual method is not supported: it always returns NULL.
-   MatrixInverse *Inverse() const override;
+   MatrixInverseMP<T> *Inverse() const override;
 
    /// Eliminates a column from the transpose matrix.
-   void EliminateRow(int row, const real_t sol, Vector &rhs);
+   void EliminateRow(int row, const T sol, Vector &rhs);
 
    /// Eliminates a row from the matrix.
    /*!
@@ -478,7 +489,7 @@ public:
 
    /** @brief Similar to EliminateCols + save the eliminated entries into
        @a Ae so that (*this) + Ae is equal to the original matrix. */
-   void EliminateCols(const Array<int> &col_marker, SparseMatrix &Ae);
+   void EliminateCols(const Array<int> &col_marker, SparseMatrixMP &Ae);
 
    /// Eliminate row @a rc and column @a rc and modify the @a rhs using @a sol.
    /** Eliminates the column @a rc to the @a rhs, deletes the row @a rc and
@@ -486,7 +497,7 @@ public:
        is assembled if and only if the element (rc,i) is assembled.
        By default, elements (rc,rc) are set to 1.0, although this behavior
        can be adjusted by changing the @a dpolicy parameter. */
-   void EliminateRowCol(int rc, const real_t sol, Vector &rhs,
+   void EliminateRowCol(int rc, const T sol, Vector &rhs,
                         DiagonalPolicy dpolicy = DIAG_ONE);
 
    /** @brief Similar to
@@ -498,7 +509,7 @@ public:
                                    DiagonalPolicy dpolicy = DIAG_ONE);
 
    /// Perform elimination and set the diagonal entry to the given value
-   void EliminateRowColDiag(int rc, real_t value);
+   void EliminateRowColDiag(int rc, T value);
 
    /// Eliminate row @a rc and column @a rc.
    void EliminateRowCol(int rc, DiagonalPolicy dpolicy = DIAG_ONE);
@@ -506,7 +517,7 @@ public:
    /** @brief Similar to EliminateRowCol(int, DiagonalPolicy) + save the
        eliminated entries into @a Ae so that (*this) + Ae is equal to the
        original matrix */
-   void EliminateRowCol(int rc, SparseMatrix &Ae,
+   void EliminateRowCol(int rc, SparseMatrixMP &Ae,
                         DiagonalPolicy dpolicy = DIAG_ONE);
 
    /** @brief Eliminate essential (Dirichlet) boundary conditions.
@@ -520,31 +531,31 @@ public:
    /// If a row contains only one diag entry of zero, set it to 1.
    void SetDiagIdentity();
    /// If a row contains only zeros, set its diagonal to 1.
-   void EliminateZeroRows(const real_t threshold = 1e-12) override;
+   void EliminateZeroRows(const T threshold = 1e-12) override;
 
    /// Gauss-Seidel forward and backward iterations over a vector x.
-   void Gauss_Seidel_forw(const Vector &x, Vector &y) const;
-   void Gauss_Seidel_back(const Vector &x, Vector &y) const;
+   void Gauss_Seidel_forw(const VectorMP<T> &x, VectorMP<T> &y) const;
+   void Gauss_Seidel_back(const VectorMP<T> &x, VectorMP<T> &y) const;
 
    /// Determine appropriate scaling for Jacobi iteration
-   real_t GetJacobiScaling() const;
+   T GetJacobiScaling() const;
    /** One scaled Jacobi iteration for the system A x = b.
        x1 = x0 + sc D^{-1} (b - A x0)  where D is the diag of A.
        Absolute values of D are used when use_abs_diag = true. */
    void Jacobi(const Vector &b, const Vector &x0, Vector &x1,
-               real_t sc, bool use_abs_diag = false) const;
+               T sc, bool use_abs_diag = false) const;
 
    /// x = sc b / A_ii. When use_abs_diag = true, |A_ii| is used.
    void DiagScale(const Vector &b, Vector &x,
-                  real_t sc = 1.0, bool use_abs_diag = false) const;
+                  T sc = 1.0, bool use_abs_diag = false) const;
 
    /** x1 = x0 + sc D^{-1} (b - A x0) where $ D_{ii} = \sum_j |A_{ij}| $. */
    void Jacobi2(const Vector &b, const Vector &x0, Vector &x1,
-                real_t sc = 1.0) const;
+                T sc = 1.0) const;
 
    /** x1 = x0 + sc D^{-1} (b - A x0) where $ D_{ii} = \sum_j A_{ij} $. */
    void Jacobi3(const Vector &b, const Vector &x0, Vector &x1,
-                real_t sc = 1.0) const;
+                T sc = 1.0) const;
 
    /** @brief Finalize the matrix initialization, switching the storage format
        from LIL to CSR. */
@@ -564,12 +575,12 @@ public:
    /** @brief Remove entries smaller in absolute value than a given tolerance
        @a tol. If @a fix_empty_rows is true, a zero value is inserted in the
        diagonal entry (for square matrices only) */
-   void Threshold(real_t tol, bool fix_empty_rows = false);
+   void Threshold(T tol, bool fix_empty_rows = false);
 
    /** Split the matrix into M x N blocks of sparse matrices in CSR format.
        The 'blocks' array is M x N (i.e. M and N are determined by its
        dimensions) and its entries are overwritten by the new blocks. */
-   void GetBlocks(Array2D<SparseMatrix *> &blocks) const;
+   void GetBlocks(Array2D<SparseMatrixMP *> &blocks) const;
 
    void GetSubMatrix(const Array<int> &rows, const Array<int> &cols,
                      DenseMatrix &subm) const;
@@ -588,24 +599,24 @@ public:
        SparseMatrix, it will be added to the sparsity pattern initialized with
        zero. If the matrix is finalized and the entry is not found, an error
        will be generated. */
-   inline real_t &SearchRow(const int col);
+   inline T &SearchRow(const int col);
    /// Add a value to an entry in the "current row". See SetColPtr().
-   inline void _Add_(const int col, const real_t a)
+   inline void _Add_(const int col, const T a)
    { SearchRow(col) += a; }
    /// Set an entry in the "current row". See SetColPtr().
-   inline void _Set_(const int col, const real_t a)
+   inline void _Set_(const int col, const T a)
    { SearchRow(col) = a; }
    /// Read the value of an entry in the "current row". See SetColPtr().
-   inline real_t _Get_(const int col) const;
+   inline T _Get_(const int col) const;
 
-   inline real_t &SearchRow(const int row, const int col);
-   inline void _Add_(const int row, const int col, const real_t a)
+   inline T &SearchRow(const int row, const int col);
+   inline void _Add_(const int row, const int col, const T a)
    { SearchRow(row, col) += a; }
-   inline void _Set_(const int row, const int col, const real_t a)
+   inline void _Set_(const int row, const int col, const T a)
    { SearchRow(row, col) = a; }
 
-   void Set(const int i, const int j, const real_t val);
-   void Add(const int i, const int j, const real_t val);
+   void Set(const int i, const int j, const T val);
+   void Add(const int i, const int j, const T val);
 
    void SetSubMatrix(const Array<int> &rows, const Array<int> &cols,
                      const DenseMatrix &subm, int skip_zeros = 1);
@@ -634,12 +645,12 @@ public:
          when the matrix is finalized.
        @warning This method breaks the const-ness when the matrix is finalized
        because it gives write access to the #J and #A arrays. */
-   int GetRow(const int row, Array<int> &cols, Vector &srow) const override;
+   int GetRow(const int row, Array<int> &cols, VectorMP<T> &srow) const override;
 
    void SetRow(const int row, const Array<int> &cols, const Vector &srow);
    void AddRow(const int row, const Array<int> &cols, const Vector &srow);
 
-   void ScaleRow(const int row, const real_t scale);
+   void ScaleRow(const int row, const T scale);
    /// this = diag(sl) * this;
    void ScaleRows(const Vector & sl);
    /// this = this * diag(sr);
@@ -647,15 +658,15 @@ public:
 
    /** @brief Add the sparse matrix 'B' to '*this'. This operation will cause an
        error if '*this' is finalized and 'B' has larger sparsity pattern. */
-   SparseMatrix &operator+=(const SparseMatrix &B);
+   SparseMatrixMP &operator+=(const SparseMatrixMP &B);
 
    /** @brief Add the sparse matrix 'B' scaled by the scalar 'a' into '*this'.
        Only entries in the sparsity pattern of '*this' are added. */
-   void Add(const real_t a, const SparseMatrix &B);
+   void Add(const T a, const SparseMatrixMP &B);
 
-   SparseMatrix &operator=(real_t a);
+   SparseMatrixMP &operator=(T a);
 
-   SparseMatrix &operator*=(real_t a);
+   SparseMatrixMP &operator*=(T a);
 
    /// Prints matrix to stream out.
    /** @note The host in synchronized when the finalized matrix is on the device. */
@@ -694,7 +705,7 @@ public:
    void PrintInfo(std::ostream &out) const;
 
    /// Returns max_{i,j} |(i,j)-(j,i)| for a finalized matrix
-   real_t IsSymmetric() const;
+   T IsSymmetric() const;
 
    /// (*this) = 1/2 ((*this) + (*this)^t)
    void Symmetrize();
@@ -702,10 +713,10 @@ public:
    /// Returns the number of the nonzero elements in the matrix
    int NumNonZeroElems() const override;
 
-   real_t MaxNorm() const;
+   T MaxNorm() const;
 
    /// Count the number of entries with |a_ij| <= tol.
-   int CountSmallElems(real_t tol) const;
+   int CountSmallElems(T tol) const;
 
    /// Count the number of entries that are NOT finite, i.e. Inf or Nan.
    int CheckFinite() const;
@@ -726,14 +737,18 @@ public:
    /// Lose the ownership of the graph (I, J) and data (A) arrays.
    void LoseData() { SetGraphOwner(false); SetDataOwner(false); }
 
-   void Swap(SparseMatrix &other);
+   void Swap(SparseMatrixMP &other);
 
    /// Destroys sparse matrix.
-   virtual ~SparseMatrix();
+   virtual ~SparseMatrixMP();
 
-   Type GetType() const { return MFEM_SPARSEMAT; }
+   using Type = OperatorBase::Type;
+
+   Type GetType() const { return OperatorBase::MFEM_SPARSEMAT; }
 };
 
+using SparseMatrix = SparseMatrixMP<real_t>;
+
 inline std::ostream& operator<<(std::ostream& os, SparseMatrix const& mat)
 {
    mat.Print(os);
@@ -741,14 +756,17 @@ inline std::ostream& operator<<(std::ostream& os, SparseMatrix const& mat)
 }
 
 /// Applies f() to each element of the matrix (after it is finalized).
-void SparseMatrixFunction(SparseMatrix &S, real_t (*f)(real_t));
+template <class T>
+void SparseMatrixFunction(SparseMatrixMP<T> &S, T (*f)(T));
 
 
 /// Transpose of a sparse matrix. A must be finalized.
-SparseMatrix *Transpose(const SparseMatrix &A);
+template <class T>
+SparseMatrixMP<T> *Transpose(const SparseMatrixMP<T> &A);
 /// Transpose of a sparse matrix. A does not need to be a CSR matrix.
-SparseMatrix *TransposeAbstractSparseMatrix (const AbstractSparseMatrix &A,
-                                             int useActualWidth);
+template <class T>
+SparseMatrixMP<T> *TransposeAbstractSparseMatrix(const AbstractSparseMatrix &A,
+                                                 int useActualWidth);
 
 /// Matrix product A.B.
 /** If @a OAB is not NULL, we assume it has the structure of A.B and store the
@@ -756,78 +774,100 @@ SparseMatrix *TransposeAbstractSparseMatrix (const AbstractSparseMatrix &A,
     the result and return a pointer to it.
 
     All matrices must be finalized. */
-SparseMatrix *Mult(const SparseMatrix &A, const SparseMatrix &B,
-                   SparseMatrix *OAB = NULL);
+template <class T>
+SparseMatrixMP<T> *Mult(const SparseMatrixMP<T> &A, const SparseMatrixMP<T> &B,
+                        SparseMatrixMP<T> *OAB = NULL);
 
 /// C = A^T B
-SparseMatrix *TransposeMult(const SparseMatrix &A, const SparseMatrix &B);
+template <class T>
+SparseMatrixMP<T> *TransposeMult(const SparseMatrixMP<T> &A,
+                                 const SparseMatrixMP<T> &B);
 
 /// Matrix product of sparse matrices. A and B do not need to be CSR matrices
-SparseMatrix *MultAbstractSparseMatrix (const AbstractSparseMatrix &A,
-                                        const AbstractSparseMatrix &B);
+template <class T>
+SparseMatrixMP<T> *MultAbstractSparseMatrix(const AbstractSparseMatrix &A,
+                                            const AbstractSparseMatrix &B);
 
 /// Matrix product A.B
-DenseMatrix *Mult(const SparseMatrix &A, DenseMatrix &B);
+template <class T>
+DenseMatrix *Mult(const SparseMatrixMP<T> &A, DenseMatrix &B);
 
 /// RAP matrix product (with R=P^T)
-DenseMatrix *RAP(const SparseMatrix &A, DenseMatrix &P);
+template <class T>
+DenseMatrix *RAP(const SparseMatrixMP<T> &A, DenseMatrix &P);
 
 /// RAP matrix product (with R=P^T)
-DenseMatrix *RAP(DenseMatrix &A, const SparseMatrix &P);
+template <class T>
+DenseMatrix *RAP(DenseMatrix &A, const SparseMatrixMP<T> &P);
 
 /** RAP matrix product (with P=R^T). ORAP is like OAB above.
     All matrices must be finalized. */
-SparseMatrix *RAP(const SparseMatrix &A, const SparseMatrix &R,
-                  SparseMatrix *ORAP = NULL);
+template <class T>
+SparseMatrixMP<T> *RAP(const SparseMatrixMP<T> &A, const SparseMatrixMP<T> &R,
+                       SparseMatrixMP<T> *ORAP = NULL);
 
 /// General RAP with given R^T, A and P
-SparseMatrix *RAP(const SparseMatrix &Rt, const SparseMatrix &A,
-                  const SparseMatrix &P);
+template <class T>
+SparseMatrixMP<T> *RAP(const SparseMatrixMP<T> &Rt, const SparseMatrixMP<T> &A,
+                       const SparseMatrixMP<T> &P);
 
 /// Matrix multiplication A^t D A. All matrices must be finalized.
-SparseMatrix *Mult_AtDA(const SparseMatrix &A, const Vector &D,
-                        SparseMatrix *OAtDA = NULL);
+template <class T>
+SparseMatrixMP<T> *Mult_AtDA(const SparseMatrixMP<T> &A, const Vector &D,
+                             SparseMatrixMP<T> *OAtDA = NULL);
 
 
 /// Matrix addition result = A + B.
-SparseMatrix * Add(const SparseMatrix & A, const SparseMatrix & B);
+template <class T>
+SparseMatrixMP<T> * Add(const SparseMatrixMP<T> & A,
+                        const SparseMatrixMP<T> & B);
 /// Matrix addition result = a*A + b*B
-SparseMatrix * Add(real_t a, const SparseMatrix & A, real_t b,
-                   const SparseMatrix & B);
+template <class T, class U>
+SparseMatrixMP<T> * Add(U a, const SparseMatrixMP<T> & A, U b,
+                        const SparseMatrixMP<T> & B);
 /// Matrix addition result = sum_i A_i
-SparseMatrix * Add(Array<SparseMatrix *> & Ai);
+template <class T>
+SparseMatrixMP<T> * Add(Array<SparseMatrixMP<T> *> & Ai);
 
 /// B += alpha * A
-void Add(const SparseMatrix &A, real_t alpha, DenseMatrix &B);
+template <class T, class U>
+void Add(const SparseMatrixMP<T> &A, U alpha, DenseMatrix &B);
 
 /// Produces a block matrix with blocks A_{ij}*B
 DenseMatrix *OuterProduct(const DenseMatrix &A, const DenseMatrix &B);
 
 /// Produces a block matrix with blocks A_{ij}*B
-SparseMatrix *OuterProduct(const DenseMatrix &A, const SparseMatrix &B);
+template <class T>
+SparseMatrixMP<T> *OuterProduct(const DenseMatrix &A,
+                                const SparseMatrixMP<T> &B);
 
 /// Produces a block matrix with blocks A_{ij}*B
-SparseMatrix *OuterProduct(const SparseMatrix &A, const DenseMatrix &B);
+template <class T>
+SparseMatrixMP<T> *OuterProduct(const SparseMatrixMP<T> &A,
+                                const DenseMatrix &B);
 
 /// Produces a block matrix with blocks A_{ij}*B
-SparseMatrix *OuterProduct(const SparseMatrix &A, const SparseMatrix &B);
+template <class T>
+SparseMatrixMP<T> *OuterProduct(const SparseMatrixMP<T> &A,
+                                const SparseMatrixMP<T> &B);
 
 
 // Inline methods
 
-inline void SparseMatrix::SetColPtr(const int row) const
+template <class T>
+inline void SparseMatrixMP<T>::SetColPtr(const int row) const
 {
    if (Rows)
    {
       if (ColPtrNode == NULL)
       {
-         ColPtrNode = new RowNode *[width];
+         ColPtrNode = new RowNode<T> *[width];
          for (int i = 0; i < width; i++)
          {
             ColPtrNode[i] = NULL;
          }
       }
-      for (RowNode *node_p = Rows[row]; node_p != NULL; node_p = node_p->Prev)
+      for (RowNode<T> *node_p = Rows[row]; node_p != NULL; node_p = node_p->Prev)
       {
          ColPtrNode[node_p->Column] = node_p;
       }
@@ -850,11 +890,12 @@ inline void SparseMatrix::SetColPtr(const int row) const
    current_row = row;
 }
 
-inline void SparseMatrix::ClearColPtr() const
+template <class T>
+inline void SparseMatrixMP<T>::ClearColPtr() const
 {
    if (Rows)
    {
-      for (RowNode *node_p = Rows[current_row]; node_p != NULL;
+      for (RowNode<T> *node_p = Rows[current_row]; node_p != NULL;
            node_p = node_p->Prev)
       {
          ColPtrNode[node_p->Column] = NULL;
@@ -869,17 +910,18 @@ inline void SparseMatrix::ClearColPtr() const
    }
 }
 
-inline real_t &SparseMatrix::SearchRow(const int col)
+template <class T>
+inline T &SparseMatrixMP<T>::SearchRow(const int col)
 {
    if (Rows)
    {
-      RowNode *node_p = ColPtrNode[col];
+      RowNode<T> *node_p = ColPtrNode[col];
       if (node_p == NULL)
       {
 #ifdef MFEM_USE_MEMALLOC
          node_p = NodesMem->Alloc();
 #else
-         node_p = new RowNode;
+         node_p = new RowNode<T>;
 #endif
          node_p->Prev = Rows[current_row];
          node_p->Column = col;
@@ -896,11 +938,12 @@ inline real_t &SparseMatrix::SearchRow(const int col)
    }
 }
 
-inline real_t SparseMatrix::_Get_(const int col) const
+template <class T>
+inline T SparseMatrixMP<T>::_Get_(const int col) const
 {
    if (Rows)
    {
-      RowNode *node_p = ColPtrNode[col];
+      RowNode<T> *node_p = ColPtrNode[col];
       return (node_p == NULL) ? 0.0 : node_p->Value;
    }
    else
@@ -910,11 +953,12 @@ inline real_t SparseMatrix::_Get_(const int col) const
    }
 }
 
-inline real_t &SparseMatrix::SearchRow(const int row, const int col)
+template <class T>
+inline T &SparseMatrixMP<T>::SearchRow(const int row, const int col)
 {
    if (Rows)
    {
-      RowNode *node_p;
+      RowNode<T> *node_p;
 
       for (node_p = Rows[row]; 1; node_p = node_p->Prev)
       {
@@ -923,7 +967,7 @@ inline real_t &SparseMatrix::SearchRow(const int row, const int col)
 #ifdef MFEM_USE_MEMALLOC
             node_p = NodesMem->Alloc();
 #else
-            node_p = new RowNode;
+            node_p = new RowNode<T>;
 #endif
             node_p->Prev = Rows[row];
             node_p->Column = col;
@@ -954,7 +998,7 @@ inline real_t &SparseMatrix::SearchRow(const int row, const int col)
 }
 
 /// Specialization of the template function Swap<> for class SparseMatrix
-template<> inline void Swap<SparseMatrix>(SparseMatrix &a, SparseMatrix &b)
+template<class T> inline void Swap(SparseMatrixMP<T> &a, SparseMatrixMP<T> &b)
 {
    a.Swap(b);
 }

linalg/sparsesmoothers.cpp

@@ -20,21 +20,23 @@
 namespace mfem
 {
 
-void SparseSmoother::SetOperator(const Operator &a)
+template <class T>
+void SparseSmootherMP<T>::SetOperator(const OperatorMP<T> &a)
 {
-   oper = dynamic_cast<const SparseMatrix*>(&a);
+   oper = dynamic_cast<const SparseMatrixMP<T>*>(&a);
    if (oper == NULL)
    {
       mfem_error("SparseSmoother::SetOperator : not a SparseMatrix!");
    }
-   height = oper->Height();
-   width = oper->Width();
+   this->height = oper->Height();
+   this->width = oper->Width();
 }
 
 /// Matrix vector multiplication with GS Smoother.
-void GSSmoother::Mult(const Vector &x, Vector &y) const
+template <class T>
+void GSSmootherMP<T>::Mult(const VectorMP<T> &x, VectorMP<T> &y) const
 {
-   if (!iterative_mode)
+   if (!this->iterative_mode)
    {
       y = 0.0;
    }
@@ -42,11 +44,11 @@ void GSSmoother::Mult(const Vector &x, Vector &y) const
    {
       if (type != 2)
       {
-         oper->Gauss_Seidel_forw(x, y);
+         this->oper->Gauss_Seidel_forw(x, y);
       }
       if (type != 1)
       {
-         oper->Gauss_Seidel_back(x, y);
+         this->oper->Gauss_Seidel_back(x, y);
       }
    }
 }
@@ -108,4 +110,7 @@ void DSmoother::Mult(const Vector &x, Vector &y) const
    }
 }
 
+template class GSSmootherMP<float>;
+template class GSSmootherMP<double>;
+
 }

linalg/sparsesmoothers.hpp

@@ -18,22 +18,26 @@
 namespace mfem
 {
 
-class SparseSmoother : public MatrixInverse
+template <class T>
+class SparseSmootherMP : public MatrixInverseMP<T>
 {
 protected:
-   const SparseMatrix *oper;
+   const SparseMatrixMP<T> *oper;
 
 public:
-   SparseSmoother() { oper = NULL; }
+   SparseSmootherMP() { oper = NULL; }
 
-   SparseSmoother(const SparseMatrix &a)
-      : MatrixInverse(a) { oper = &a; }
+   SparseSmootherMP(const SparseMatrixMP<T> &a)
+      : MatrixInverseMP<T>(a) { oper = &a; }
 
-   void SetOperator(const Operator &a) override;
+   void SetOperator(const OperatorMP<T> &a) override;
 };
 
+using SparseSmoother = SparseSmootherMP<real_t>;
+
 /// Data type for Gauss-Seidel smoother of sparse matrix
-class GSSmoother : public SparseSmoother
+template <class T>
+class GSSmootherMP : public SparseSmootherMP<T>
 {
 protected:
    int type; // 0, 1, 2 - symmetric, forward, backward
@@ -41,16 +45,19 @@ protected:
 
 public:
    /// Create GSSmoother.
-   GSSmoother(int t = 0, int it = 1) { type = t; iterations = it; }
+   GSSmootherMP(int t = 0, int it = 1) { type = t; iterations = it; }
 
    /// Create GSSmoother.
-   GSSmoother(const SparseMatrix &a, int t = 0, int it = 1)
-      : SparseSmoother(a) { type = t; iterations = it; }
+   GSSmootherMP(const SparseMatrixMP<T> &a, int t = 0, int it = 1)
+      : SparseSmootherMP<T>(a) { type = t; iterations = it; }
 
    /// Matrix vector multiplication with GS Smoother.
-   void Mult(const Vector &x, Vector &y) const override;
+   void Mult(const VectorMP<T> &x, VectorMP<T> &y) const override;
 };
 
+using GSSmoother = GSSmootherMP<real_t>;
+
+
 /// Data type for scaled Jacobi-type smoother of sparse matrix
 class DSmoother : public SparseSmoother
 {

linalg/vector.hpp

@@ -40,7 +40,8 @@ namespace mfem
 
 /** Count the number of entries in an array of doubles for which isfinite
     is false, i.e. the entry is a NaN or +/-Inf. */
-inline int CheckFinite(const real_t *v, const int n);
+template <class T>
+inline int CheckFinite(const T *v, const int n);
 
 /// Define a shortcut for std::numeric_limits<double>::infinity()
 #ifndef __CYGWIN__
@@ -77,60 +78,84 @@ inline real_t rand_real()
 #endif
 }
 
+template <class T>
+class VectorMP;
+
+template <class T>
+void add(const VectorMP<T> &v1, const VectorMP<T> &v2, VectorMP<T> &v);
+
+template <class T, class U>
+void add(const VectorMP<T> &v1, U alpha, const VectorMP<T> &v2, VectorMP<T> &v);
+
+template <class T, class U>
+void add(const U a, const VectorMP<T> &x, const VectorMP<T> &y, VectorMP<T> &z);
+
+template <class T, class U>
+void add(const U a, const VectorMP<T> &x,
+         const U b, const VectorMP<T> &y, VectorMP<T> &z);
+
+template <class T>
+void subtract(const VectorMP<T> &x, const VectorMP<T> &y, VectorMP<T> &z);
+
+template <class T, class U>
+void subtract(const U a, const VectorMP<T> &x, const VectorMP<T> &y,
+              VectorMP<T> &z);
+
 /// Vector data type.
-class Vector
+template <class T>
+class VectorMP
 {
 protected:
 
-   Memory<real_t> data;
+   Memory<T> data;
    int size;
 
 public:
 
    /** Default constructor for Vector. Sets size = 0, and calls Memory::Reset on
        data through Memory<double>'s default constructor. */
-   Vector(): size(0) { }
+   VectorMP(): size(0) { }
 
    /// Copy constructor. Allocates a new data array and copies the data.
-   Vector(const Vector &);
+   VectorMP(const VectorMP<T> &);
 
    /// Move constructor. "Steals" data from its argument.
-   Vector(Vector&& v);
+   VectorMP(VectorMP<T>&& v);
 
    /// @brief Creates vector of size s.
    /// @warning Entries are not initialized to zero!
-   explicit Vector(int s);
+   explicit VectorMP(int s);
 
    /// Creates a vector referencing an array of doubles, owned by someone else.
    /** The pointer @a data_ can be NULL. The data array can be replaced later
        with SetData(). */
-   Vector(real_t *data_, int size_)
+   VectorMP(T *data_, int size_)
    { data.Wrap(data_, size_, false); size = size_; }
 
    /** @brief Create a Vector referencing a sub-vector of the Vector @a base
        starting at the given offset, @a base_offset, and size @a size_. */
-   Vector(Vector &base, int base_offset, int size_)
+   VectorMP(VectorMP<T> &base, int base_offset, int size_)
       : data(base.data, base_offset, size_), size(size_) { }
 
    /// Create a Vector of size @a size_ using MemoryType @a mt.
-   Vector(int size_, MemoryType mt)
+   VectorMP(int size_, MemoryType mt)
       : data(size_, mt), size(size_) { }
 
    /** @brief Create a Vector of size @a size_ using host MemoryType @a h_mt and
        device MemoryType @a d_mt. */
-   Vector(int size_, MemoryType h_mt, MemoryType d_mt)
+   VectorMP(int size_, MemoryType h_mt, MemoryType d_mt)
       : data(size_, h_mt, d_mt), size(size_) { }
 
    /// Create a vector from a statically sized C-style array of convertible type
    template <typename CT, int N>
-   explicit Vector(const CT (&values)[N]) : Vector(N)
+   explicit VectorMP(const CT (&values)[N]) : VectorMP(N)
    { std::copy(values, values + N, begin()); }
 
    /// Create a vector using a braced initializer list
    template <typename CT, typename std::enable_if<
-                std::is_convertible<CT,real_t>::value,bool>::type = true>
-   explicit Vector(std::initializer_list<CT> values) :
-      Vector(static_cast<int> (values.size()))
+                std::is_convertible<CT,T>::value,bool>::type = true>
+   explicit VectorMP(std::initializer_list<CT> values) :
+      VectorMP(static_cast<int> (values.size()))
    { std::copy(values.begin(), values.end(), begin()); }
 
    /// Enable execution of Vector operations using the mfem::Device.
@@ -170,7 +195,7 @@ public:
    void SetSize(int s, MemoryType mt);
 
    /// Resize the vector to size @a s using the MemoryType of @a v.
-   void SetSize(int s, const Vector &v) { SetSize(s, v.GetMemory().GetMemoryType()); }
+   void SetSize(int s, const VectorMP<T> &v) { SetSize(s, v.GetMemory().GetMemoryType()); }
 
    /// Update \ref Capacity() to @a res (if less than current), keeping existing entries.
    void Reserve(int res);
@@ -181,20 +206,20 @@ public:
 
    /// Set the Vector data.
    /// @warning This method should be called only when OwnsData() is false.
-   void SetData(real_t *d) { data.Wrap(d, data.Capacity(), false); }
+   void SetData(T *d) { data.Wrap(d, data.Capacity(), false); }
 
    /// Set the Vector data and size.
    /** The Vector does not assume ownership of the new data. The new size is
        also used as the new Capacity().
        @warning This method should be called only when OwnsData() is false.
        @sa NewDataAndSize(). */
-   void SetDataAndSize(real_t *d, int s) { data.Wrap(d, s, false); size = s; }
+   void SetDataAndSize(T *d, int s) { data.Wrap(d, s, false); size = s; }
 
    /// Set the Vector data and size, deleting the old data, if owned.
    /** The Vector does not assume ownership of the new data. The new size is
        also used as the new Capacity().
        @sa SetDataAndSize(). */
-   void NewDataAndSize(real_t *d, int s)
+   void NewDataAndSize(T *d, int s)
    {
       data.Delete();
       SetDataAndSize(d, s);
@@ -209,14 +234,14 @@ public:
        the Vector object takes ownership of all pointers owned by @a mem.
 
        @sa NewDataAndSize(). */
-   inline void NewMemoryAndSize(const Memory<real_t> &mem, int s, bool own_mem);
+   inline void NewMemoryAndSize(const Memory<T> &mem, int s, bool own_mem);
 
    /// Reset the Vector to be a reference to a sub-vector of @a base.
-   inline void MakeRef(Vector &base, int offset, int size);
+   inline void MakeRef(VectorMP<T> &base, int offset, int size);
 
    /** @brief Reset the Vector to be a reference to a sub-vector of @a base
        without changing its current size. */
-   inline void MakeRef(Vector &base, int offset);
+   inline void MakeRef(VectorMP<T> &base, int offset);
 
    /// Set the Vector data (host pointer) ownership flag.
    void MakeDataOwner() const { data.SetHostPtrOwner(true); }
@@ -240,72 +265,72 @@ public:
    /// Return a pointer to the beginning of the Vector data.
    /** @warning This method should be used with caution as it gives write access
        to the data of const-qualified Vector%s. */
-   inline real_t *GetData() const
-   { return const_cast<real_t*>((const real_t*)data); }
+   inline T *GetData() const
+   { return const_cast<T*>((const T*)data); }
 
    /// Conversion to `double *`. Deprecated.
-   MFEM_DEPRECATED inline operator real_t *() { return data; }
+   MFEM_DEPRECATED inline operator T *() { return data; }
 
    /// Conversion to `const double *`. Deprecated.
-   MFEM_DEPRECATED inline operator const real_t *() const { return data; }
+   MFEM_DEPRECATED inline operator const T *() const { return data; }
 
    /// STL-like begin.
-   inline real_t *begin() { return data; }
+   inline T *begin() { return data; }
 
    /// STL-like end.
-   inline real_t *end() { return data + size; }
+   inline T *end() { return data + size; }
 
    /// STL-like begin (const version).
-   inline const real_t *begin() const { return data; }
+   inline const T *begin() const { return data; }
 
    /// STL-like end (const version).
-   inline const real_t *end() const { return data + size; }
+   inline const T *end() const { return data + size; }
 
    /// Return a reference to the Memory object used by the Vector.
-   Memory<real_t> &GetMemory() { return data; }
+   Memory<T> &GetMemory() { return data; }
 
    /** @brief Return a reference to the Memory object used by the Vector, const
        version. */
-   const Memory<real_t> &GetMemory() const { return data; }
+   const Memory<T> &GetMemory() const { return data; }
 
    /// Update the memory location of the vector to match @a v.
-   void SyncMemory(const Vector &v) const { GetMemory().Sync(v.GetMemory()); }
+   void SyncMemory(const VectorMP<T> &v) const { GetMemory().Sync(v.GetMemory()); }
 
    /// Update the alias memory location of the vector to match @a v.
-   void SyncAliasMemory(const Vector &v) const
+   void SyncAliasMemory(const VectorMP<T> &v) const
    { GetMemory().SyncAlias(v.GetMemory(),Size()); }
 
    /// Read the Vector data (host pointer) ownership flag.
    inline bool OwnsData() const { return data.OwnsHostPtr(); }
 
    /// Changes the ownership of the data; after the call the Vector is empty
-   inline void StealData(real_t **p)
+   inline void StealData(T **p)
    { *p = data; data.Reset(); size = 0; }
 
    /// Changes the ownership of the data; after the call the Vector is empty
-   inline real_t *StealData() { real_t *p; StealData(&p); return p; }
+   inline T *StealData() { T *p; StealData(&p); return p; }
 
    /// Access Vector entries. Index i = 0 .. size-1.
-   real_t &Elem(int i);
+   T &Elem(int i);
 
    /// Read only access to Vector entries. Index i = 0 .. size-1.
-   const real_t &Elem(int i) const;
+   const T &Elem(int i) const;
 
    /// Access Vector entries using () for 0-based indexing.
    /** @note If MFEM_DEBUG is enabled, bounds checking is performed. */
-   inline real_t &operator()(int i);
+   inline T &operator()(int i);
 
    /// Read only access to Vector entries using () for 0-based indexing.
    /** @note If MFEM_DEBUG is enabled, bounds checking is performed. */
-   inline const real_t &operator()(int i) const;
+   inline const T &operator()(int i) const;
 
    /// Access Vector entries using [] for 0-based indexing.
    /** @note If MFEM_DEBUG is enabled, bounds checking is performed. */
-   inline real_t &operator[](int i) { return (*this)(i); }
+   inline T &operator[](int i) { return (*this)(i); }
 
    /// Read only access to Vector entries using [] for 0-based indexing.
    /** @note If MFEM_DEBUG is enabled, bounds checking is performed. */
-   inline const real_t &operator[](int i) const { return (*this)(i); }
+   inline const T &operator[](int i) const { return (*this)(i); }
 
    /// Dot product with a `double *` array.
    /// This function always executes on the CPU. A HostRead() will be called if
@@ -313,54 +338,54 @@ public:
    /// To optionally execute on the device:
    /// Vector tmp(v, Size());
    /// res = (*this) * tmp;
-   real_t operator*(const real_t *v) const;
+   T operator*(const T *v) const;
 
    /// Return the inner-product.
-   real_t operator*(const Vector &v) const;
+   T operator*(const VectorMP<T> &v) const;
 
    /// Copy Size() entries from @a v.
-   Vector &operator=(const real_t *v);
+   VectorMP &operator=(const T *v);
 
    /// Copy assignment.
    /** @note Defining this method overwrites the implicitly defined copy
        assignment operator. */
-   Vector &operator=(const Vector &v);
+   VectorMP &operator=(const VectorMP<T> &v);
 
    /// Move assignment
-   Vector &operator=(Vector&& v);
+   VectorMP &operator=(VectorMP<T>&& v);
 
    /// Redefine '=' for vector = constant.
-   Vector &operator=(real_t value);
+   VectorMP &operator=(T value);
 
-   Vector &operator*=(real_t c);
+   VectorMP &operator*=(T c);
 
    /// Component-wise scaling: (*this)(i) *= v(i)
-   Vector &operator*=(const Vector &v);
+   VectorMP &operator*=(const VectorMP<T> &v);
 
-   Vector &operator/=(real_t c);
+   VectorMP &operator/=(T c);
 
    /// Component-wise division: (*this)(i) /= v(i)
-   Vector &operator/=(const Vector &v);
+   VectorMP &operator/=(const VectorMP<T> &v);
 
-   Vector &operator-=(real_t c);
+   VectorMP &operator-=(T c);
 
-   Vector &operator-=(const Vector &v);
+   VectorMP &operator-=(const VectorMP<T> &v);
 
-   Vector &operator+=(real_t c);
+   VectorMP &operator+=(T c);
 
-   Vector &operator+=(const Vector &v);
+   VectorMP &operator+=(const VectorMP<T> &v);
 
    /// (*this) += a * Va
-   Vector &Add(const real_t a, const Vector &Va);
+   VectorMP &Add(const T a, const VectorMP<T> &Va);
 
    /// (*this) = a * x
-   Vector &Set(const real_t a, const Vector &x);
+   VectorMP &Set(const T a, const VectorMP<T> &x);
 
    /// (*this)[i + offset] = v[i]
-   void SetVector(const Vector &v, int offset);
+   void SetVector(const VectorMP<T> &v, int offset);
 
    /// (*this)[i + offset] += v[i]
-   void AddSubVector(const Vector &v, int offset);
+   void AddSubVector(const VectorMP<T> &v, int offset);
 
    /// (*this) = -(*this)
    void Neg();
@@ -372,53 +397,57 @@ public:
    void Abs();
 
    /// (*this)(i) = pow((*this)(i), p)
-   void Pow(const real_t p);
+   void Pow(const T p);
 
    /// Swap the contents of two Vectors
    /** Implemented without using move assignment, avoiding Destroy() calls. */
-   inline void Swap(Vector &other);
+   inline void Swap(VectorMP<T> &other);
 
    /// Set v = v1 + v2.
-   friend void add(const Vector &v1, const Vector &v2, Vector &v);
+   friend void add<T>(const VectorMP<T> &v1, const VectorMP<T> &v2,
+                      VectorMP<T> &v);
 
    /// Set v = v1 + alpha * v2.
-   friend void add(const Vector &v1, real_t alpha, const Vector &v2, Vector &v);
+   friend void add<T>(const VectorMP<T> &v1, T alpha, const VectorMP<T> &v2,
+                      VectorMP<T> &v);
 
    /// z = a * (x + y)
-   friend void add(const real_t a, const Vector &x, const Vector &y, Vector &z);
+   friend void add<T>(const T a, const VectorMP<T> &x, const VectorMP<T> &y,
+                      VectorMP<T> &z);
 
    /// z = a * x + b * y
-   friend void add(const real_t a, const Vector &x,
-                   const real_t b, const Vector &y, Vector &z);
+   friend void add<T>(const T a, const VectorMP<T> &x,
+                      const T b, const VectorMP<T> &y, VectorMP<T> &z);
 
    /// Set v = v1 - v2.
-   friend void subtract(const Vector &v1, const Vector &v2, Vector &v);
+   friend void subtract<T>(const VectorMP<T> &v1, const VectorMP<T> &v2,
+                           VectorMP<T> &v);
 
    /// z = a * (x - y)
-   friend void subtract(const real_t a, const Vector &x,
-                        const Vector &y, Vector &z);
+   friend void subtract<T>(const T a, const VectorMP<T> &x,
+                           const VectorMP<T> &y, VectorMP<T> &z);
 
    /// Computes cross product of this vector with another 3D vector.
    /// vout = this x vin.
-   void cross3D(const Vector &vin, Vector &vout) const;
+   void cross3D(const VectorMP<T> &vin, VectorMP<T> &vout) const;
 
    /// v = median(v,lo,hi) entrywise.  Implementation assumes lo <= hi.
-   void median(const Vector &lo, const Vector &hi);
+   void median(const VectorMP<T> &lo, const VectorMP<T> &hi);
 
    /// Extract entries listed in @a dofs to the output Vector @a elemvect.
    /** Negative dof values cause the -dof-1 position in @a elemvect to receive
        the -val in from this Vector. */
-   void GetSubVector(const Array<int> &dofs, Vector &elemvect) const;
+   void GetSubVector(const Array<int> &dofs, VectorMP<T> &elemvect) const;
 
    /// Extract entries listed in @a dofs to the output array @a elem_data.
    /** Negative dof values cause the -dof-1 position in @a elem_data to receive
        the -val in from this Vector. */
-   void GetSubVector(const Array<int> &dofs, real_t *elem_data) const;
+   void GetSubVector(const Array<int> &dofs, T *elem_data) const;
 
    /// Set the entries listed in @a dofs to the given @a value.
    /** Negative dof values cause the -dof-1 position in this Vector to receive
        the -value. */
-   void SetSubVector(const Array<int> &dofs, const real_t value);
+   void SetSubVector(const Array<int> &dofs, const T value);
 
    /// Set the entries listed in @a dofs to the given @a value (always on host).
    /** Negative dof values cause the -dof-1 position in this Vector to receive
@@ -427,36 +456,36 @@ public:
        As opposed to SetSubVector(const Array<int>&, const real_t), this
        function will execute only on host, even if the vector or the @a dofs
        array have the device flag set. */
-   void SetSubVectorHost(const Array<int> &dofs, const real_t value);
+   void SetSubVectorHost(const Array<int> &dofs, const T value);
 
    /** @brief Set the entries listed in @a dofs to the values given in the @a
        elemvect Vector. Negative dof values cause the -dof-1 position in this
        Vector to receive the -val from @a elemvect. */
-   void SetSubVector(const Array<int> &dofs, const Vector &elemvect);
+   void SetSubVector(const Array<int> &dofs, const VectorMP<T> &elemvect);
 
    /** @brief Set the entries listed in @a dofs to the values given the @a ,
        elem_data array. Negative dof values cause the -dof-1 position in this
        Vector to receive the -val from @a elem_data. */
-   void SetSubVector(const Array<int> &dofs, real_t *elem_data);
+   void SetSubVector(const Array<int> &dofs, T *elem_data);
 
    /** @brief Add elements of the @a elemvect Vector to the entries listed in @a
        dofs. Negative dof values cause the -dof-1 position in this Vector to add
        the -val from @a elemvect. */
-   void AddElementVector(const Array<int> & dofs, const Vector & elemvect);
+   void AddElementVector(const Array<int> & dofs, const VectorMP<T> & elemvect);
 
    /** @brief Add elements of the @a elem_data array to the entries listed in @a
        dofs. Negative dof values cause the -dof-1 position in this Vector to add
        the -val from @a elem_data. */
-   void AddElementVector(const Array<int> & dofs, real_t *elem_data);
+   void AddElementVector(const Array<int> & dofs, T *elem_data);
 
    /** @brief Add @a times the elements of the @a elemvect Vector to the entries
        listed in @a dofs. Negative dof values cause the -dof-1 position in this
        Vector to add the -a*val from @a elemvect. */
-   void AddElementVector(const Array<int> & dofs, const real_t a,
-                         const Vector & elemvect);
+   void AddElementVector(const Array<int> & dofs, const T a,
+                         const VectorMP<T> & elemvect);
 
    /// Set all vector entries NOT in the @a dofs Array to the given @a val.
-   void SetSubVectorComplement(const Array<int> &dofs, const real_t val);
+   void SetSubVectorComplement(const Array<int> &dofs, const T val);
 
    /// Prints vector to stream out.
    void Print(std::ostream &out = mfem::out, int width = 8) const;
@@ -487,61 +516,63 @@ public:
    /// Set random values in the vector.
    void Randomize(int seed = 0);
    /// Returns the l2 norm of the vector.
-   real_t Norml2() const;
+   T Norml2() const;
    /// Returns the l_infinity norm of the vector.
-   real_t Normlinf() const;
+   T Normlinf() const;
    /// Returns the l_1 norm of the vector.
-   real_t Norml1() const;
+   T Norml1() const;
    /// Returns the l_p norm of the vector.
-   real_t Normlp(real_t p) const;
+   T Normlp(T p) const;
    /// Returns the maximal element of the vector.
-   real_t Max() const;
+   T Max() const;
    /// Returns the minimal element of the vector.
-   real_t Min() const;
+   T Min() const;
    /// Return the sum of the vector entries
-   real_t Sum() const;
+   T Sum() const;
    /// Compute the square of the Euclidean distance to another vector.
-   inline real_t DistanceSquaredTo(const real_t *p) const;
+   inline T DistanceSquaredTo(const T *p) const;
    /// Compute the square of the Euclidean distance to another vector.
-   inline real_t DistanceSquaredTo(const Vector &p) const;
+   inline T DistanceSquaredTo(const VectorMP<T> &p) const;
    /// Compute the Euclidean distance to another vector.
-   inline real_t DistanceTo(const real_t *p) const;
+   inline T DistanceTo(const T *p) const;
    /// Compute the Euclidean distance to another vector.
-   inline real_t DistanceTo(const Vector &p) const;
+   inline T DistanceTo(const VectorMP<T> &p) const;
 
    /** @brief Count the number of entries in the Vector for which isfinite
        is false, i.e. the entry is a NaN or +/-Inf. */
    int CheckFinite() const { return mfem::CheckFinite(HostRead(), size); }
 
    /// Destroys vector.
-   virtual ~Vector();
+   virtual ~VectorMP<T>();
 
    /// Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
-   virtual const real_t *Read(bool on_dev = true) const
+   virtual const T *Read(bool on_dev = true) const
    { return mfem::Read(data, size, on_dev); }
 
    /// Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), false).
-   virtual const real_t *HostRead() const
+   virtual const T *HostRead() const
    { return mfem::Read(data, size, false); }
 
    /// Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
-   virtual real_t *Write(bool on_dev = true)
+   virtual T *Write(bool on_dev = true)
    { return mfem::Write(data, size, on_dev); }
 
    /// Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), false).
-   virtual real_t *HostWrite()
+   virtual T *HostWrite()
    { return mfem::Write(data, size, false); }
 
    /// Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), on_dev).
-   virtual real_t *ReadWrite(bool on_dev = true)
+   virtual T *ReadWrite(bool on_dev = true)
    { return mfem::ReadWrite(data, size, on_dev); }
 
    /// Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), false).
-   virtual real_t *HostReadWrite()
+   virtual T *HostReadWrite()
    { return mfem::ReadWrite(data, size, false); }
 
 };
 
+using Vector = VectorMP<real_t>;
+
 // Inline methods
 
 template <typename T>
@@ -550,7 +581,8 @@ inline T ZeroSubnormal(T val)
    return (std::fpclassify(val) == FP_SUBNORMAL) ? 0.0 : val;
 }
 
-inline bool IsFinite(const real_t &val)
+template <class T>
+inline bool IsFinite(const T &val)
 {
    // isfinite didn't appear in a standard until C99, and later C++11. It wasn't
    // standard in C89 or C++98. PGI as of 14.7 still defines it as a macro.
@@ -561,7 +593,8 @@ inline bool IsFinite(const real_t &val)
 #endif
 }
 
-inline int CheckFinite(const real_t *v, const int n)
+template <class T>
+inline int CheckFinite(const T *v, const int n)
 {
    int bad = 0;
    for (int i = 0; i < n; i++)
@@ -571,7 +604,8 @@ inline int CheckFinite(const real_t *v, const int n)
    return bad;
 }
 
-inline Vector::Vector(int s)
+template <class T>
+inline VectorMP<T>::VectorMP(int s)
 {
    MFEM_ASSERT(s>=0,"Unexpected negative size.");
    size = s;
@@ -581,7 +615,8 @@ inline Vector::Vector(int s)
    }
 }
 
-inline void Vector::SetSize(int s)
+template <class T>
+inline void VectorMP<T>::SetSize(int s)
 {
    if (s == size)
    {
@@ -601,7 +636,8 @@ inline void Vector::SetSize(int s)
    data.UseDevice(use_dev);
 }
 
-inline void Vector::SetSize(int s, MemoryType mt)
+template <class T>
+inline void VectorMP<T>::SetSize(int s, MemoryType mt)
 {
    if (mt == data.GetMemoryType())
    {
@@ -630,11 +666,12 @@ inline void Vector::SetSize(int s, MemoryType mt)
    data.UseDevice(use_dev);
 }
 
-inline void Vector::Reserve(int res)
+template <class T>
+inline void VectorMP<T>::Reserve(int res)
 {
    if (res > Capacity())
    {
-      Memory<real_t> p(res, data.GetMemoryType());
+      Memory<T> p(res, data.GetMemoryType());
       p.CopyFrom(data, size);
       p.UseDevice(data.UseDevice());
       data.Delete();
@@ -642,8 +679,9 @@ inline void Vector::Reserve(int res)
    }
 }
 
-inline void Vector::NewMemoryAndSize(const Memory<real_t> &mem, int s,
-                                     bool own_mem)
+template <class T>
+inline void VectorMP<T>::NewMemoryAndSize(const Memory<T> &mem, int s,
+                                          bool own_mem)
 {
    data.Delete();
    size = s;
@@ -657,20 +695,23 @@ inline void Vector::NewMemoryAndSize(const Memory<real_t> &mem, int s,
    }
 }
 
-inline void Vector::MakeRef(Vector &base, int offset, int s)
+template <class T>
+inline void VectorMP<T>::MakeRef(VectorMP<T> &base, int offset, int s)
 {
    data.Delete();
    size = s;
    data.MakeAlias(base.GetMemory(), offset, s);
 }
 
-inline void Vector::MakeRef(Vector &base, int offset)
+template <class T>
+inline void VectorMP<T>::MakeRef(VectorMP<T> &base, int offset)
 {
    data.Delete();
    data.MakeAlias(base.GetMemory(), offset, size);
 }
 
-inline void Vector::Destroy()
+template <class T>
+inline void VectorMP<T>::Destroy()
 {
    const bool use_dev = data.UseDevice();
    data.Delete();  // calls data.Reset(h_mt) as well
@@ -678,7 +719,8 @@ inline void Vector::Destroy()
    data.UseDevice(use_dev);
 }
 
-inline real_t &Vector::operator()(int i)
+template <class T>
+inline T &VectorMP<T>::operator()(int i)
 {
    MFEM_ASSERT(data && i >= 0 && i < size,
                "index [" << i << "] is out of range [0," << size << ")");
@@ -686,7 +728,8 @@ inline real_t &Vector::operator()(int i)
    return data[i];
 }
 
-inline const real_t &Vector::operator()(int i) const
+template <class T>
+inline const T &VectorMP<T>::operator()(int i) const
 {
    MFEM_ASSERT(data && i >= 0 && i < size,
                "index [" << i << "] is out of range [0," << size << ")");
@@ -694,7 +737,8 @@ inline const real_t &Vector::operator()(int i) const
    return data[i];
 }
 
-inline void Vector::Swap(Vector &other)
+template <class T>
+inline void VectorMP<T>::Swap(VectorMP<T> &other)
 {
    mfem::Swap(data, other.data);
    mfem::Swap(size, other.size);
@@ -702,19 +746,22 @@ inline void Vector::Swap(Vector &other)
 
 /** @brief Swap of Vector objects for use with standard library algorithms.
     Also, used by mfem::Swap(). */
-inline void swap(Vector &a, Vector &b)
+template <class T>
+inline void swap(VectorMP<T> &a, VectorMP<T> &b)
 {
    a.Swap(b);
 }
 
-inline Vector::~Vector()
+template <class T>
+inline VectorMP<T>::~VectorMP()
 {
    data.Delete();
 }
 
-inline real_t DistanceSquared(const real_t *x, const real_t *y, const int n)
+template <class T>
+inline T DistanceSquared(const T *x, const T *y, const int n)
 {
-   real_t d = 0.0;
+   T d = 0.0;
 
    for (int i = 0; i < n; i++)
    {
@@ -724,43 +771,50 @@ inline real_t DistanceSquared(const real_t *x, const real_t *y, const int n)
    return d;
 }
 
-inline real_t Distance(const real_t *x, const real_t *y, const int n)
+template <class T>
+inline T Distance(const T *x, const T *y, const int n)
 {
    return std::sqrt(DistanceSquared(x, y, n));
 }
 
-inline real_t Distance(const Vector &x, const Vector &y)
+template <class T>
+inline T Distance(const VectorMP<T> &x, const VectorMP<T> &y)
 {
    return x.DistanceTo(y);
 }
 
-inline real_t Vector::DistanceSquaredTo(const real_t *p) const
+template <class T>
+inline T VectorMP<T>::DistanceSquaredTo(const T *p) const
 {
-   return DistanceSquared(data, p, size);
+   return DistanceSquared<T>(data, p, size);
 }
 
-inline real_t Vector::DistanceSquaredTo(const Vector &p) const
+template <class T>
+inline T VectorMP<T>::DistanceSquaredTo(const VectorMP<T> &p) const
 {
    MFEM_ASSERT(p.Size() == Size(), "Incompatible vector sizes.");
-   return DistanceSquared(data, p.data, size);
+   return DistanceSquared<T>(data, p.data, size);
 }
 
-inline real_t Vector::DistanceTo(const real_t *p) const
+template <class T>
+inline T VectorMP<T>::DistanceTo(const T *p) const
 {
-   return Distance(data, p, size);
+   return Distance<T>(data, p, size);
 }
 
-inline real_t Vector::DistanceTo(const Vector &p) const
+template <class T>
+inline T VectorMP<T>::DistanceTo(const VectorMP<T> &p) const
 {
    MFEM_ASSERT(p.Size() == Size(), "Incompatible vector sizes.");
-   return Distance(data, p.data, size);
+   return Distance<T>(data, p.data, size);
 }
 
 /// Returns the inner product of x and y
 /** In parallel this computes the inner product of the local vectors,
     producing different results on each MPI rank.
 */
-inline real_t InnerProduct(const Vector &x, const Vector &y)
+template <class T>
+inline T InnerProduct(const VectorMP<T> &x, const VectorMP<T> &y)
 {
    return x * y;
 }
@@ -770,11 +824,25 @@ inline real_t InnerProduct(const Vector &x, const Vector &y)
 /** In parallel this computes the inner product of the global vectors,
     producing identical results on each MPI rank.
 */
-inline real_t InnerProduct(MPI_Comm comm, const Vector &x, const Vector &y)
+template <class T>
+inline T InnerProduct(MPI_Comm comm, const VectorMP<T> &x, const VectorMP<T> &y)
 {
-   real_t loc_prod = x * y;
-   real_t glb_prod;
-   MPI_Allreduce(&loc_prod, &glb_prod, 1, MFEM_MPI_REAL_T, MPI_SUM, comm);
+   T loc_prod = x * y;
+   T glb_prod;
+
+   if (std::is_same<T, double>::value)
+   {
+      MPI_Allreduce(&loc_prod, &glb_prod, 1, MPI_DOUBLE, MPI_SUM, comm);
+   }
+   else if (std::is_same<T, float>::value)
+   {
+      MPI_Allreduce(&loc_prod, &glb_prod, 1, MPI_FLOAT, MPI_SUM, comm);
+   }
+   else
+   {
+      MFEM_ABORT("Floating point type not supported");
+   }
+
    return glb_prod;
 }
 #endif

miniapps/meshing/minimal-surface.cpp

@@ -389,7 +389,7 @@ public:
             add(*nodes, delta, *nodes);
          }
          // x = lambda*nodes + (1-lambda)*x
-         add(lambda, *nodes, (1.0-lambda), x, x);
+         add(lambda, *nodes, (real_t)(1.0-lambda), x, x);
          return Converged(rnorm);
       }
 

miniapps/meshing/pminimal-surface.cpp

@@ -396,7 +396,7 @@ public:
             add(*nodes, delta, *nodes);
          }
          // x = lambda*nodes + (1-lambda)*x
-         add(lambda, *nodes, (1.0-lambda), x, x);
+         add(lambda, *nodes, (real_t)(1.0-lambda), x, x);
          return Converged(rnorm);
       }