add optional MemoryType to the SparseMatrix copy-ctor

There are things that hurt performance in laster mfem@master commits,
not sure what yet.

This includes:
- feature/artv3/cusparse-Spmv
- feature/tomstitt/temp-mem-type
- patches to tmop.cpp and tmop_tools.cpp to address memory issues
- patch to mem_manager/device to use mfem's default allcator instead of
the host umpire one

config/defaults.mk

@@ -125,6 +125,7 @@ MFEM_USE_GINKGO        = NO
 MFEM_USE_GNUTLS        = NO
 MFEM_USE_NETCDF        = NO
 MFEM_USE_PETSC         = NO
+MFEM_USE_SLEPC         = NO
 MFEM_USE_MPFR          = NO
 MFEM_USE_SIDRE         = NO
 MFEM_USE_CONDUIT       = NO
@@ -276,6 +277,20 @@ ifeq ($(PETSC_FOUND),YES)
       -L$(abspath $(PETSC_DIR))/lib -lpetsc $(PETSC_LIB)
 endif
 
+SLEPC_DIR := $(MFEM_DIR)/../slepc
+SLEPC_VARS := $(SLEPC_DIR)/lib/slepc/conf/slepc_variables
+SLEPC_FOUND := $(if $(wildcard $(SLEPC_VARS)),YES,)
+SLEPC_INC_VAR = SLEPC_INCLUDE
+SLEPC_LIB_VAR = SLEPC_EXTERNAL_LIB
+ifeq ($(SLEPC_FOUND),YES)
+   SLEPC_OPT := $(shell sed -n "s/$(SLEPC_INC_VAR) *= *//p" $(SLEPC_VARS))
+   # Some additional external libraries might be defined in this file
+   -include ${SLEPC_DIR}/${PETSC_ARCH}/lib/slepc/conf/slepcvariables
+   SLEPC_LIB := $(shell sed -n "s/$(SLEPC_LIB_VAR) *= *//p" $(SLEPC_VARS))
+   SLEPC_LIB := -Wl,-rpath,$(abspath $(SLEPC_DIR))/$(PETSC_ARCH)/lib\
+      -L$(abspath $(SLEPC_DIR))/$(PETSC_ARCH)/lib -lslepc $(SLEPC_LIB)
+endif
+
 # MPFR library configuration
 MPFR_OPT =
 MPFR_LIB = -lmpfr
@@ -324,9 +339,9 @@ GSLIB_DIR = @MFEM_DIR@/../gslib/build
 GSLIB_OPT = -I$(GSLIB_DIR)/include
 GSLIB_LIB = -L$(GSLIB_DIR)/lib -lgs
 
-# CUDA library configuration (currently not needed)
+# CUDA library configuration
 CUDA_OPT =
-CUDA_LIB =
+CUDA_LIB = -lcusparse
 
 # HIP library configuration (currently not needed)
 HIP_OPT =

fem/bilininteg_diffusion_pa.cpp

@@ -176,42 +176,41 @@ static void PADiffusionSetup3D(const int Q1D,
    auto J = Reshape(j.Read(), NQ, 3, 3, NE);
    auto C = const_c ? Reshape(c.Read(), 1, 1) : Reshape(c.Read(), NQ, NE);
    auto D = Reshape(d.Write(), NQ, 6, NE);
-   MFEM_FORALL(e, NE,
+   MFEM_FORALL(eq, NE*NQ,
    {
-      for (int q = 0; q < NQ; ++q)
-      {
-         const double J11 = J(q,0,0,e);
-         const double J21 = J(q,1,0,e);
-         const double J31 = J(q,2,0,e);
-         const double J12 = J(q,0,1,e);
-         const double J22 = J(q,1,1,e);
-         const double J32 = J(q,2,1,e);
-         const double J13 = J(q,0,2,e);
-         const double J23 = J(q,1,2,e);
-         const double J33 = J(q,2,2,e);
-         const double detJ = J11 * (J22 * J33 - J32 * J23) -
-         /* */               J21 * (J12 * J33 - J32 * J13) +
-         /* */               J31 * (J12 * J23 - J22 * J13);
-         const double coeff = const_c ? C(0,0) : C(q,e);
-         const double c_detJ = W[q] * coeff / detJ;
-         // adj(J)
-         const double A11 = (J22 * J33) - (J23 * J32);
-         const double A12 = (J32 * J13) - (J12 * J33);
-         const double A13 = (J12 * J23) - (J22 * J13);
-         const double A21 = (J31 * J23) - (J21 * J33);
-         const double A22 = (J11 * J33) - (J13 * J31);
-         const double A23 = (J21 * J13) - (J11 * J23);
-         const double A31 = (J21 * J32) - (J31 * J22);
-         const double A32 = (J31 * J12) - (J11 * J32);
-         const double A33 = (J11 * J22) - (J12 * J21);
-         // detJ J^{-1} J^{-T} = (1/detJ) adj(J) adj(J)^T
-         D(q,0,e) = c_detJ * (A11*A11 + A12*A12 + A13*A13); // 1,1
-         D(q,1,e) = c_detJ * (A11*A21 + A12*A22 + A13*A23); // 2,1
-         D(q,2,e) = c_detJ * (A11*A31 + A12*A32 + A13*A33); // 3,1
-         D(q,3,e) = c_detJ * (A21*A21 + A22*A22 + A23*A23); // 2,2
-         D(q,4,e) = c_detJ * (A21*A31 + A22*A32 + A23*A33); // 3,2
-         D(q,5,e) = c_detJ * (A31*A31 + A32*A32 + A33*A33); // 3,3
-      }
+      const int e = eq / NQ;
+      const int q = eq % NQ;
+      const double J11 = J(q,0,0,e);
+      const double J21 = J(q,1,0,e);
+      const double J31 = J(q,2,0,e);
+      const double J12 = J(q,0,1,e);
+      const double J22 = J(q,1,1,e);
+      const double J32 = J(q,2,1,e);
+      const double J13 = J(q,0,2,e);
+      const double J23 = J(q,1,2,e);
+      const double J33 = J(q,2,2,e);
+      const double detJ = J11 * (J22 * J33 - J32 * J23) -
+      /* */               J21 * (J12 * J33 - J32 * J13) +
+      /* */               J31 * (J12 * J23 - J22 * J13);
+      const double coeff = const_c ? C(0,0) : C(q,e);
+      const double c_detJ = W[q] * coeff / detJ;
+      // adj(J)
+      const double A11 = (J22 * J33) - (J23 * J32);
+      const double A12 = (J32 * J13) - (J12 * J33);
+      const double A13 = (J12 * J23) - (J22 * J13);
+      const double A21 = (J31 * J23) - (J21 * J33);
+      const double A22 = (J11 * J33) - (J13 * J31);
+      const double A23 = (J21 * J13) - (J11 * J23);
+      const double A31 = (J21 * J32) - (J31 * J22);
+      const double A32 = (J31 * J12) - (J11 * J32);
+      const double A33 = (J11 * J22) - (J12 * J21);
+      // detJ J^{-1} J^{-T} = (1/detJ) adj(J) adj(J)^T
+      D(q,0,e) = c_detJ * (A11*A11 + A12*A12 + A13*A13); // 1,1
+      D(q,1,e) = c_detJ * (A11*A21 + A12*A22 + A13*A23); // 2,1
+      D(q,2,e) = c_detJ * (A11*A31 + A12*A32 + A13*A33); // 3,1
+      D(q,3,e) = c_detJ * (A21*A21 + A22*A22 + A23*A23); // 2,2
+      D(q,4,e) = c_detJ * (A21*A31 + A22*A32 + A23*A33); // 3,2
+      D(q,5,e) = c_detJ * (A31*A31 + A32*A32 + A33*A33); // 3,3
    });
 }
 

fem/bilininteg_mass_pa.cpp

@@ -25,6 +25,7 @@ namespace mfem
 
 void MassIntegrator::SetupPA(const FiniteElementSpace &fes, const bool force)
 {
+
    // Assuming the same element type
    fespace = &fes;
    Mesh *mesh = fes.GetMesh();
@@ -51,21 +52,30 @@ void MassIntegrator::SetupPA(const FiniteElementSpace &fes, const bool force)
    dofs1D = maps->ndof;
    quad1D = maps->nqpt;
    pa_data.SetSize(ne*nq, Device::GetDeviceMemoryType());
-   Vector coeff;
+   Vector *coeff{nullptr};
+   bool own_coeff{true};
    if (Q == nullptr)
    {
-      coeff.SetSize(1);
-      coeff(0) = 1.0;
+      coeff = new Vector;
+      coeff->SetSize(1);
+      (*coeff)(0) = 1.0;
    }
    else if (ConstantCoefficient* cQ = dynamic_cast<ConstantCoefficient*>(Q))
    {
-      coeff.SetSize(1);
-      coeff(0) = cQ->constant;
+      coeff = new Vector;
+      coeff->SetSize(1);
+      (*coeff)(0) = 1.0;
+   }
+   else if (QuadratureCoefficient* cQ = dynamic_cast<QuadratureCoefficient*>(Q))
+   {
+      coeff = cQ->Data();
+      own_coeff = false;
    }
    else
    {
-      coeff.SetSize(nq * ne);
-      auto C = Reshape(coeff.HostWrite(), nq, ne);
+      coeff = new Vector;
+      coeff->SetSize(nq * ne);
+      auto C = Reshape(coeff->HostWrite(), nq, ne);
       for (int e = 0; e < ne; ++e)
       {
          ElementTransformation& T = *fes.GetElementTransformation(e);
@@ -80,11 +90,11 @@ void MassIntegrator::SetupPA(const FiniteElementSpace &fes, const bool force)
    {
       const int NE = ne;
       const int NQ = nq;
-      const bool const_c = coeff.Size() == 1;
+      const bool const_c = coeff->Size() == 1;
       auto w = ir->GetWeights().Read();
       auto J = Reshape(geom->J.Read(), NQ,2,2,NE);
       auto C =
-         const_c ? Reshape(coeff.Read(), 1,1) : Reshape(coeff.Read(), NQ,NE);
+         const_c ? Reshape(coeff->Read(), 1,1) : Reshape(coeff->Read(), NQ,NE);
       auto v = Reshape(pa_data.Write(), NQ, NE);
       MFEM_FORALL(e, NE,
       {
@@ -104,11 +114,11 @@ void MassIntegrator::SetupPA(const FiniteElementSpace &fes, const bool force)
    {
       const int NE = ne;
       const int NQ = nq;
-      const bool const_c = coeff.Size() == 1;
+      const bool const_c = coeff->Size() == 1;
       auto W = ir->GetWeights().Read();
       auto J = Reshape(geom->J.Read(), NQ,3,3,NE);
       auto C =
-         const_c ? Reshape(coeff.Read(), 1,1) : Reshape(coeff.Read(), NQ,NE);
+         const_c ? Reshape(coeff->Read(), 1,1) : Reshape(coeff->Read(), NQ,NE);
       auto v = Reshape(pa_data.Write(), NQ,NE);
       MFEM_FORALL(e, NE,
       {
@@ -125,6 +135,8 @@ void MassIntegrator::SetupPA(const FiniteElementSpace &fes, const bool force)
          }
       });
    }
+
+   if (own_coeff) { delete coeff; }
 }
 
 void MassIntegrator::AssemblePA(const FiniteElementSpace &fes)

fem/coefficient.cpp

@@ -12,6 +12,7 @@
 // Implementation of Coefficient class
 
 #include "fem.hpp"
+#include "../linalg/dtensor.hpp"
 
 #include <cmath>
 #include <limits>
@@ -21,6 +22,13 @@ namespace mfem
 
 using namespace std;
 
+double QuadratureCoefficient::Eval(ElementTransformation & T,
+                                   const IntegrationPoint & ip)
+{
+   auto coeff = mfem::Reshape(qData->HostRead(), nip, NE);
+   return coeff(ip.index, T.ElementNo);
+}
+
 double PWConstCoefficient::Eval(ElementTransformation & T,
                                 const IntegrationPoint & ip)
 {

fem/coefficient.hpp

@@ -84,6 +84,33 @@ public:
    { return (constant); }
 };
 
+
+/// class for quadrature coefficient
+class QuadratureCoefficient : public Coefficient
+{
+
+private:
+   const int nip;
+   const int NE;
+public:
+   Vector *qData{nullptr};
+
+   //Set external data
+   QuadratureCoefficient(Vector *Data, int in_nip, int in_NE)
+      : qData(Data), nip(in_nip), NE(in_NE)
+   { }
+
+   virtual double Eval(ElementTransformation &T,
+                       const IntegrationPoint &ip);
+
+   Vector *Data()
+   {
+      return qData;
+   }
+
+};
+
+/// class for piecewise constant coefficient
 /** @brief A piecewise constant coefficient with the constants keyed
     off the element attribute numbers. */
 class PWConstCoefficient : public Coefficient

fem/pgridfunc.cpp

@@ -214,7 +214,7 @@ void ParGridFunction::ExchangeFaceNbrData()
    ParMesh *pmesh = pfes->GetParMesh();
 
    face_nbr_data.SetSize(pfes->GetFaceNbrVSize());
-   Vector send_data(pfes->send_face_nbr_ldof.Size_of_connections());
+   send_data.SetSize(pfes->send_face_nbr_ldof.Size_of_connections());
 
    int *send_offset = pfes->send_face_nbr_ldof.GetI();
    const int *d_send_ldof = mfem::Read(pfes->send_face_nbr_ldof.GetJMemory(),

fem/pgridfunc.hpp

@@ -38,6 +38,11 @@ protected:
        initialized by ExchangeFaceNbrData(). */
    Vector face_nbr_data;
 
+   /** @brief Vector used as an MPI buffer to send face-neighbor data
+        in ExchangeFaceNbrData() to neighboring processors. */
+   //TODO: Use temporary memory to avoid CUDA malloc allocation cost.
+   Vector send_data;
+
    void ProjectBdrCoefficient(Coefficient *coeff[], VectorCoefficient *vcoeff,
                               Array<int> &attr);
 

fem/quadinterpolator.cpp

@@ -150,6 +150,7 @@ void QuadratureInterpolator::Eval3D(
    const int nq = maps.nqpt;
    const int ND = T_ND ? T_ND : nd;
    const int NQ = T_NQ ? T_NQ : nq;
+   const int NMAX = NQ > ND ? NQ : ND;
    const int VDIM = T_VDIM ? T_VDIM : vdim;
    MFEM_VERIFY(ND <= MAX_ND3D, "");
    MFEM_VERIFY(NQ <= MAX_NQ3D, "");
@@ -160,22 +161,24 @@ void QuadratureInterpolator::Eval3D(
    auto val = Reshape(q_val.Write(), NQ, VDIM, NE);
    auto der = Reshape(q_der.Write(), NQ, VDIM, 3, NE);
    auto det = Reshape(q_det.Write(), NQ, NE);
-   MFEM_FORALL(e, NE,
+   MFEM_FORALL_2D(e, NE, NMAX, 1, 1,
    {
       const int ND = T_ND ? T_ND : nd;
       const int NQ = T_NQ ? T_NQ : nq;
       const int VDIM = T_VDIM ? T_VDIM : vdim;
       constexpr int max_ND = T_ND ? T_ND : MAX_ND3D;
       constexpr int max_VDIM = T_VDIM ? T_VDIM : MAX_VDIM3D;
-      double s_E[max_VDIM*max_ND];
-      for (int d = 0; d < ND; d++)
+      MFEM_SHARED double s_E[max_VDIM*max_ND];
+      MFEM_FOREACH_THREAD(d, x, ND)
       {
          for (int c = 0; c < VDIM; c++)
          {
             s_E[c+d*VDIM] = E(d,c,e);
          }
       }
-      for (int q = 0; q < NQ; ++q)
+      MFEM_SYNC_THREAD;
+
+      MFEM_FOREACH_THREAD(q, x, NQ)
       {
          if (eval_flags & VALUES)
          {

fem/quadinterpolator_face.cpp

@@ -495,8 +495,9 @@ void FaceQuadratureInterpolator::Mult(
    }
 }
 
-void FaceQuadratureInterpolator::Values(
-   const Vector &e_vec, Vector &q_val) const
+
+void FaceQuadratureInterpolator::Values(const Vector &e_vec,
+                                        Vector &q_val) const
 {
    Vector q_der, q_det, q_nor;
    Mult(e_vec, VALUES, q_val, q_der, q_det, q_nor);

fem/tmop.cpp

@@ -1058,6 +1058,7 @@ void DiscreteAdaptTC::SetDiscreteTargetBase(const GridFunction &tspec_)
    // make a copy of tspec->tspec_temp, increase its size, and
    // copy data from tspec_temp -> tspec, then add new entries
    Vector tspec_temp = tspec;
+   tspec.UseDevice(true);
    tspec.SetSize(ncomp*dof_cnt);
 
    for (int i = 0; i < tspec_temp.Size(); i++)
@@ -1214,6 +1215,7 @@ void DiscreteAdaptTC::ComputeElementTargets(int e_id, const FiniteElement &fe,
          Array<int> dofs;
          DenseMatrix D_rho(dim), Q_phi(dim), R_theta(dim);
          tspec_fesv->GetElementVDofs(e_id, dofs);
+         tspec.UseDevice(true);
          tspec.GetSubVector(dofs, tspec_vals);
 
          for (int i = 0; i < ir.GetNPoints(); i++)

fem/tmop_tools.cpp

@@ -33,10 +33,11 @@ void AdvectorCG::ComputeAtNewPosition(const Vector &new_nodes,
    const int pnt_cnt = new_field.Size()/ncomp;
 
    new_field = field0;
-
+   new_field.HostReadWrite();
+   Vector new_field_temp;
    for (int i = 0; i < ncomp; i++)
    {
-      Vector new_field_temp(new_field.GetData()+i*pnt_cnt, pnt_cnt);
+      new_field_temp.MakeRef(new_field, i*pnt_cnt, pnt_cnt);
       ComputeAtNewPositionScalar(new_nodes, new_field_temp);
    }
 
@@ -94,6 +95,7 @@ void AdvectorCG::ComputeAtNewPositionScalar(const Vector &new_nodes,
    double v_max = 0.0;
    const int s = new_field.Size();
 
+   u.HostReadWrite();
    for (int i = 0; i < s; i++)
    {
       double vel = 0.;
@@ -149,6 +151,7 @@ void AdvectorCG::ComputeAtNewPositionScalar(const Vector &new_nodes,
 #endif
 
    // Trim the overshoots and undershoots.
+   new_field.HostReadWrite();
    for (int i = 0; i < s; i++)
    {
       if (new_field(i) < glob_minv) { new_field(i) = glob_minv; }

general/device.cpp

@@ -61,15 +61,11 @@ Device Device::device_singleton;
 bool Device::device_env = false;
 bool Device::mem_host_env = false;
 bool Device::mem_device_env = false;
+#ifdef MFEM_USE_UMPIRE
+bool Device::use_umpire = true;
+#endif
 
-Device::Device() : mode(Device::SEQUENTIAL),
-   backends(Backend::CPU),
-   destroy_mm(false),
-   mpi_gpu_aware(false),
-   host_mem_type(MemoryType::HOST),
-   host_mem_class(MemoryClass::HOST),
-   device_mem_type(MemoryType::HOST),
-   device_mem_class(MemoryClass::HOST)
+Device::Device()
 {
    if (getenv("MFEM_MEMORY") && !mem_host_env && !mem_device_env)
    {
@@ -135,7 +131,7 @@ Device::Device() : mode(Device::SEQUENTIAL),
       {
          MFEM_ABORT("Unknown memory backend!");
       }
-      mm.Configure(host_mem_type, device_mem_type);
+      mm.Configure(host_mem_type, device_mem_type, device_mem_type);
    }
 
    if (getenv("MFEM_DEVICE"))
@@ -165,6 +161,8 @@ Device::~Device()
    Get().host_mem_class = MemoryClass::HOST;
    Get().device_mem_type = MemoryType::HOST;
    Get().device_mem_class = MemoryClass::HOST;
+   Get().device_temp_mem_type = MemoryType::HOST;
+   Get().device_temp_mem_class = MemoryClass::HOST;
 }
 
 void Device::Configure(const std::string &device, const int dev)
@@ -260,6 +258,10 @@ void Device::Print(std::ostream &out)
    if (Device::Allows(Backend::DEVICE_MASK))
    {
       out << ',' << MemoryTypeName[static_cast<int>(device_mem_type)];
+      if (device_temp_mem_type != device_mem_type)
+      {
+         out << ',' << MemoryTypeName[static_cast<int>(device_temp_mem_type)];
+      }
    }
    out << std::endl;
 }
@@ -272,7 +274,8 @@ void Device::UpdateMemoryTypeAndClass()
 
 #ifdef MFEM_USE_UMPIRE
    // If MFEM has been compiled with Umpire support, use it as the default
-   if (!mem_host_env) { host_mem_type = MemoryType::HOST_UMPIRE; }
+   // TODO TMS: temporary
+   //if (!mem_host_env && use_umpire) { host_mem_type = MemoryType::HOST_UMPIRE; }
 #endif
 
    // Enable the device memory type
@@ -296,11 +299,16 @@ void Device::UpdateMemoryTypeAndClass()
          }
          else
          {
-#ifndef MFEM_USE_UMPIRE
-            device_mem_type = MemoryType::DEVICE;
-#else
-            device_mem_type = MemoryType::DEVICE_UMPIRE;
+#ifdef MFEM_USE_UMPIRE
+            if (use_umpire)
+            {
+               device_mem_type = MemoryType::DEVICE_UMPIRE;
+            }
+            else
 #endif
+            {
+               device_mem_type = MemoryType::DEVICE;
+            }
          }
       }
       device_mem_class = MemoryClass::DEVICE;
@@ -320,8 +328,20 @@ void Device::UpdateMemoryTypeAndClass()
       device_mem_type = MemoryType::DEVICE_DEBUG;
    }
 
+   // Setup device_temp_mem_{type,class}
+   switch (device_mem_type)
+   {
+      case MemoryType::DEVICE_UMPIRE:
+         device_temp_mem_type = MemoryType::DEVICE_TEMP_UMPIRE;
+         break;
+      default:
+         device_temp_mem_type = device_mem_type;
+         break;
+   }
+   device_temp_mem_class = device_mem_class;
+
    // Update the memory manager with the new settings
-   mm.Configure(host_mem_type, device_mem_type);
+   mm.Configure(host_mem_type, device_mem_type, device_temp_mem_type);
 }
 
 void Device::Enable()

general/device.hpp

@@ -119,20 +119,26 @@ private:
 
    static bool device_env, mem_host_env, mem_device_env;
    static Device device_singleton;
+#ifdef MFEM_USE_UMPIRE
+   static bool use_umpire;
+#endif
 
-   MODES mode;
+   MODES mode{Device::SEQUENTIAL};
    int dev = 0; ///< Device ID of the configured device.
    int ngpu = -1; ///< Number of detected devices; -1: not initialized.
-   unsigned long backends; ///< Bitwise-OR of all configured backends.
+   unsigned long backends{Backend::CPU}; ///< Bitwise-OR of all configured backends.
    /// Set to true during configuration, except in 'device_singleton'.
-   bool destroy_mm;
-   bool mpi_gpu_aware;
+   bool destroy_mm{false};
+   bool mpi_gpu_aware{false};
 
-   MemoryType host_mem_type;      ///< Current Host MemoryType
-   MemoryClass host_mem_class;    ///< Current Host MemoryClass
+   MemoryType host_mem_type{MemoryType::HOST};      ///< Current Host MemoryType
+   MemoryClass host_mem_class{MemoryClass::HOST};    ///< Current Host MemoryClass
 
-   MemoryType device_mem_type;    ///< Current Device MemoryType
-   MemoryClass device_mem_class;  ///< Current Device MemoryClass
+   MemoryType device_mem_type{MemoryType::HOST};    ///< Current Device MemoryType
+   MemoryClass device_mem_class{MemoryClass::HOST};  ///< Current Device MemoryClass
+
+   MemoryType device_temp_mem_type{MemoryType::HOST};    ///< Current Device MemoryType
+   MemoryClass device_temp_mem_class{MemoryClass::HOST};  ///< Current Device MemoryClass
 
    char *device_option = NULL;
    Device(Device const&);
@@ -173,14 +179,6 @@ public:
        @note This object should be destroyed after all other MFEM objects that
        use the Device are destroyed. */
    Device(const std::string &device, const int dev = 0)
-      : mode(Device::SEQUENTIAL),
-        backends(Backend::CPU),
-        destroy_mm(false),
-        mpi_gpu_aware(false),
-        host_mem_type(MemoryType::HOST),
-        host_mem_class(MemoryClass::HOST),
-        device_mem_type(MemoryType::HOST),
-        device_mem_class(MemoryClass::HOST)
    { Configure(device, dev); }
 
    /// Destructor.
@@ -260,10 +258,24 @@ public:
    /** @deprecated Use GetDeviceMemoryClass() instead. */
    static inline MemoryClass GetMemoryClass() { return Get().device_mem_class; }
 
+   /** @brief Get the current Device Temporary MemoryType. This is the MemoryType used by
+       MFEM classes when allocating temporary memory to be used with device kernels.
+   */
+   static inline MemoryType GetDeviceTempMemoryType() { return Get().device_temp_mem_type; }
+
+   /** @brief Get the current Device Temporary MemoryClass. This is the MemoryClass used
+       by MFEM device kernels when they need to access temporary Memory objects. */
+   static inline MemoryClass GetDeviceTempMemoryClass() { return Get().device_temp_mem_class; }
+
    static void SetGPUAwareMPI(const bool force = true)
    { Get().mpi_gpu_aware = force; }
 
    static bool GetGPUAwareMPI() { return Get().mpi_gpu_aware; }
+
+#ifdef MFEM_USE_UMPIRE
+   static bool UseUmpire() { return Get().use_umpire; }
+   static void UseUmpire(bool use) { Get().use_umpire = use; }
+#endif
 };
 
 

general/mem_manager.cpp

@@ -67,15 +67,19 @@ MemoryType MemoryManager::GetDualMemoryType_(MemoryType mt)
 {
    switch (mt)
    {
-      case MemoryType::HOST:           return MemoryType::DEVICE;
+      // TODO TMS: temporary
+      case MemoryType::HOST:           return MemoryType::DEVICE_UMPIRE;
       case MemoryType::HOST_32:        return MemoryType::DEVICE;
       case MemoryType::HOST_64:        return MemoryType::DEVICE;
       case MemoryType::HOST_DEBUG:     return MemoryType::DEVICE_DEBUG;
-      case MemoryType::HOST_UMPIRE:    return MemoryType::DEVICE_UMPIRE;
+      //case MemoryType::HOST_UMPIRE:    return MemoryType::DEVICE_UMPIRE;
       case MemoryType::MANAGED:        return MemoryType::MANAGED;
       case MemoryType::DEVICE:         return MemoryType::HOST;
       case MemoryType::DEVICE_DEBUG:   return MemoryType::HOST_DEBUG;
-      case MemoryType::DEVICE_UMPIRE:  return MemoryType::HOST_UMPIRE;
+      //case MemoryType::DEVICE_UMPIRE:  return MemoryType::HOST_UMPIRE;
+      case MemoryType::DEVICE_UMPIRE:  return MemoryType::HOST;
+      //case MemoryType::DEVICE_TEMP_UMPIRE:  return MemoryType::HOST_UMPIRE;
+      case MemoryType::DEVICE_TEMP_UMPIRE:  return MemoryType::HOST;
       default: mfem_error("Unknown memory type!");
    }
    MFEM_VERIFY(false,"");
@@ -88,6 +92,9 @@ static void MFEM_VERIFY_TYPES(const MemoryType h_mt, const MemoryType d_mt)
    MFEM_ASSERT(IsDeviceMemory(d_mt),"");
    const bool sync =
       (h_mt == MemoryType::HOST_UMPIRE && d_mt == MemoryType::DEVICE_UMPIRE) ||
+      (h_mt == MemoryType::HOST_UMPIRE && d_mt == MemoryType::DEVICE_TEMP_UMPIRE) ||
+      (h_mt == MemoryType::HOST && d_mt == MemoryType::DEVICE_UMPIRE) ||
+      (h_mt == MemoryType::HOST && d_mt == MemoryType::DEVICE_TEMP_UMPIRE) ||
       (h_mt == MemoryType::HOST_DEBUG && d_mt == MemoryType::DEVICE_DEBUG) ||
       (h_mt == MemoryType::MANAGED && d_mt == MemoryType::MANAGED) ||
       (h_mt == MemoryType::HOST_64 && d_mt == MemoryType::DEVICE) ||
@@ -461,48 +468,96 @@ public:
 #ifndef MFEM_USE_UMPIRE
 class UmpireHostMemorySpace : public NoHostMemorySpace { };
 class UmpireDeviceMemorySpace : public NoDeviceMemorySpace { };
+class UmpireDeviceTempMemorySpace : public NoDeviceMemorySpace { };
 #else
+
+// TODO TMS: replace with um.hasAllocatorId(int) when it exists
+bool UmpireHasId(const umpire::ResourceManager & rm, int id)
+{
+   const auto & ids = rm.getAllocatorIds();
+   return std::find(ids.begin(), ids.end(), id) != ids.end();
+}
+
 /// The Umpire host memory space
 class UmpireHostMemorySpace : public HostMemorySpace
 {
 private:
-   const char *name;
    umpire::ResourceManager &rm;
    umpire::Allocator h_allocator;
-   umpire::strategy::AllocationStrategy *strat;
+   bool owns_allocator{false};
 public:
-   ~UmpireHostMemorySpace() { h_allocator.release(); }
-   UmpireHostMemorySpace():
-      HostMemorySpace(),
-      name(mm.GetUmpireAllocatorHostName()),
-      rm(umpire::ResourceManager::getInstance()),
-      h_allocator(rm.isAllocator(name)? rm.getAllocator(name):
-                  rm.makeAllocator<umpire::strategy::DynamicPool>
-                  (name, rm.getAllocator("HOST"))),
-      strat(h_allocator.getAllocationStrategy()) { }
+   // TODO: this only releases unused memory
+   ~UmpireHostMemorySpace() { if (owns_allocator) { h_allocator.release(); } }
+   UmpireHostMemorySpace(): HostMemorySpace(),
+      rm(umpire::ResourceManager::getInstance())
+   {
+      const int id = MemoryManager::GetUmpireHostAllocatorId();
+      if (!UmpireHasId(rm, id))
+      {
+         h_allocator = rm.makeAllocator<umpire::strategy::DynamicPool>("MFEM_HOST",
+                                                                       rm.getAllocator("HOST"));
+         owns_allocator = true;
+      }
+      else
+      {
+         h_allocator = rm.getAllocator(id);
+      }
+      MemoryManager::SetUmpireHostAllocatorId(id);
+   }
    void Alloc(void **ptr, size_t bytes) { *ptr = h_allocator.allocate(bytes); }
    void Dealloc(void *ptr) { h_allocator.deallocate(ptr); }
    void Insert(void *ptr, size_t bytes)
-   { rm.registerAllocation(ptr, {ptr, bytes, strat}); }
+   { mfem_error("UmpireHostMemorySpace::Insert is unsupported"); }
 };
 
 /// The Umpire device memory space
 #ifdef MFEM_USE_CUDA
-class UmpireDeviceMemorySpace : public DeviceMemorySpace
+class UmpireDeviceMemorySpaceImpl : public DeviceMemorySpace
 {
+public:
+   enum class AllocatorType { TEMPORARY, PERMANENT };
 private:
-   const char *name;
    umpire::ResourceManager &rm;
    umpire::Allocator d_allocator;
+   bool owns_allocator{false};
+
+   int SetupAllocator(int possible_id, const char * allocator_name)
+   {
+      if (!UmpireHasId(rm, possible_id))
+      {
+         d_allocator = rm.makeAllocator<umpire::strategy::DynamicPool>(allocator_name,
+                                                                       rm.getAllocator("DEVICE"));
+         owns_allocator = true;
+      }
+      else
+      {
+         d_allocator = rm.getAllocator(possible_id);
+      }
+
+      return d_allocator.getId();
+   }
 public:
-   ~UmpireDeviceMemorySpace() { d_allocator.release(); }
-   UmpireDeviceMemorySpace():
-      DeviceMemorySpace(),
-      name(mm.GetUmpireAllocatorDeviceName()),
-      rm(umpire::ResourceManager::getInstance()),
-      d_allocator(rm.isAllocator(name)? rm.getAllocator(name):
-                  rm.makeAllocator<umpire::strategy::DynamicPool>
-                  (name, rm.getAllocator("DEVICE"))) { }
+   // TODO: this only releases unused memory
+   ~UmpireDeviceMemorySpaceImpl() { if (owns_allocator) { d_allocator.release(); } }
+   UmpireDeviceMemorySpaceImpl(AllocatorType t): DeviceMemorySpace(),
+      rm(umpire::ResourceManager::getInstance())
+   {
+      switch (t)
+      {
+         case AllocatorType::PERMANENT:
+            MemoryManager::SetUmpireDeviceAllocatorId(SetupAllocator(
+                                                         MemoryManager::GetUmpireDeviceAllocatorId(),
+                                                         "MFEM_DEVICE"));
+            break;
+         case AllocatorType::TEMPORARY:
+            MemoryManager::SetUmpireDeviceTempAllocatorId(SetupAllocator(
+                                                             MemoryManager::GetUmpireDeviceTempAllocatorId(),
+                                                             "MFEM_DEVICE_TEMPORARY"));
+            break;
+         default:
+            mfem_error("Unknown Umpire AllocatorType");
+      }
+   }
    void Alloc(Memory &base) { base.d_ptr = d_allocator.allocate(base.bytes); }
    void Dealloc(Memory &base) { d_allocator.deallocate(base.d_ptr); }
    void *HtoD(void *dst, const void *src, size_t bytes)
@@ -536,8 +591,23 @@ public:
       //rm.copy(dst, const_cast<void*>(src), bytes); return dst;
    }
 };
+
+class UmpireDeviceMemorySpace : public UmpireDeviceMemorySpaceImpl
+{
+public:
+   UmpireDeviceMemorySpace() : UmpireDeviceMemorySpaceImpl(
+         AllocatorType::PERMANENT) {}
+};
+
+class UmpireDeviceTempMemorySpace : public UmpireDeviceMemorySpaceImpl
+{
+public:
+   UmpireDeviceTempMemorySpace() : UmpireDeviceMemorySpaceImpl(
+         AllocatorType::TEMPORARY) {}
+};
 #else
 class UmpireDeviceMemorySpace : public NoDeviceMemorySpace { };
+class UmpireDeviceTempMemorySpace : public NoDeviceMemorySpace { };
 #endif // MFEM_USE_CUDA
 #endif // MFEM_USE_UMPIRE
 
@@ -568,7 +638,7 @@ public:
       host[static_cast<int>(MT::HOST_64)] = new Aligned64HostMemorySpace();
       // HOST_DEBUG is delayed, as it reroutes signals
       host[static_cast<int>(MT::HOST_DEBUG)] = nullptr;
-      host[static_cast<int>(MT::HOST_UMPIRE)] = new UmpireHostMemorySpace();
+      host[static_cast<int>(MT::HOST_UMPIRE)] = nullptr;
       host[static_cast<int>(MT::MANAGED)] = new UvmHostMemorySpace();
 
       // Filling the device memory backends, shifting with the device size
@@ -610,8 +680,12 @@ public:
 private:
    HostMemorySpace* NewHostCtrl(const MemoryType mt)
    {
-      if (mt == MT::HOST_DEBUG) { return new MmuHostMemorySpace(); }
-      MFEM_ABORT("Unknown host memory controller!");
+      switch (mt)
+      {
+         case MT::HOST_DEBUG: return new MmuHostMemorySpace();
+         case MT::HOST_UMPIRE: return new UmpireHostMemorySpace();
+         default: MFEM_ABORT("Unknown host memory controller!");
+      }
       return nullptr;
    }
 
@@ -620,6 +694,7 @@ private:
       switch (mt)
       {
          case MT::DEVICE_UMPIRE: return new UmpireDeviceMemorySpace();
+         case MT::DEVICE_TEMP_UMPIRE: return new UmpireDeviceTempMemorySpace();
          case MT::DEVICE_DEBUG: return new MmuDeviceMemorySpace();
          case MT::DEVICE:
          {
@@ -760,7 +835,7 @@ bool MemoryManager::MemoryClassCheck_(MemoryClass mc, void *h_ptr,
    const bool known = mm.IsKnown(h_ptr);
    const bool alias = mm.IsAlias(h_ptr);
    const bool check = known || ((flags & Mem::ALIAS) && alias);
-   MFEM_VERIFY(check,"");
+   MFEM_VERIFY(check,"Unknown host pointer: " << h_ptr);
    const internal::Memory &mem =
       (flags & Mem::ALIAS) ?
       *maps->aliases.at(h_ptr).mem : maps->memories.at(h_ptr);
@@ -783,6 +858,7 @@ bool MemoryManager::MemoryClassCheck_(MemoryClass mc, void *h_ptr,
          MFEM_VERIFY(d_mt == MemoryType::DEVICE ||
                      d_mt == MemoryType::DEVICE_DEBUG ||
                      d_mt == MemoryType::DEVICE_UMPIRE ||
+                     d_mt == MemoryType::DEVICE_TEMP_UMPIRE ||
                      d_mt == MemoryType::MANAGED,"");
          return true;
       }
@@ -1262,22 +1338,15 @@ MemoryManager::MemoryManager() { Init(); }
 MemoryManager::~MemoryManager() { if (exists) { Destroy(); } }
 
 void MemoryManager::Configure(const MemoryType host_mt,
-                              const MemoryType device_mt)
+                              const MemoryType device_mt,
+                              const MemoryType device_tmt)
 {
    Init();
    host_mem_type = host_mt;
    device_mem_type = device_mt;
+   device_temp_mem_type = device_tmt;
 }
 
-#ifdef MFEM_USE_UMPIRE
-void MemoryManager::SetUmpireAllocatorNames(const char *h_name,
-                                            const char *d_name)
-{
-   h_umpire_name = h_name;
-   d_umpire_name = d_name;
-}
-#endif
-
 void MemoryManager::Destroy()
 {
    MFEM_VERIFY(exists, "MemoryManager has already been destroyed!");
@@ -1381,12 +1450,14 @@ MemoryManager mm;
 bool MemoryManager::exists = false;
 
 #ifdef MFEM_USE_UMPIRE
-const char* MemoryManager::h_umpire_name = "HOST";
-const char* MemoryManager::d_umpire_name = "DEVICE";
+int MemoryManager::h_umpire_id = -1;
+int MemoryManager::d_umpire_id = -1;
+int MemoryManager::d_umpire_temp_id = -1;
 #endif
 
 MemoryType MemoryManager::host_mem_type = MemoryType::HOST;
 MemoryType MemoryManager::device_mem_type = MemoryType::HOST;
+MemoryType MemoryManager::device_temp_mem_type = MemoryType::HOST;
 
 const char *MemoryTypeName[MemoryTypeSize] =
 {
@@ -1403,11 +1474,14 @@ const char *MemoryTypeName[MemoryTypeSize] =
 #endif
    "device-debug",
 #if defined(MFEM_USE_CUDA)
-   "cuda-umpire"
+   "cuda-umpire",
+   "cuda-umpire-temp"
 #elif defined(MFEM_USE_HIP)
-   "hip-umpire"
+   "hip-umpire",
+   "hip-umpire-temp"
 #else
-   "device-umpire"
+   "device-umpire",
+   "device-umpire-temp"
 #endif
 };
 

general/mem_manager.hpp

@@ -27,18 +27,19 @@ namespace mfem
 /// Memory types supported by MFEM.
 enum class MemoryType
 {
-   HOST,           ///< Host memory; using new[] and delete[]
-   HOST_32,        ///< Host memory; aligned at 32 bytes
-   HOST_64,        ///< Host memory; aligned at 64 bytes
-   HOST_DEBUG,     ///< Host memory; allocated from a "host-debug" pool
-   HOST_UMPIRE,    ///< Host memory; using Umpire
-   MANAGED,        /**< Managed memory; using CUDA or HIP *MallocManaged
-                        and *Free */
-   DEVICE,         ///< Device memory; using CUDA or HIP *Malloc and *Free
-   DEVICE_DEBUG,   /**< Pseudo-device memory; allocated on host from a
-                        "device-debug" pool */
-   DEVICE_UMPIRE,  ///< Device memory; using Umpire
-   SIZE            ///< Number of host and device memory types
+   HOST,                ///< Host memory; using new[] and delete[]
+   HOST_32,             ///< Host memory; aligned at 32 bytes
+   HOST_64,             ///< Host memory; aligned at 64 bytes
+   HOST_DEBUG,          ///< Host memory; allocated from a "host-debug" pool
+   HOST_UMPIRE,         ///< Host memory; using Umpire
+   MANAGED,             /**< Managed memory; using CUDA or HIP *MallocManaged
+                             and *Free */
+   DEVICE,              ///< Device memory; using CUDA or HIP *Malloc and *Free
+   DEVICE_DEBUG,        /**< Pseudo-device memory; allocated on host from a
+                             "device-debug" pool */
+   DEVICE_UMPIRE,       ///< Device memory; using Umpire
+   DEVICE_TEMP_UMPIRE,  ///< Temporary Device memory; using Umpire
+   SIZE                 ///< Number of host and device memory types
 };
 
 /// Static casts to 'int' and sizes of some useful memory types.
@@ -61,7 +62,7 @@ enum class MemoryClass
                                  HOST_UMPIRE, MANAGED } */
    HOST_32, ///< Memory types: { HOST_32, HOST_64, HOST_DEBUG }
    HOST_64, ///< Memory types: { HOST_64, HOST_DEBUG }
-   DEVICE,  ///< Memory types: { DEVICE, DEVICE_DEBUG, DEVICE_UMPIRE, MANAGED }
+   DEVICE,  ///< Memory types: { DEVICE, DEVICE_DEBUG, DEVICE_UMPIRE, DEVICE_TEMP_UMPIRE, MANAGED }
    MANAGED  ///< Memory types: { MANAGED }
 };
 
@@ -450,16 +451,16 @@ private:
 
    template <std::size_t align_bytes, bool dummy = true> struct Alloc
    {
+#if __cplusplus < 201703L
       static inline T *New(std::size_t)
       {
-#if __cplusplus < 201703L
          // Generate an error in debug mode
          MFEM_ASSERT(false, "overaligned type cannot use MemoryType::HOST");
          return nullptr;
-#else
-         return new T[size];
-#endif
       }
+#else
+      static inline T *New(std::size_t size) { return new T[size]; }
+#endif
    };
 
 #if __cplusplus < 201703L
@@ -489,6 +490,9 @@ private:
    /// Device memory type set during the Setup.
    static MemoryType device_mem_type;
 
+   /// Device temporary memory type set during the Setup.
+   static MemoryType device_temp_mem_type;
+
    /// Allow to detect if a global memory manager instance exists.
    static bool exists;
 
@@ -497,8 +501,9 @@ private:
 
    /// Host and device allocator names for Umpire.
 #ifdef MFEM_USE_UMPIRE
-   static const char *h_umpire_name;
-   static const char *d_umpire_name;
+   static int h_umpire_id;
+   static int d_umpire_id;
+   static int d_umpire_temp_id;
 #endif
 
 private: // Static methods used by the Memory<T> class
@@ -621,15 +626,21 @@ public:
    /// Initialize the memory manager.
    void Init();
 
-   /// Configure the Memory manager with given default host and device types
+   /// Configure the Memory manager with given default host, device, and device temporary types
    /// This method will be called when configuring a device.
-   void Configure(const MemoryType h_mt, const MemoryType d_mt);
+   void Configure(const MemoryType h_mt, const MemoryType d_mt,
+                  const MemoryType d_tmt);
 
 #ifdef MFEM_USE_UMPIRE
-   /// Set the host and device UMpire allocator names
-   void SetUmpireAllocatorNames(const char *h_name, const char *d_name);
-   const char *GetUmpireAllocatorHostName() { return h_umpire_name; }
-   const char *GetUmpireAllocatorDeviceName() { return d_umpire_name; }
+   /// Set the host and device Umpire allocator ids
+   static void SetUmpireHostAllocatorId(int h_id) { h_umpire_id = h_id; }
+   static void SetUmpireDeviceAllocatorId(int d_id) { d_umpire_id = d_id; }
+   static void SetUmpireDeviceTempAllocatorId(int d_id) { d_umpire_temp_id = d_id; }
+
+   /// Get the host and device Umpire allocator ids
+   static int GetUmpireHostAllocatorId() { return h_umpire_id; }
+   static int GetUmpireDeviceAllocatorId() { return d_umpire_id; }
+   static int GetUmpireDeviceTempAllocatorId() { return d_umpire_temp_id; }
 #endif
 
    /// Free all the device memories
@@ -654,6 +665,7 @@ public:
 
    static MemoryType GetHostMemoryType() { return host_mem_type; }
    static MemoryType GetDeviceMemoryType() { return device_mem_type; }
+   static MemoryType GetDeviceTempMemoryType() { return device_temp_mem_type; }
 };
 
 

linalg/densemat.cpp

@@ -3595,4 +3595,114 @@ void BatchLUSolve(const DenseTensor &Mlu, const Array<int> &P, Vector &X)
 
 }
 
+void BatchLUFactor(Vector &Minv,const int m,const int NE, Array<int> &P)
+{
+   P.SetSize(m*NE);
+   auto data_all = mfem::Reshape(Minv.ReadWrite(), m, m, NE);
+   auto piv_all = mfem::Reshape(P.Write(), m, NE);
+
+   MFEM_FORALL(e, NE,
+   {
+
+      double *data = &data_all(0,0,e);
+      int *ipiv = &piv_all(0,e);
+      for (int i = 0; i < m; i++)
+      {
+
+         // pivoting
+         {
+            int piv = i;
+            double a = fabs(data[piv+i*m]);
+            for (int j = i+1; j < m; j++)
+            {
+               const double b = fabs(data[j+i*m]);
+               if (b > a)
+               {
+                  a = b;
+                  piv = j;
+               }
+            }
+            ipiv[i] = piv;
+            if (piv != i)
+            {
+               // swap rows i and piv in both L and U parts
+               for (int j = 0; j < m; j++)
+               {
+                  mfem::kernels::internal::Swap<double>(data[i+j*m], data[piv+j*m]);
+               }
+            }
+         }//pivot end
+
+         //Q: How to check for errors?
+         //if (abs(data[i + i*m]) <= TOL)
+         //{
+         //return false; // failed
+         //}
+
+         const double a_ii_inv = 1.0 / data[i+i*m];
+         for (int j = i+1; j < m; j++)
+         {
+            data[j+i*m] *= a_ii_inv;
+         }
+
+         for (int k = i+1; k < m; k++)
+         {
+            const double a_ik = data[i+k*m];
+            for (int j = i+1; j < m; j++)
+            {
+               data[j+k*m] -= a_ik * data[j+i*m];
+            }
+         }
+
+      }//m loop
+
+   });
+
+}
+
+void BatchLUSolve(Vector &Minv, int m, int NE,
+                  Array<int> &P, Vector &X)
+{
+
+   auto data_all = mfem::Reshape(Minv.Read(), m, m, NE);
+   auto piv_all = mfem::Reshape(P.Read(), m, NE);
+   auto x_all = mfem::Reshape(X.ReadWrite(), m, NE);
+
+   MFEM_FORALL(e, NE,
+   {
+
+      const double *data = &data_all(0,0,e);
+      const int *ipiv = &piv_all(0,e);
+      double *x = &x_all(0,e);
+
+      // X <- P X
+      for (int i = 0; i < m; i++)
+      {
+         mfem::kernels::internal::Swap<double>(x[i], x[ipiv[i]]);
+      }
+
+      // X <- L^{-1} X
+      for (int j = 0; j < m; j++)
+      {
+         const double x_j = x[j];
+         for (int i = j+1; i < m; i++)
+         {
+            x[i] -= data[i+j*m] * x_j;
+         }
+      }
+
+      // X <- U^{-1} X
+      for (int j = m-1; j >= 0; j--)
+      {
+         const double x_j = ( x[j] /= data[j+j*m] );
+         for (int i = 0; i < j; i++)
+         {
+            x[i] -= data[i+j*m] * x_j;
+         }
+      }
+   });
+
+}
+
+
 } // namespace mfem

linalg/densemat.hpp

@@ -884,6 +884,10 @@ void BatchLUFactor(DenseTensor &Mlu, Array<int> &P, const double TOL = 0.0);
     dimension m x n. */
 void BatchLUSolve(const DenseTensor &Mlu, const Array<int> &P, Vector &X);
 
+void BatchLUFactor(Vector &Minv,int m,int NE, Array<int> &P);
+
+void BatchLUSolve(Vector &Minv, int m, int NE,
+                  Array<int> &P, Vector &X);
 
 // Inline methods
 

linalg/sparsemat.cpp

@@ -28,6 +28,25 @@ namespace mfem
 
 using namespace std;
 
+#ifdef MFEM_USE_CUDA
+int SparseMatrix::SparseMatrixCount = 0;
+cusparseHandle_t SparseMatrix::handle;
+size_t SparseMatrix::bufferSize = 0;
+void * SparseMatrix::dBuffer = nullptr;
+#endif
+
+void SparseMatrix::InitCuSparse()
+{
+   /* Initialize CuSparse library */
+#ifdef MFEM_USE_CUDA
+   SparseMatrixCount++;
+   if (SparseMatrixCount == 1 && Device::Allows(Backend::CUDA_MASK))
+   {
+      cusparseCreate(&handle);
+   }
+#endif
+}
+
 SparseMatrix::SparseMatrix(int nrows, int ncols)
    : AbstractSparseMatrix(nrows, (ncols >= 0) ? ncols : nrows),
      Rows(new RowNode *[nrows]),
@@ -50,6 +69,8 @@ SparseMatrix::SparseMatrix(int nrows, int ncols)
 #ifdef MFEM_USE_MEMALLOC
    NodesMem = new RowNodeAlloc;
 #endif
+
+   InitCuSparse();
 }
 
 SparseMatrix::SparseMatrix(int *i, int *j, double *data, int m, int n)
@@ -67,6 +88,8 @@ SparseMatrix::SparseMatrix(int *i, int *j, double *data, int m, int n)
 #ifdef MFEM_USE_MEMALLOC
    NodesMem = NULL;
 #endif
+
+   InitCuSparse();
 }
 
 SparseMatrix::SparseMatrix(int *i, int *j, double *data, int m, int n,
@@ -98,6 +121,8 @@ SparseMatrix::SparseMatrix(int *i, int *j, double *data, int m, int n,
          A[i] = 0.0;
       }
    }
+
+   InitCuSparse();
 }
 
 SparseMatrix::SparseMatrix(int nrows, int ncols, int rowsize)
@@ -119,9 +144,11 @@ SparseMatrix::SparseMatrix(int nrows, int ncols, int rowsize)
    {
       I[i] = i * rowsize;
    }
+
+   InitCuSparse();
 }
 
-SparseMatrix::SparseMatrix(const SparseMatrix &mat, bool copy_graph)
+SparseMatrix::SparseMatrix(const SparseMatrix &mat, bool copy_graph, MemoryType mt)
    : AbstractSparseMatrix(mat.Height(), mat.Width())
 {
    if (mat.Finalized())
@@ -129,8 +156,8 @@ SparseMatrix::SparseMatrix(const SparseMatrix &mat, bool copy_graph)
       const int nnz = mat.I[height];
       if (copy_graph)
       {
-         I.New(height+1, mat.I.GetMemoryType());
-         J.New(nnz, mat.J.GetMemoryType());
+         I.New(height+1, mt == MemoryType::SIZE ? mat.I.GetMemoryType() : mt);
+         J.New(nnz, mt == MemoryType::SIZE ? mat.J.GetMemoryType() : mt);
          I.CopyFrom(mat.I, height+1);
          J.CopyFrom(mat.J, nnz);
       }
@@ -141,7 +168,7 @@ SparseMatrix::SparseMatrix(const SparseMatrix &mat, bool copy_graph)
          I.ClearOwnerFlags();
          J.ClearOwnerFlags();
       }
-      A.New(nnz, mat.A.GetMemoryType());
+      A.New(nnz, mt == MemoryType::SIZE ? mat.A.GetMemoryType() : mt);
       A.CopyFrom(mat.A, nnz);
 
       Rows = NULL;
@@ -184,6 +211,8 @@ SparseMatrix::SparseMatrix(const SparseMatrix &mat, bool copy_graph)
    ColPtrNode = NULL;
    At = NULL;
    isSorted = mat.isSorted;
+
+   InitCuSparse();
 }
 
 SparseMatrix::SparseMatrix(const Vector &v)
@@ -211,6 +240,8 @@ SparseMatrix::SparseMatrix(const Vector &v)
       J[r] = r;
       A[r] = v[r];
    }
+
+   InitCuSparse();
 }
 
 SparseMatrix& SparseMatrix::operator=(const SparseMatrix &rhs)
@@ -250,6 +281,16 @@ void SparseMatrix::SetEmpty()
    NodesMem = NULL;
 #endif
    isSorted = false;
+
+#ifdef MFEM_USE_CUDA
+   if (initBuffers)
+   {
+      cusparseDestroySpMat(matA_descr);
+      cusparseDestroyDnVec(vecX_descr);
+      cusparseDestroyDnVec(vecY_descr);
+      initBuffers = false;
+   }
+#endif
 }
 
 int SparseMatrix::RowSize(const int i) const
@@ -494,24 +535,29 @@ void SparseMatrix::GetDiag(Vector & d) const
 
    d.SetSize(height);
 
-   int j, end;
-   for (int i = 0; i < height; i++)
-   {
+   auto I = this->ReadI();
+   auto J = this->ReadJ();
+   auto A = this->ReadData();
+   auto dd = d.Write();
 
-      end = I[i+1];
-      for (j = I[i]; j < end; j++)
+   MFEM_FORALL(i, height,
+   {
+      const int begin = I[i];
+      const int end = I[i+1];
+      int j;
+      for (j = begin; j < end; j++)
       {
          if (J[j] == i)
          {
-            d[i] = A[j];
+            dd[i] = A[j];
             break;
          }
       }
       if (j == end)
       {
-         d[i] = 0.;
+         dd[i] = 0.;
       }
-   }
+   });
 }
 
 /// Produces a DenseMatrix from a SparseMatrix
@@ -587,16 +633,72 @@ void SparseMatrix::AddMult(const Vector &x, Vector &y, const double a) const
    auto d_A = Read(A, nnz);
    auto d_x = x.Read();
    auto d_y = y.ReadWrite();
-   MFEM_FORALL(i, height,
+
+   //Skip if matrix has no non-zeros
+   if (nnz == 0) {return;}
+   if (Device::Allows(Backend::CUDA_MASK) && useCuSparse)
    {
-      double d = 0.0;
-      const int end = d_I[i+1];
-      for (int j = d_I[i]; j < end; j++)
+#ifdef MFEM_USE_CUDA
+      const double alpha = a;
+      const double beta  = 1.0;
+
+      //Setup descriptors
+      if (!initBuffers)
       {
-         d += d_A[j] * d_x[d_J[j]];
+         /* Setup matrix descriptor */
+         cusparseCreateCsr(&matA_descr,Height(), Width(), J.Capacity(),
+                           const_cast<int *>(d_I),
+                           const_cast<int *>(d_J), const_cast<double *>(d_A), CUSPARSE_INDEX_32I,
+                           CUSPARSE_INDEX_32I, CUSPARSE_INDEX_BASE_ZERO, CUDA_R_64F);
+
+         /*Create handles for input/output vectors */
+         cusparseCreateDnVec(&vecX_descr, x.Size(), const_cast<double *>(d_x),
+                             CUDA_R_64F);
+         cusparseCreateDnVec(&vecY_descr, y.Size(), d_y, CUDA_R_64F);
+
+         initBuffers = true;
       }
-      d_y[i] += a * d;
-   });
+
+      /*Allocate space for kernel. Buffer is shared between different sparsemats */
+      size_t newBufferSize = 0;
+      cusparseSpMV_bufferSize(handle, CUSPARSE_OPERATION_NON_TRANSPOSE, &alpha,
+                              matA_descr,
+                              vecX_descr, &beta, vecY_descr, CUDA_R_64F,
+                              CUSPARSE_CSRMV_ALG1, &newBufferSize);
+
+      //Check if need to resize
+      if (newBufferSize > bufferSize)
+      {
+         bufferSize = newBufferSize;
+         if (dBuffer != NULL) { CuMemFree(dBuffer); }
+         CuMemAlloc(&dBuffer, bufferSize);
+      }
+
+      //Update input/output vectors
+      cusparseDnVecSetValues(vecX_descr, const_cast<double *>(d_x));
+      cusparseDnVecSetValues(vecY_descr, d_y);
+
+      // Y = alpha A * X + beta * Y
+      cusparseSpMV(handle, CUSPARSE_OPERATION_NON_TRANSPOSE, &alpha, matA_descr,
+                   vecX_descr, &beta, vecY_descr, CUDA_R_64F, CUSPARSE_CSRMV_ALG1, dBuffer);
+#endif
+   }
+   else
+   {
+      //Native version
+      MFEM_FORALL(i, height,
+      {
+         double d = 0.0;
+         const int end = d_I[i+1];
+         for (int j = d_I[i]; j < end; j++)
+         {
+            d += d_A[j] * d_x[d_J[j]];
+         }
+         d_y[i] += a * d;
+      });
+
+   }
+
 #else
    const double *Ap = A, *xp = x.GetData();
    double *yp = y.GetData();
@@ -2145,31 +2247,46 @@ void SparseMatrix::DiagScale(const Vector &b, Vector &x, double sc) const
 {
    MFEM_VERIFY(Finalized(), "Matrix must be finalized.");
 
+   const int nnz = J.Capacity();
+
+   const bool use_dev = b.UseDevice() || x.UseDevice();
+
+   auto bp = b.Read(use_dev);
+   auto xp = x.Write(use_dev);
+
+   auto Ap = Read(A, nnz);
+   auto Ip = Read(I, height+1);
+   auto Jp = Read(J, nnz);
+
    bool scale = (sc != 1.0);
-   for (int i = 0, j = 0; i < height; i++)
+   MFEM_FORALL(i, height,
    {
-      int end = I[i+1];
-      for ( ; true; j++)
+      int end = Ip[i+1];
+      for (int j = Ip[i]; true; j++)
       {
-         MFEM_VERIFY(j != end, "Couldn't find diagonal in row. i = " << i
-                     << ", j = " << j
-                     << ", I[i+1] = " << end );
-         if (J[j] == i)
+         if (j == end)
          {
-            MFEM_VERIFY(std::abs(A[j]) > 0.0, "Diagonal " << j << " must be nonzero");
+            //MFEM_ABORT_KERNEL("Diagonal not found in SparseMatrix::DiagScale");
+         }
+         if (Jp[j] == i)
+         {
+            if (!(std::abs(Ap[j]) > 0.0))
+            {
+               //MFEM_ABORT_KERNEL("Zero diagonal in SparseMatrix::DiagScale");
+            }
+
             if (scale)
             {
-               x(i) = sc * b(i) / A[j];
+               xp[i] = sc * bp[i] / Ap[j];
             }
             else
             {
-               x(i) = b(i) / A[j];
+               xp[i] = bp[i] / Ap[j];
             }
             break;
          }
       }
-      j = end;
-   }
+   });
    return;
 }
 
@@ -2749,6 +2866,10 @@ void SparseMatrix::Print(std::ostream & out, int _width) const
       return;
    }
 
+   // HostRead forces synchronization
+   HostReadI();
+   HostReadJ();
+   HostReadData();
    for (i = 0; i < height; i++)
    {
       out << "[row " << i << "]\n";
@@ -2938,6 +3059,16 @@ void SparseMatrix::Destroy()
    delete NodesMem;
 #endif
    delete At;
+
+#ifdef MFEM_USE_CUDA
+   if (initBuffers)
+   {
+      cusparseDestroySpMat(matA_descr);
+      cusparseDestroyDnVec(vecX_descr);
+      cusparseDestroyDnVec(vecY_descr);
+      initBuffers = false;
+   }
+#endif
 }
 
 int SparseMatrix::ActualWidth() const

linalg/sparsemat.hpp

@@ -21,6 +21,12 @@
 #include "../general/globals.hpp"
 #include "densemat.hpp"
 
+#ifdef MFEM_USE_CUDA
+#include <cusparse.h>
+#include <library_types.h>
+#include "../general/cuda.hpp"
+#endif
+
 namespace mfem
 {
 
@@ -80,9 +86,33 @@ protected:
    void Destroy();   // Delete all owned data
    void SetEmpty();  // Init all entries with empty values
 
+   bool useCuSparse{true}; //Use CuSparse if available
+
+   // Initialize CuSparse
+   void InitCuSparse();
+
+#ifdef MFEM_USE_CUDA
+   cusparseStatus_t status;
+   static cusparseHandle_t handle;
+   cusparseMatDescr_t descr=0;
+   static size_t bufferSize;
+   static void *dBuffer;
+   mutable bool initBuffers{false};
+
+   static int SparseMatrixCount;
+   mutable cusparseSpMatDescr_t matA_descr;
+   mutable cusparseDnVecDescr_t vecX_descr;
+   mutable cusparseDnVecDescr_t vecY_descr;
+#endif
+
 public:
    /// Create an empty SparseMatrix.
-   SparseMatrix() { SetEmpty(); }
+   SparseMatrix()
+   {
+      SetEmpty();
+
+      InitCuSparse();
+   }
 
    /** @brief Create a sparse matrix with flexible sparsity structure using a
        row-wise linked list (LIL) format. */
@@ -113,11 +143,14 @@ public:
    /** If @a mat is finalized and @a copy_graph is false, the #I and #J arrays
        will use a shallow copy (copy the pointers only) without transferring
        ownership. */
-   SparseMatrix(const SparseMatrix &mat, bool copy_graph = true);
+   SparseMatrix(const SparseMatrix &mat, bool copy_graph = true, MemoryType mt = MemoryType::SIZE);
 
    /// Create a SparseMatrix with diagonal @a v, i.e. A = Diag(v)
    SparseMatrix(const Vector & v);
 
+   // Runtime option to use CuSparse
+   // Only valid when using a CUDA backend
+   void UseCuSparse(bool _useCuSparse = true) { useCuSparse = _useCuSparse;}
 
    /// Assignment operator: deep copy
    SparseMatrix& operator=(const SparseMatrix &rhs);
@@ -573,11 +606,28 @@ public:
    void Swap(SparseMatrix &other);
 
    /// Destroys sparse matrix.
-   virtual ~SparseMatrix() { Destroy(); }
+   virtual ~SparseMatrix()
+   {
+      Destroy();
+#ifdef MFEM_USE_CUDA
+      if (handle && SparseMatrixCount==1 && Device::Allows(Backend::CUDA_MASK))
+      {
+         cusparseDestroy(handle);
+         CuMemFree(dBuffer);
+      }
+      SparseMatrixCount--;
+#endif
+   }
 
    Type GetType() const { return MFEM_SPARSEMAT; }
 };
 
+inline std::ostream& operator<<(std::ostream& os, SparseMatrix const& mat)
+{
+   mat.Print(os);
+   return os;
+}
+
 /// Applies f() to each element of the matrix (after it is finalized).
 void SparseMatrixFunction(SparseMatrix &S, double (*f)(double));
 

mesh/mesh.cpp

@@ -10524,17 +10524,17 @@ GeometricFactors::GeometricFactors(const Mesh *mesh, const IntegrationRule &ir,
    unsigned eval_flags = 0;
    if (flags & GeometricFactors::COORDINATES)
    {
-      X.SetSize(vdim*NQ*NE);
+      X.SetSize(vdim*NQ*NE, Device::GetDeviceTempMemoryType());
       eval_flags |= QuadratureInterpolator::VALUES;
    }
    if (flags & GeometricFactors::JACOBIANS)
    {
-      J.SetSize(dim*vdim*NQ*NE);
+      J.SetSize(dim*vdim*NQ*NE, Device::GetDeviceTempMemoryType());
       eval_flags |= QuadratureInterpolator::DERIVATIVES;
    }
    if (flags & GeometricFactors::DETERMINANTS)
    {
-      detJ.SetSize(NQ*NE);
+      detJ.SetSize(NQ*NE, Device::GetDeviceTempMemoryType());
       eval_flags |= QuadratureInterpolator::DETERMINANTS;
    }
 
@@ -10544,7 +10544,7 @@ GeometricFactors::GeometricFactors(const Mesh *mesh, const IntegrationRule &ir,
    qi->SetOutputLayout(QVectorLayout::byNODES);
    if (elem_restr)
    {
-      Vector Enodes(vdim*ND*NE);
+      Vector Enodes(vdim*ND*NE, Device::GetDeviceTempMemoryType());
       elem_restr->Mult(*nodes, Enodes);
       qi->Mult(Enodes, eval_flags, X, J, detJ);
    }

mesh/pmesh.cpp

@@ -2393,9 +2393,18 @@ void ParMesh::GetGhostFaceTransformation(
 }
 
 FaceElementTransformations *ParMesh::
-GetSharedFaceTransformations(int sf, bool fill2)
+GetSharedFaceTransformations(int sf, bool fill2, bool direct)
 {
-   int FaceNo = GetSharedFace(sf);
+   //int FaceNo = GetSharedFace(sf);
+   int FaceNo;
+   if (direct)
+   {
+      FaceNo = sf;
+   }
+   else
+   {
+      FaceNo = GetSharedFace(sf);
+   }
 
    FaceInfo &face_info = faces_info[FaceNo];
 

mesh/pmesh.hpp

@@ -291,9 +291,9 @@ public:
 
    /** Get the FaceElementTransformations for the given shared face (edge 2D).
        In the returned object, 1 and 2 refer to the local and the neighbor
-       elements, respectively. */
+       elements, respectively. Use direct if sf is the face number */
    FaceElementTransformations *
-   GetSharedFaceTransformations(int sf, bool fill2 = true);
+   GetSharedFaceTransformations(int sf, bool fill2 = true, bool direct = false);
 
    /// Return the number of shared faces (3D), edges (2D), vertices (1D)
    int GetNSharedFaces() const;