more template instantiations

add more specializations

fem/bilininteg.hpp

@@ -2174,7 +2174,7 @@ public:
 /** Class for local mass matrix assembling a(u,v) := (Q u, v) */
 class MassIntegrator: public BilinearFormIntegrator
 {
-protected:
+public:
 #ifndef MFEM_THREAD_SAFE
    Vector shape, te_shape;
 #endif
@@ -2941,6 +2941,7 @@ protected:
    const DofToQuad *maps;             ///< Not owned
    const FaceGeometricFactors *geom;  ///< Not owned
    int dim, nf, nq, dofs1D, quad1D;
+   Vector *external_vel = nullptr;
 
 private:
    Vector shape1, shape2;
@@ -2951,12 +2952,13 @@ public:
    { rho = NULL; u = &u_; alpha = a; beta = 0.5*a; }
 
    /// Construct integrator with rho = 1.
-   DGTraceIntegrator(VectorCoefficient &u_, double a, double b)
-   { rho = NULL; u = &u_; alpha = a; beta = b; }
+   DGTraceIntegrator(VectorCoefficient &u_, double a, double b,
+                     Vector *external_vel_ = nullptr)
+   { rho = NULL; u = &u_; alpha = a; beta = b; external_vel = external_vel_;}
 
    DGTraceIntegrator(Coefficient &rho_, VectorCoefficient &u_,
                      double a, double b)
-   { rho = &rho_; u = &u_; alpha = a; beta = b; }
+   { rho = &rho_; u = &u_; alpha = a; beta = b;}
 
    using BilinearFormIntegrator::AssembleFaceMatrix;
    virtual void AssembleFaceMatrix(const FiniteElement &el1,

fem/bilininteg_dgtrace_pa.cpp

@@ -159,6 +159,7 @@ void DGTraceIntegrator::SetupPA(const FiniteElementSpace &fes, FaceType type)
    quad1D = maps->nqpt;
    pa_data.SetSize(symmDims * nq * nf, Device::GetMemoryType());
    Vector vel;
+
    if (VectorConstantCoefficient *c_u = dynamic_cast<VectorConstantCoefficient*>
                                         (u))
    {
@@ -293,9 +294,21 @@ void DGTraceIntegrator::SetupPA(const FiniteElementSpace &fes, FaceType type)
       }
       MFEM_VERIFY(f_ind==nf, "Incorrect number of faces.");
    }
-   PADGTraceSetup(dim, dofs1D, quad1D, nf, ir->GetWeights(),
-                  geom->detJ, geom->normal, r, vel,
-                  alpha, beta, pa_data);
+
+   if (external_vel)
+   {
+
+      PADGTraceSetup(dim, dofs1D, quad1D, nf, ir->GetWeights(),
+                     geom->detJ, geom->normal, r, *external_vel,
+                     alpha, beta, pa_data);
+   }
+   else
+   {
+
+      PADGTraceSetup(dim, dofs1D, quad1D, nf, ir->GetWeights(),
+                     geom->detJ, geom->normal, r, vel,
+                     alpha, beta, pa_data);
+   }
 }
 
 void DGTraceIntegrator::AssemblePAInteriorFaces(const FiniteElementSpace& fes)
@@ -676,13 +689,16 @@ static void PADGTraceApply(const int dim,
                            const Vector &x,
                            Vector &y)
 {
+
    if (dim == 2)
    {
       switch ((D1D << 4 ) | Q1D)
       {
          case 0x22: return PADGTraceApply2D<2,2>(NF,B,Bt,op,x,y);
          case 0x33: return PADGTraceApply2D<3,3>(NF,B,Bt,op,x,y);
+         case 0x34: return PADGTraceApply2D<3,4>(NF,B,Bt,op,x,y);
          case 0x44: return PADGTraceApply2D<4,4>(NF,B,Bt,op,x,y);
+         case 0x46: return PADGTraceApply2D<4,6>(NF,B,Bt,op,x,y);
          case 0x55: return PADGTraceApply2D<5,5>(NF,B,Bt,op,x,y);
          case 0x66: return PADGTraceApply2D<6,6>(NF,B,Bt,op,x,y);
          case 0x77: return PADGTraceApply2D<7,7>(NF,B,Bt,op,x,y);
@@ -696,8 +712,8 @@ static void PADGTraceApply(const int dim,
       switch ((D1D << 4 ) | Q1D)
       {
          case 0x23: return SmemPADGTraceApply3D<2,3,1>(NF,B,Bt,op,x,y);
-         case 0x34: return SmemPADGTraceApply3D<3,4,2>(NF,B,Bt,op,x,y);
-         case 0x45: return SmemPADGTraceApply3D<4,5,2>(NF,B,Bt,op,x,y);
+         case 0x34: return SmemPADGTraceApply3D<3,4,1>(NF,B,Bt,op,x,y);
+         case 0x45: return SmemPADGTraceApply3D<4,5,1>(NF,B,Bt,op,x,y);
          case 0x56: return SmemPADGTraceApply3D<5,6,1>(NF,B,Bt,op,x,y);
          case 0x67: return SmemPADGTraceApply3D<6,7,1>(NF,B,Bt,op,x,y);
          case 0x78: return SmemPADGTraceApply3D<7,8,1>(NF,B,Bt,op,x,y);
@@ -1111,7 +1127,9 @@ static void PADGTraceApplyTranspose(const int dim,
       {
          case 0x22: return PADGTraceApplyTranspose2D<2,2>(NF,B,Bt,op,x,y);
          case 0x33: return PADGTraceApplyTranspose2D<3,3>(NF,B,Bt,op,x,y);
+         case 0x34: return PADGTraceApplyTranspose2D<3,4>(NF,B,Bt,op,x,y);
          case 0x44: return PADGTraceApplyTranspose2D<4,4>(NF,B,Bt,op,x,y);
+         case 0x46: return PADGTraceApplyTranspose2D<4,6>(NF,B,Bt,op,x,y);
          case 0x55: return PADGTraceApplyTranspose2D<5,5>(NF,B,Bt,op,x,y);
          case 0x66: return PADGTraceApplyTranspose2D<6,6>(NF,B,Bt,op,x,y);
          case 0x77: return PADGTraceApplyTranspose2D<7,7>(NF,B,Bt,op,x,y);

fem/bilininteg_mass_pa.hpp

@@ -0,0 +1,627 @@
+// Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef MFEM_BILININTEG_MASS_PA_HPP
+#define MFEM_BILININTEG_MASS_PA_HPP
+
+#include "../config/config.hpp"
+#include "../general/forall.hpp"
+#include "../linalg/dtensor.hpp"
+
+namespace mfem
+{
+
+template <bool ACCUMULATE = true>
+MFEM_HOST_DEVICE inline
+void PAMassApply2D_Element(const int e,
+                           const int NE,
+                           const double *b_,
+                           const double *bt_,
+                           const double *d_,
+                           const double *x_,
+                           double *y_,
+                           const int d1d = 0,
+                           const int q1d = 0)
+{
+   const int D1D = d1d;
+   const int Q1D = q1d;
+   auto B = ConstDeviceMatrix(b_, Q1D, D1D);
+   auto Bt = ConstDeviceMatrix(bt_, D1D, Q1D);
+   auto D = ConstDeviceCube(d_, Q1D, Q1D, NE);
+   auto X = ConstDeviceCube(x_, D1D, D1D, NE);
+   auto Y = DeviceCube(y_, D1D, D1D, NE);
+
+   if (!ACCUMULATE)
+   {
+      for (int dy = 0; dy < D1D; ++dy)
+      {
+         for (int dx = 0; dx < D1D; ++dx)
+         {
+            Y(dx, dy, e) = 0.0;
+         }
+      }
+   }
+
+   constexpr int max_D1D = MAX_D1D;
+   constexpr int max_Q1D = MAX_Q1D;
+   double sol_xy[max_Q1D][max_Q1D];
+   for (int qy = 0; qy < Q1D; ++qy)
+   {
+      for (int qx = 0; qx < Q1D; ++qx)
+      {
+         sol_xy[qy][qx] = 0.0;
+      }
+   }
+   for (int dy = 0; dy < D1D; ++dy)
+   {
+      double sol_x[max_Q1D];
+      for (int qy = 0; qy < Q1D; ++qy)
+      {
+         sol_x[qy] = 0.0;
+      }
+      for (int dx = 0; dx < D1D; ++dx)
+      {
+         const double s = X(dx,dy,e);
+         for (int qx = 0; qx < Q1D; ++qx)
+         {
+            sol_x[qx] += B(qx,dx)* s;
+         }
+      }
+      for (int qy = 0; qy < Q1D; ++qy)
+      {
+         const double d2q = B(qy,dy);
+         for (int qx = 0; qx < Q1D; ++qx)
+         {
+            sol_xy[qy][qx] += d2q * sol_x[qx];
+         }
+      }
+   }
+   for (int qy = 0; qy < Q1D; ++qy)
+   {
+      for (int qx = 0; qx < Q1D; ++qx)
+      {
+         sol_xy[qy][qx] *= D(qx,qy,e);
+      }
+   }
+   for (int qy = 0; qy < Q1D; ++qy)
+   {
+      double sol_x[max_D1D];
+      for (int dx = 0; dx < D1D; ++dx)
+      {
+         sol_x[dx] = 0.0;
+      }
+      for (int qx = 0; qx < Q1D; ++qx)
+      {
+         const double s = sol_xy[qy][qx];
+         for (int dx = 0; dx < D1D; ++dx)
+         {
+            sol_x[dx] += Bt(dx,qx) * s;
+         }
+      }
+      for (int dy = 0; dy < D1D; ++dy)
+      {
+         const double q2d = Bt(dy,qy);
+         for (int dx = 0; dx < D1D; ++dx)
+         {
+            Y(dx,dy,e) += q2d * sol_x[dx];
+         }
+      }
+   }
+}
+
+template<int T_D1D, int T_Q1D, int T_NBZ, bool ACCUMULATE = true>
+MFEM_HOST_DEVICE inline
+void SmemPAMassApply2D_Element(const int e,
+                               const int NE,
+                               const double *b_,
+                               const double *d_,
+                               const double *x_,
+                               double *y_,
+                               int d1d = 0,
+                               int q1d = 0)
+{
+   const int D1D = T_D1D ? T_D1D : d1d;
+   const int Q1D = T_Q1D ? T_Q1D : q1d;
+   constexpr int NBZ = T_NBZ ? T_NBZ : 1;
+
+   constexpr int MQ1 = T_Q1D ? T_Q1D : MAX_Q1D;
+   constexpr int MD1 = T_D1D ? T_D1D : MAX_D1D;
+   constexpr int MDQ = (MQ1 > MD1) ? MQ1 : MD1;
+
+   auto b = ConstDeviceMatrix(b_, Q1D, D1D);
+   auto D = ConstDeviceCube(d_, Q1D, Q1D, NE);
+   auto x = ConstDeviceCube(x_, D1D, D1D, NE);
+   auto Y = DeviceCube(y_, D1D, D1D, NE);
+
+   const int tidz = MFEM_THREAD_ID(z);
+
+   MFEM_SHARED double BBt[MQ1*MD1];
+   double (*B)[MD1] = (double (*)[MD1]) BBt;
+   double (*Bt)[MQ1] = (double (*)[MQ1]) BBt;
+   MFEM_SHARED double sm0[NBZ][MDQ*MDQ];
+   MFEM_SHARED double sm1[NBZ][MDQ*MDQ];
+   double (*X)[MD1] = (double (*)[MD1]) (sm0 + tidz);
+   double (*DQ)[MQ1] = (double (*)[MQ1]) (sm1 + tidz);
+   double (*QQ)[MQ1] = (double (*)[MQ1]) (sm0 + tidz);
+   double (*QD)[MD1] = (double (*)[MD1]) (sm1 + tidz);
+
+
+   MFEM_FOREACH_THREAD(dy,y,D1D)
+   {
+      MFEM_FOREACH_THREAD(dx,x,D1D)
+      {
+         X[dy][dx] = x(dx,dy,e);
+      }
+   }
+   if (tidz == 0)
+   {
+      MFEM_FOREACH_THREAD(dy,y,D1D)
+      {
+         MFEM_FOREACH_THREAD(q,x,Q1D)
+         {
+            B[q][dy] = b(q,dy);
+         }
+      }
+   }
+   MFEM_SYNC_THREAD;
+   MFEM_FOREACH_THREAD(dy,y,D1D)
+   {
+      MFEM_FOREACH_THREAD(qx,x,Q1D)
+      {
+         double dq = 0.0;
+         for (int dx = 0; dx < D1D; ++dx)
+         {
+            dq += X[dy][dx] * B[qx][dx];
+         }
+         DQ[dy][qx] = dq;
+      }
+   }
+   MFEM_SYNC_THREAD;
+   MFEM_FOREACH_THREAD(qy,y,Q1D)
+   {
+      MFEM_FOREACH_THREAD(qx,x,Q1D)
+      {
+         double qq = 0.0;
+         for (int dy = 0; dy < D1D; ++dy)
+         {
+            qq += DQ[dy][qx] * B[qy][dy];
+         }
+         QQ[qy][qx] = qq * D(qx, qy, e);
+      }
+   }
+   MFEM_SYNC_THREAD;
+   if (tidz == 0)
+   {
+      MFEM_FOREACH_THREAD(dy,y,D1D)
+      {
+         MFEM_FOREACH_THREAD(q,x,Q1D)
+         {
+            Bt[dy][q] = b(q,dy);
+         }
+      }
+   }
+   MFEM_SYNC_THREAD;
+   MFEM_FOREACH_THREAD(qy,y,Q1D)
+   {
+      MFEM_FOREACH_THREAD(dx,x,D1D)
+      {
+         double dq = 0.0;
+         for (int qx = 0; qx < Q1D; ++qx)
+         {
+            dq += QQ[qy][qx] * Bt[dx][qx];
+         }
+         QD[qy][dx] = dq;
+      }
+   }
+   MFEM_SYNC_THREAD;
+   MFEM_FOREACH_THREAD(dy,y,D1D)
+   {
+      MFEM_FOREACH_THREAD(dx,x,D1D)
+      {
+         double dd = 0.0;
+         for (int qy = 0; qy < Q1D; ++qy)
+         {
+            dd += (QD[qy][dx] * Bt[dy][qy]);
+         }
+         if (ACCUMULATE)
+         {
+            Y(dx, dy, e) += dd;
+         }
+         else
+         {
+            Y(dx, dy, e) = dd;
+         }
+      }
+   }
+}
+
+template <bool ACCUMULATE = true>
+MFEM_HOST_DEVICE inline
+void PAMassApply3D_Element(const int e,
+                           const int NE,
+                           const double *b_,
+                           const double *bt_,
+                           const double *d_,
+                           const double *x_,
+                           double *y_,
+                           const int d1d,
+                           const int q1d)
+{
+   const int D1D = d1d;
+   const int Q1D = q1d;
+   auto B = ConstDeviceMatrix(b_, Q1D, D1D);
+   auto Bt = ConstDeviceMatrix(bt_, D1D, Q1D);
+   auto D = DeviceTensor<4,const double>(d_, Q1D, Q1D, Q1D, NE);
+   auto X = DeviceTensor<4,const double>(x_, D1D, D1D, D1D, NE);
+   auto Y = DeviceTensor<4,double>(y_, D1D, D1D, D1D, NE);
+
+   if (!ACCUMULATE)
+   {
+      for (int dz = 0; dz < D1D; ++dz)
+      {
+         for (int dy = 0; dy < D1D; ++dy)
+         {
+            for (int dx = 0; dx < D1D; ++dx)
+            {
+               Y(dx, dy, dz, e) = 0.0;
+            }
+         }
+      }
+   }
+
+   constexpr int max_D1D = MAX_D1D;
+   constexpr int max_Q1D = MAX_Q1D;
+   double sol_xyz[max_Q1D][max_Q1D][max_Q1D];
+   for (int qz = 0; qz < Q1D; ++qz)
+   {
+      for (int qy = 0; qy < Q1D; ++qy)
+      {
+         for (int qx = 0; qx < Q1D; ++qx)
+         {
+            sol_xyz[qz][qy][qx] = 0.0;
+         }
+      }
+   }
+   for (int dz = 0; dz < D1D; ++dz)
+   {
+      double sol_xy[max_Q1D][max_Q1D];
+      for (int qy = 0; qy < Q1D; ++qy)
+      {
+         for (int qx = 0; qx < Q1D; ++qx)
+         {
+            sol_xy[qy][qx] = 0.0;
+         }
+      }
+      for (int dy = 0; dy < D1D; ++dy)
+      {
+         double sol_x[max_Q1D];
+         for (int qx = 0; qx < Q1D; ++qx)
+         {
+            sol_x[qx] = 0;
+         }
+         for (int dx = 0; dx < D1D; ++dx)
+         {
+            const double s = X(dx,dy,dz,e);
+            for (int qx = 0; qx < Q1D; ++qx)
+            {
+               sol_x[qx] += B(qx,dx) * s;
+            }
+         }
+         for (int qy = 0; qy < Q1D; ++qy)
+         {
+            const double wy = B(qy,dy);
+            for (int qx = 0; qx < Q1D; ++qx)
+            {
+               sol_xy[qy][qx] += wy * sol_x[qx];
+            }
+         }
+      }
+      for (int qz = 0; qz < Q1D; ++qz)
+      {
+         const double wz = B(qz,dz);
+         for (int qy = 0; qy < Q1D; ++qy)
+         {
+            for (int qx = 0; qx < Q1D; ++qx)
+            {
+               sol_xyz[qz][qy][qx] += wz * sol_xy[qy][qx];
+            }
+         }
+      }
+   }
+   for (int qz = 0; qz < Q1D; ++qz)
+   {
+      for (int qy = 0; qy < Q1D; ++qy)
+      {
+         for (int qx = 0; qx < Q1D; ++qx)
+         {
+            sol_xyz[qz][qy][qx] *= D(qx,qy,qz,e);
+         }
+      }
+   }
+   for (int qz = 0; qz < Q1D; ++qz)
+   {
+      double sol_xy[max_D1D][max_D1D];
+      for (int dy = 0; dy < D1D; ++dy)
+      {
+         for (int dx = 0; dx < D1D; ++dx)
+         {
+            sol_xy[dy][dx] = 0;
+         }
+      }
+      for (int qy = 0; qy < Q1D; ++qy)
+      {
+         double sol_x[max_D1D];
+         for (int dx = 0; dx < D1D; ++dx)
+         {
+            sol_x[dx] = 0;
+         }
+         for (int qx = 0; qx < Q1D; ++qx)
+         {
+            const double s = sol_xyz[qz][qy][qx];
+            for (int dx = 0; dx < D1D; ++dx)
+            {
+               sol_x[dx] += Bt(dx,qx) * s;
+            }
+         }
+         for (int dy = 0; dy < D1D; ++dy)
+         {
+            const double wy = Bt(dy,qy);
+            for (int dx = 0; dx < D1D; ++dx)
+            {
+               sol_xy[dy][dx] += wy * sol_x[dx];
+            }
+         }
+      }
+      for (int dz = 0; dz < D1D; ++dz)
+      {
+         const double wz = Bt(dz,qz);
+         for (int dy = 0; dy < D1D; ++dy)
+         {
+            for (int dx = 0; dx < D1D; ++dx)
+            {
+               Y(dx,dy,dz,e) += wz * sol_xy[dy][dx];
+            }
+         }
+      }
+   }
+}
+
+template<int T_D1D, int T_Q1D, bool ACCUMULATE = true>
+MFEM_HOST_DEVICE inline
+void SmemPAMassApply3D_Element(const int e,
+                               const int NE,
+                               const double *b_,
+                               const double *d_,
+                               const double *x_,
+                               double *y_,
+                               const int d1d = 0,
+                               const int q1d = 0)
+{
+   constexpr int D1D = T_D1D ? T_D1D : d1d;
+   constexpr int Q1D = T_Q1D ? T_Q1D : q1d;
+   constexpr int MQ1 = T_Q1D ? T_Q1D : MAX_Q1D;
+   constexpr int MD1 = T_D1D ? T_D1D : MAX_D1D;
+   constexpr int MDQ = (MQ1 > MD1) ? MQ1 : MD1;
+
+   auto b = ConstDeviceMatrix(b_, Q1D, D1D);
+   auto d = DeviceTensor<4,const double>(d_, Q1D, Q1D, Q1D, NE);
+   auto x = DeviceTensor<4,const double>(x_, D1D, D1D, D1D, NE);
+   auto y = DeviceTensor<4,double>(y_, D1D, D1D, D1D, NE);
+
+   MFEM_SHARED double sDQ[MQ1*MD1];
+   double (*B)[MD1] = (double (*)[MD1]) sDQ;
+   double (*Bt)[MQ1] = (double (*)[MQ1]) sDQ;
+   MFEM_SHARED double sm0[MDQ*MDQ*MDQ];
+   MFEM_SHARED double sm1[MDQ*MDQ*MDQ];
+   double (*X)[MD1][MD1]   = (double (*)[MD1][MD1]) sm0;
+   double (*DDQ)[MD1][MQ1] = (double (*)[MD1][MQ1]) sm1;
+   double (*DQQ)[MQ1][MQ1] = (double (*)[MQ1][MQ1]) sm0;
+   double (*QQQ)[MQ1][MQ1] = (double (*)[MQ1][MQ1]) sm1;
+   double (*QQD)[MQ1][MD1] = (double (*)[MQ1][MD1]) sm0;
+   double (*QDD)[MD1][MD1] = (double (*)[MD1][MD1]) sm1;
+   MFEM_FOREACH_THREAD(dy,y,D1D)
+   {
+      MFEM_FOREACH_THREAD(dx,x,D1D)
+      {
+         MFEM_UNROLL(MD1)
+         for (int dz = 0; dz < D1D; ++dz)
+         {
+            X[dz][dy][dx] = x(dx,dy,dz,e);
+         }
+      }
+      MFEM_FOREACH_THREAD(dx,x,Q1D)
+      {
+         B[dx][dy] = b(dx,dy);
+      }
+   }
+   MFEM_SYNC_THREAD;
+   MFEM_FOREACH_THREAD(dy,y,D1D)
+   {
+      MFEM_FOREACH_THREAD(qx,x,Q1D)
+      {
+         double u[D1D];
+         MFEM_UNROLL(MD1)
+         for (int dz = 0; dz < D1D; dz++)
+         {
+            u[dz] = 0;
+         }
+         MFEM_UNROLL(MD1)
+         for (int dx = 0; dx < D1D; ++dx)
+         {
+            MFEM_UNROLL(MD1)
+            for (int dz = 0; dz < D1D; ++dz)
+            {
+               u[dz] += X[dz][dy][dx] * B[qx][dx];
+            }
+         }
+         MFEM_UNROLL(MD1)
+         for (int dz = 0; dz < D1D; ++dz)
+         {
+            DDQ[dz][dy][qx] = u[dz];
+         }
+      }
+   }
+   MFEM_SYNC_THREAD;
+   MFEM_FOREACH_THREAD(qy,y,Q1D)
+   {
+      MFEM_FOREACH_THREAD(qx,x,Q1D)
+      {
+         double u[D1D];
+         MFEM_UNROLL(MD1)
+         for (int dz = 0; dz < D1D; dz++)
+         {
+            u[dz] = 0;
+         }
+         MFEM_UNROLL(MD1)
+         for (int dy = 0; dy < D1D; ++dy)
+         {
+            MFEM_UNROLL(MD1)
+            for (int dz = 0; dz < D1D; dz++)
+            {
+               u[dz] += DDQ[dz][dy][qx] * B[qy][dy];
+            }
+         }
+         MFEM_UNROLL(MD1)
+         for (int dz = 0; dz < D1D; dz++)
+         {
+            DQQ[dz][qy][qx] = u[dz];
+         }
+      }
+   }
+   MFEM_SYNC_THREAD;
+   MFEM_FOREACH_THREAD(qy,y,Q1D)
+   {
+      MFEM_FOREACH_THREAD(qx,x,Q1D)
+      {
+         double u[Q1D];
+         MFEM_UNROLL(MQ1)
+         for (int qz = 0; qz < Q1D; qz++)
+         {
+            u[qz] = 0;
+         }
+         MFEM_UNROLL(MD1)
+         for (int dz = 0; dz < D1D; ++dz)
+         {
+            MFEM_UNROLL(MQ1)
+            for (int qz = 0; qz < Q1D; qz++)
+            {
+               u[qz] += DQQ[dz][qy][qx] * B[qz][dz];
+            }
+         }
+         MFEM_UNROLL(MQ1)
+         for (int qz = 0; qz < Q1D; qz++)
+         {
+            QQQ[qz][qy][qx] = u[qz] * d(qx,qy,qz,e);
+         }
+      }
+   }
+   MFEM_SYNC_THREAD;
+   MFEM_FOREACH_THREAD(di,y,D1D)
+   {
+      MFEM_FOREACH_THREAD(q,x,Q1D)
+      {
+         Bt[di][q] = b(q,di);
+      }
+   }
+   MFEM_SYNC_THREAD;
+   MFEM_FOREACH_THREAD(qy,y,Q1D)
+   {
+      MFEM_FOREACH_THREAD(dx,x,D1D)
+      {
+         double u[Q1D];
+         MFEM_UNROLL(MQ1)
+         for (int qz = 0; qz < Q1D; ++qz)
+         {
+            u[qz] = 0;
+         }
+         MFEM_UNROLL(MQ1)
+         for (int qx = 0; qx < Q1D; ++qx)
+         {
+            MFEM_UNROLL(MQ1)
+            for (int qz = 0; qz < Q1D; ++qz)
+            {
+               u[qz] += QQQ[qz][qy][qx] * Bt[dx][qx];
+            }
+         }
+         MFEM_UNROLL(MQ1)
+         for (int qz = 0; qz < Q1D; ++qz)
+         {
+            QQD[qz][qy][dx] = u[qz];
+         }
+      }
+   }
+   MFEM_SYNC_THREAD;
+   MFEM_FOREACH_THREAD(dy,y,D1D)
+   {
+      MFEM_FOREACH_THREAD(dx,x,D1D)
+      {
+         double u[Q1D];
+         MFEM_UNROLL(MQ1)
+         for (int qz = 0; qz < Q1D; ++qz)
+         {
+            u[qz] = 0;
+         }
+         MFEM_UNROLL(MQ1)
+         for (int qy = 0; qy < Q1D; ++qy)
+         {
+            MFEM_UNROLL(MQ1)
+            for (int qz = 0; qz < Q1D; ++qz)
+            {
+               u[qz] += QQD[qz][qy][dx] * Bt[dy][qy];
+            }
+         }
+         MFEM_UNROLL(MQ1)
+         for (int qz = 0; qz < Q1D; ++qz)
+         {
+            QDD[qz][dy][dx] = u[qz];
+         }
+      }
+   }
+   MFEM_SYNC_THREAD;
+   MFEM_FOREACH_THREAD(dy,y,D1D)
+   {
+      MFEM_FOREACH_THREAD(dx,x,D1D)
+      {
+         double u[D1D];
+         MFEM_UNROLL(MD1)
+         for (int dz = 0; dz < D1D; ++dz)
+         {
+            u[dz] = 0;
+         }
+         MFEM_UNROLL(MQ1)
+         for (int qz = 0; qz < Q1D; ++qz)
+         {
+            MFEM_UNROLL(MD1)
+            for (int dz = 0; dz < D1D; ++dz)
+            {
+               u[dz] += QDD[qz][dy][dx] * Bt[dz][qz];
+            }
+         }
+         MFEM_UNROLL(MD1)
+         for (int dz = 0; dz < D1D; ++dz)
+         {
+            if (ACCUMULATE)
+            {
+               y(dx,dy,dz,e) += u[dz];
+            }
+            else
+            {
+               y(dx,dy,dz,e) = u[dz];
+            }
+         }
+      }
+   }
+   MFEM_SYNC_THREAD;
+}
+
+} // namespace mfem
+
+#endif

fem/dgmassinv.cpp

@@ -0,0 +1,641 @@
+// Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "dgmassinv.hpp"
+#include "bilinearform.hpp"
+#include "dgmassinv_kernels.hpp"
+#include "../general/forall.hpp"
+
+#ifdef MFEM_USE_CUDA
+#include <cublas.h>
+#include <cusolverDn.h>
+#endif
+
+namespace mfem
+{
+
+void DGMassInverse::Setup()
+{
+   M->Assemble();
+
+   diag_inv.SetSize(height);
+   M->AssembleDiagonal(diag_inv);
+
+   MakeReciprocal(diag_inv.Size(), diag_inv.ReadWrite());
+}
+
+DGMassInverse::DGMassInverse(FiniteElementSpace &fes_orig, Coefficient *coeff,
+                             const IntegrationRule *ir,
+                             int btype)
+   : Solver(fes_orig.GetTrueVSize()),
+     fec(fes_orig.GetMaxElementOrder(),
+         fes_orig.GetMesh()->Dimension(),
+         btype,
+         fes_orig.GetFE(0)->GetMapType()),
+     fes(fes_orig.GetMesh(), &fec)
+{
+   MFEM_VERIFY(fes.IsDGSpace(), "Space must be DG.");
+   MFEM_VERIFY(!fes.IsVariableOrder(), "Variable orders not supported.");
+
+   const int btype_orig =
+      static_cast<const L2_FECollection*>(fes_orig.FEColl())->GetBasisType();
+
+   if (btype_orig == btype)
+   {
+      // No change of basis required
+      d2q = nullptr;
+   }
+   else
+   {
+      // original basis to solver basis
+      const auto mode = DofToQuad::TENSOR;
+      d2q = &fes_orig.GetFE(0)->GetDofToQuad(fes.GetFE(0)->GetNodes(), mode);
+
+      int n = d2q->ndof;
+      Array<double> B_inv = d2q->B; // deep copy
+      Array<int> ipiv(n);
+      // solver basis to original
+      LUFactors lu(B_inv.HostReadWrite(), ipiv.HostWrite());
+      lu.Factor(n);
+      B_.SetSize(n*n);
+      lu.GetInverseMatrix(n, B_.HostWrite());
+      Bt_.SetSize(n*n);
+      DenseMatrix B_matrix(B_.HostReadWrite(), n, n);
+      DenseMatrix Bt_matrix(Bt_.HostWrite(), n, n);
+      Bt_matrix.Transpose(B_matrix);
+   }
+
+   if (coeff) { m = new MassIntegrator(*coeff, ir); }
+   else { m = new MassIntegrator(ir); }
+
+   M = new BilinearForm(&fes);
+   M->AddDomainIntegrator(m);
+   M->UseExternalIntegrators();
+   M->SetAssemblyLevel(AssemblyLevel::PARTIAL);
+   Setup();
+
+   // Workspace vectors used for CG
+   r_.SetSize(height);
+   d_.SetSize(height);
+   z_.SetSize(height);
+   // Only need transformed RHS if basis is different
+   if (btype_orig != btype) { b2_.SetSize(height); }
+}
+
+DGMassInverse::DGMassInverse(FiniteElementSpace &fes_, Coefficient &coeff,
+                             int btype)
+   : DGMassInverse(fes_, &coeff, nullptr, btype) { }
+
+DGMassInverse::DGMassInverse(FiniteElementSpace &fes_, Coefficient &coeff,
+                             const IntegrationRule &ir, int btype)
+   : DGMassInverse(fes_, &coeff, &ir, btype) { }
+
+DGMassInverse::DGMassInverse(FiniteElementSpace &fes_,
+                             const IntegrationRule &ir, int btype)
+   : DGMassInverse(fes_, nullptr, &ir, btype) { }
+
+DGMassInverse::DGMassInverse(FiniteElementSpace &fes_, int btype)
+   : DGMassInverse(fes_, nullptr, nullptr, btype) { }
+
+void DGMassInverse::SetOperator(const Operator &op)
+{
+   MFEM_ABORT("SetOperator not supported with DGMassInverse.")
+}
+
+void DGMassInverse::SetRelTol(const double rel_tol_) { rel_tol = rel_tol_; }
+
+void DGMassInverse::SetAbsTol(const double abs_tol_) { abs_tol = abs_tol_; }
+
+void DGMassInverse::SetMaxIter(const double max_iter_) { max_iter = max_iter_; }
+
+DGMassInverse::~DGMassInverse()
+{
+   delete M;
+   delete m;
+}
+
+template<int DIM, int D1D, int Q1D>
+void DGMassInverse::DGMassCGIteration(const Vector &b_, Vector &u_) const
+{
+   MFEM_PERF_SCOPE("DG Mass Inverse");
+
+   const int NE = fes.GetNE();
+   const int d1d = m->dofs1D;
+   const int q1d = m->quad1D;
+
+   const int ND = pow(d1d, DIM);
+
+   const auto B = m->maps->B.Read();
+   const auto Bt = m->maps->Bt.Read();
+   const auto pa_data = m->pa_data.Read();
+   const auto dinv = diag_inv.Read();
+   auto r = r_.Write();
+   auto d = d_.Write();
+   auto z = z_.Write();
+   auto u = u_.ReadWrite();
+
+   const double RELTOL = rel_tol;
+   const double ABSTOL = abs_tol;
+   const double MAXIT = max_iter;
+
+   const bool it_mode = iterative_mode;
+   const bool change_basis = (d2q != nullptr);
+
+   // b is the right-hand side (if no change of basis, this just points to the
+   // incoming RHS vector, if we have to change basis, this points to the
+   // internal b2 vector where we put the transformed RHS)
+   const double *b;
+   // the following are non-null if we have to change basis
+   double *b2 = nullptr; // non-const access to b2
+   const double *b_orig = nullptr; // RHS vector in "original" basis
+   const double *d2q_B = nullptr; // matrix to transform initial guess
+   const double *q2d_B = nullptr; // matrix to transform solution
+   const double *q2d_Bt = nullptr; // matrix to transform RHS
+   if (change_basis)
+   {
+      d2q_B = d2q->B.Read();
+      q2d_B = B_.Read();
+      q2d_Bt = Bt_.Read();
+
+      b2 = b2_.Write();
+      b_orig = b_.Read();
+      b = b2;
+   }
+   else
+   {
+      b = b_.Read();
+   }
+
+   const int NB = Q1D ? Q1D : 1; // block size
+
+   // printf("  El.     It.    (Br,r)\n");
+   // printf("=============================\n");
+   MFEM_FORALL_2D(e, NE, NB, NB, 1,
+   {
+      // Perform change of basis if needed
+      if (change_basis)
+      {
+         // Transform RHS
+         DGMassBasis<DIM,D1D,MAX_D1D>(e, NE, q2d_Bt, b_orig, b2, d1d);
+         if (it_mode)
+         {
+            // Transform initial guess
+            // Double check that "in-place" eval is OK here
+            DGMassBasis<DIM,D1D,MAX_D1D>(e, NE, d2q_B, u, u, d1d);
+         }
+      }
+
+      const int tid = MFEM_THREAD_ID(x) + NB*MFEM_THREAD_ID(y);
+      // int final_iter;
+      // double final_norm;
+      // bool converged;
+
+      if (it_mode)
+      {
+         DGMassApply<DIM,D1D,Q1D>(e, NE, B, Bt, pa_data, u, r, d1d, q1d);
+         DGMassAxpy(e, NE, ND, 1.0, b, -1.0, r, r); // r = b - r
+      }
+      else
+      {
+         const int BX = MFEM_THREAD_SIZE(x);
+         const int BY = MFEM_THREAD_SIZE(y);
+         const int bxy = BX*BY;
+         const auto B = ConstDeviceMatrix(b, ND, NE);
+         auto U = DeviceMatrix(u, ND, NE);
+         auto R = DeviceMatrix(r, ND, NE);
+         for (int i = tid; i < ND; i += bxy)
+         {
+            U(i, e) = 0.0;
+            R(i, e) = B(i, e);
+         }
+         MFEM_SYNC_THREAD;
+      }
+
+      // TODO: get rid of extra memory usage for z
+      DGMassPreconditioner(e, NE, ND, dinv, r, z);
+      DGMassAxpy(e, NE, ND, 1.0, z, 0.0, z, d); // d = z
+
+      double nom = DGMassDot<NB>(e, NE, ND, d, r);
+      // MFEM_ASSERT(IsFinite(nom), "nom = " << nom);
+
+      if (nom < 0.0)
+      {
+         return; // Not positive definite...
+      }
+      double r0 = fmax(nom*RELTOL*RELTOL, ABSTOL*ABSTOL);
+      if (nom <= r0)
+      {
+         // converged = true;
+         // final_iter = 0;
+         // final_norm = sqrt(nom);
+         // return;
+         return;
+      }
+
+      DGMassApply<DIM,D1D,Q1D>(e, NE, B, Bt, pa_data, d, z, d1d, q1d);
+      double den = DGMassDot<NB>(e, NE, ND, z, d);
+      if (den <= 0.0)
+      {
+         const double d2 = DGMassDot<NB>(e, NE, ND, d, d);
+         if (d2 > 0.0 && tid == 0) { printf("Not positive definite.\n"); }
+         if (den == 0.0)
+         {
+            // converged = false;
+            // final_iter = 0;
+            // final_norm = sqrt(nom);
+            // return;
+            return;
+         }
+      }
+
+      // start iteration
+      int i = 1;
+      while (true)
+      {
+         const double alpha = nom/den;
+         DGMassAxpy(e, NE, ND, 1.0, u, alpha, d, u); // u = u + alpha*d
+         DGMassAxpy(e, NE, ND, 1.0, r, -alpha, z, r); // r = r - alpha*A*d
+
+         DGMassPreconditioner(e, NE, ND, dinv, r, z);
+
+         double betanom = DGMassDot<NB>(e, NE, ND, r, z);
+         if (betanom < 0.0)
+         {
+            if (tid == 0) { printf("Not positive definite.\n"); }
+            // converged = false;
+            // final_iter = i;
+            // return;
+            break;
+         }
+
+         // if (tid == 0) { printf(" %4d    %4d    %10.6e\n", e, i, betanom); }
+
+         if (betanom <= r0)
+         {
+            // converged = true;
+            // final_iter = i;
+            // return;
+            break;
+         }
+
+         if (++i > MAXIT) { break; }
+
+         const double beta = betanom/nom;
+         DGMassAxpy(e, NE, ND, 1.0, z, beta, d, d); // d = z + beta*d
+         DGMassApply<DIM,D1D,Q1D>(e, NE, B, Bt, pa_data, d, z, d1d, q1d); // z = A d
+         den = DGMassDot<NB>(e, NE, ND, d, z);
+         if (den <= 0.0)
+         {
+            const double d2 = DGMassDot<NB>(e, NE, ND, d, d);
+            if (d2 > 0.0 && tid == 0) { printf("Not positive definite.\n"); }
+            if (den == 0.0)
+            {
+               // final_iter = i;
+               // return;
+               break;
+            }
+         }
+         nom = betanom;
+      }
+
+      if (change_basis)
+      {
+         // Double check that "in-place" eval is OK here
+         DGMassBasis<DIM,D1D,MAX_D1D>(e, NE, q2d_B, u, u, d1d);
+      }
+   });
+}
+
+void DGMassInverse::Mult(const Vector &Mu, Vector &u) const
+{
+   // Dispatch to templated version based on dim, d1d, and q1d.
+   const int dim = fes.GetMesh()->Dimension();
+   const int d1d = m->dofs1D;
+   const int q1d = m->quad1D;
+
+   const int id = (d1d << 4) | q1d;
+
+   // printf("dim = %d id = 0x%x\n", dim, id);
+
+   if (dim == 2)
+   {
+      switch (id)
+      {
+         case 0x11: return DGMassCGIteration<2,1,1>(Mu, u);
+         case 0x22: return DGMassCGIteration<2,2,2>(Mu, u);
+         case 0x23: return DGMassCGIteration<2,2,3>(Mu, u);
+         case 0x33: return DGMassCGIteration<2,3,3>(Mu, u);
+         case 0x35: return DGMassCGIteration<2,3,5>(Mu, u);
+         case 0x44: return DGMassCGIteration<2,4,4>(Mu, u);
+         case 0x46: return DGMassCGIteration<2,4,6>(Mu, u);
+         case 0x55: return DGMassCGIteration<2,5,5>(Mu, u);
+         case 0x57: return DGMassCGIteration<2,5,7>(Mu, u);
+         case 0x66: return DGMassCGIteration<2,6,6>(Mu, u);
+         case 0x68: return DGMassCGIteration<2,6,8>(Mu, u);
+         default:
+            printf("dim = %d, id = 0x%x\n", dim, id);
+            MFEM_ABORT("Fallback");
+            return DGMassCGIteration<2>(Mu, u);
+      }
+   }
+   else if (dim == 3)
+   {
+      switch (id)
+      {
+         case 0x22: return DGMassCGIteration<3,2,2>(Mu, u);
+         case 0x23: return DGMassCGIteration<3,2,3>(Mu, u);
+         case 0x33: return DGMassCGIteration<3,3,3>(Mu, u);
+         case 0x34: return DGMassCGIteration<3,3,4>(Mu, u);
+         case 0x35: return DGMassCGIteration<3,3,5>(Mu, u);
+         case 0x44: return DGMassCGIteration<3,4,4>(Mu, u);
+         case 0x45: return DGMassCGIteration<3,4,5>(Mu, u);
+         case 0x46: return DGMassCGIteration<3,4,6>(Mu, u);
+         case 0x48: return DGMassCGIteration<3,4,8>(Mu, u);
+         case 0x55: return DGMassCGIteration<3,5,5>(Mu, u);
+         case 0x56: return DGMassCGIteration<3,5,6>(Mu, u);
+         case 0x57: return DGMassCGIteration<3,5,7>(Mu, u);
+         case 0x58: return DGMassCGIteration<3,5,8>(Mu, u);
+         case 0x66: return DGMassCGIteration<3,6,6>(Mu, u);
+         case 0x67: return DGMassCGIteration<3,6,7>(Mu, u);
+         default:
+            printf("dim = %d, id = 0x%x\n", dim, id);
+            MFEM_ABORT("Fallback");
+            return DGMassCGIteration<3>(Mu, u);
+      }
+   }
+}
+
+#ifdef MFEM_USE_CUDA
+class CuSolver
+{
+protected:
+   cusolverDnHandle_t handle = nullptr;
+   CuSolver()
+   {
+      cusolverStatus_t status = cusolverDnCreate(&handle);
+      MFEM_VERIFY(status == CUSOLVER_STATUS_SUCCESS,
+                  "Cannot initialize CuSolver.");
+   }
+   ~CuSolver()
+   {
+      cusolverDnDestroy(handle);
+   }
+   static CuSolver &Instance()
+   {
+      static CuSolver instance;
+      return instance;
+   }
+public:
+   static cusolverDnHandle_t Handle()
+   {
+      return Instance().handle;
+   }
+};
+
+class CuBLAS
+{
+protected:
+   cublasHandle_t handle = nullptr;
+   CuBLAS()
+   {
+      cublasStatus_t status = cublasCreate(&handle);
+      MFEM_VERIFY(status == CUBLAS_STATUS_SUCCESS, "Cannot initialize cuBLAS.");
+   }
+   ~CuBLAS()
+   {
+      cublasDestroy(handle);
+   }
+   static CuBLAS &Instance()
+   {
+      static CuBLAS instance;
+      return instance;
+   }
+public:
+   static cublasHandle_t Handle()
+   {
+      return Instance().handle;
+   }
+};
+#endif
+
+DGMassInverse_Direct::DGMassInverse_Direct(FiniteElementSpace &fes_,
+                                           BatchSolverMode mode_)
+   : DGMassInverse_Direct(fes_, nullptr, nullptr, mode_) { }
+
+DGMassInverse_Direct::DGMassInverse_Direct(FiniteElementSpace &fes_,
+                                           Coefficient &coeff,
+                                           BatchSolverMode mode_)
+   : DGMassInverse_Direct(fes_, &coeff, nullptr, mode_) { }
+
+DGMassInverse_Direct::DGMassInverse_Direct(FiniteElementSpace &fes_,
+                                           Coefficient &coeff,
+                                           const IntegrationRule &ir,
+                                           BatchSolverMode mode_)
+   : DGMassInverse_Direct(fes_, &coeff, &ir, mode_) { }
+
+DGMassInverse_Direct::DGMassInverse_Direct(FiniteElementSpace &fes_,
+                                           Coefficient *coeff,
+                                           const IntegrationRule *ir,
+                                           BatchSolverMode mode_)
+   : Solver(fes_.GetTrueVSize()), fes(fes_), mode(mode_)
+{
+   const int ne = fes.GetNE();
+   const int elem_dofs = fes.GetFE(0)->GetDof();
+
+   blocks.SetSize(ne*elem_dofs*elem_dofs);
+
+   if (coeff) { m = new MassIntegrator(*coeff, ir); }
+   else { m = new MassIntegrator(ir); }
+
+   Setup();
+}
+
+void DGMassInverse_Direct::Setup()
+{
+   const int ne = fes.GetNE();
+   const int elem_dofs = fes.GetFE(0)->GetDof();
+   MFEM_CONTRACT_VAR(ne);
+   MFEM_CONTRACT_VAR(elem_dofs);
+
+   m->AssembleEA(fes, blocks, false);
+
+   tensor.UseExternalData(NULL, elem_dofs, elem_dofs, ne);
+   tensor.GetMemory().MakeAlias(blocks.GetMemory(), 0, blocks.Size());
+
+   if ((mode == BatchSolverMode::CUSOLVER || mode == BatchSolverMode::CUBLAS)
+       && !Device::Allows(Backend::CUDA))
+   {
+      MFEM_ABORT("Unsupported mode. Use BatchSolverMode::NATIVE.");
+   }
+
+   if (mode == BatchSolverMode::NATIVE || !Device::Allows(Backend::CUDA))
+   {
+      BatchLUFactor(tensor, ipiv);
+   }
+   else if (mode == BatchSolverMode::CUSOLVER)
+   {
+#ifdef MFEM_USE_CUDA
+      vector_array.SetSize(ne);
+      matrix_array.SetSize(ne);
+      double *ptr_base = blocks.ReadWrite();
+      for (int i = 0; i < ne; ++i)
+      {
+         matrix_array[i] = ptr_base + i*elem_dofs*elem_dofs;
+      }
+      info_array.SetSize(ne);
+
+      cusolverStatus_t status = cusolverDnDpotrfBatched(
+                                   CuSolver::Handle(),
+                                   CUBLAS_FILL_MODE_LOWER,
+                                   elem_dofs,
+                                   matrix_array.ReadWrite(),
+                                   elem_dofs,
+                                   info_array.Write(),
+                                   ne);
+      MFEM_VERIFY(status == CUSOLVER_STATUS_SUCCESS, "");
+#else
+      MFEM_ABORT("CUDA must be enabled.");
+#endif
+   }
+   else if (mode == BatchSolverMode::CUBLAS)
+   {
+#ifdef MFEM_USE_CUDA
+      ipiv.SetSize(ne*elem_dofs);
+      info_array.SetSize(ne);
+
+      Vector tmp(blocks);
+      Array<double*> tmp_array(ne);
+      double *ptr_base = tmp.ReadWrite();
+      for (int i = 0; i < ne; ++i)
+      {
+         tmp_array[i] = ptr_base + i*elem_dofs*elem_dofs;
+      }
+
+      matrix_array.SetSize(ne);
+      ptr_base = blocks.ReadWrite();
+      for (int i = 0; i < ne; ++i)
+      {
+         matrix_array[i] = ptr_base + i*elem_dofs*elem_dofs;
+      }
+
+      cublasStatus_t status = cublasDgetrfBatched(
+                                 CuBLAS::Handle(),
+                                 elem_dofs,
+                                 tmp_array.ReadWrite(),
+                                 elem_dofs,
+                                 ipiv.Write(),
+                                 info_array.Write(),
+                                 ne);
+      MFEM_VERIFY(status == CUBLAS_STATUS_SUCCESS, "");
+
+
+      status = cublasDgetriBatched(
+                  CuBLAS::Handle(),
+                  elem_dofs,
+                  tmp_array.ReadWrite(),
+                  elem_dofs,
+                  ipiv.ReadWrite(),
+                  matrix_array.ReadWrite(),
+                  elem_dofs,
+                  info_array.Write(),
+                  ne);
+      MFEM_VERIFY(status == CUBLAS_STATUS_SUCCESS, "");
+#else
+      MFEM_ABORT("CUDA must be enabled.");
+#endif
+   }
+}
+
+void DGMassInverse_Direct::Mult(const Vector &Mu, Vector &u) const
+{
+   MFEM_PERF_SCOPE("DG Mass Inverse");
+   if (mode == BatchSolverMode::CUBLAS)
+   {
+#ifdef MFEM_USE_CUDA
+      const int n = tensor.SizeI();
+      const int nblocks = tensor.SizeK();
+
+      const double alpha = 1.0, beta = 0.0;
+
+      // cublasStatus_t status = cublasDgemvStridedBatched(
+      //    CuBLAS::Handle(),
+      //    CUBLAS_OP_N,
+      //    &alpha,
+      //    blocks.Read(), n, n*n,
+      //    Mu.Read(), 1, n,
+      //    &beta,
+      //    u.Write(), 1, n,
+      //    nblocks);
+      cublasStatus_t status = cublasDgemmStridedBatched(
+                                 CuBLAS::Handle(),
+                                 CUBLAS_OP_N, CUBLAS_OP_N,
+                                 n, 1, n,
+                                 &alpha,
+                                 blocks.Read(), n, n*n,
+                                 Mu.Read(), n, n,
+                                 &beta,
+                                 u.Write(), n, n,
+                                 nblocks);
+#else
+      MFEM_ABORT("CUDA must be enabled.");
+#endif
+   }
+   else
+   {
+      u = Mu;
+      Solve(u);
+   }
+}
+
+void DGMassInverse_Direct::Solve(Vector &u) const
+{
+   if (mode == BatchSolverMode::NATIVE)
+   {
+      BatchLUSolve(tensor, ipiv, u);
+   }
+   else if (mode == BatchSolverMode::CUSOLVER)
+   {
+#ifdef MFEM_USE_CUDA
+      const int n = tensor.SizeI();
+      const int nblocks = tensor.SizeK();
+      double *ptr_base = u.ReadWrite();
+      auto u_ptr = Reshape(vector_array.Write(), nblocks);
+
+      MFEM_FORALL(i, nblocks, u_ptr[i] = ptr_base + i*n; );
+
+      cusolverStatus_t status = cusolverDnDpotrsBatched(
+                                   CuSolver::Handle(),
+                                   CUBLAS_FILL_MODE_LOWER,
+                                   n,
+                                   1,
+                                   matrix_array.ReadWrite(),
+                                   n,
+                                   u_ptr,
+                                   n,
+                                   info_array.Write(),
+                                   nblocks);
+#else
+      MFEM_ABORT("CUDA must be enabled.");
+#endif
+   }
+}
+
+void DGMassInverse_Direct::SetOperator(const Operator &op)
+{
+   MFEM_ABORT("Not supported.");
+}
+
+DGMassInverse_Direct::~DGMassInverse_Direct()
+{
+   delete m;
+}
+
+} // namespace mfem

fem/dgmassinv.hpp

@@ -0,0 +1,187 @@
+// Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef MFEM_DG_MASS_INVERSE
+#define MFEM_DG_MASS_INVERSE
+
+#include "../linalg/operator.hpp"
+#include "fespace.hpp"
+
+namespace mfem
+{
+
+/// @brief Solver for the discontinuous Galerkin mass matrix.
+///
+/// This class performs a @a local (diagonally preconditioned) conjugate
+/// gradient iteration for each element. Optionally, a change of basis is
+/// performed to iterate on a better-conditioned system.
+///
+/// This class fully supports execution on device.
+class DGMassInverse : public Solver
+{
+protected:
+   DG_FECollection fec; ///< FE collection in requested basis.
+   FiniteElementSpace fes; ///< FE space in requested basis.
+   const DofToQuad *d2q; ///< Change of basis. Not owned.
+   Array<double> B_; ///< Inverse of change of basis.
+   Array<double> Bt_; ///< Inverse of change of basis, transposed.
+   class BilinearForm *M;
+   class MassIntegrator *m; ///< Owned.
+   Vector diag_inv; ///< Jacobi preconditioner.
+   double rel_tol = 1e-12; ///< Relative CG tolerance.
+   double abs_tol = 1e-12; ///< Absolute CG tolerance.
+   int max_iter = 100; ///< Maximum number of CG iterations;
+
+   /// @name Intermediate vectors needed for CG three-term recurrence.
+   ///@{
+   mutable Vector r_, d_, z_, b2_;
+   ///@}
+
+   /// @brief Protected constructor, used internally.
+   ///
+   /// Custom coefficient and integration rule are used in @a coeff and @a ir
+   /// are non-NULL.
+   DGMassInverse(FiniteElementSpace &fes_, Coefficient *coeff,
+                 const IntegrationRule *ir, int btype);
+public:
+   /// @brief Construct the DG inverse mass operator for @a fes_.
+   ///
+   /// The basis type @a btype determines which basis should be used internally
+   /// in the solver. This <b>does not</b> have to be the same basis as @a fes_.
+   /// The best choice is typically BasisType::GaussLegendre because it is
+   /// well-preconditioned by its diagonal.
+   ///
+   /// The solution and right-hand side used for the solver are not affected by
+   /// this basis (they correspond to the basis of @a fes_). @a btype is only
+   /// used internally, and only has an effect on the convergence rate.
+   DGMassInverse(FiniteElementSpace &fes_, int btype=BasisType::GaussLegendre);
+   /// @brief Construct the DG inverse mass operator for @a fes_ with
+   /// Coefficient @a coeff.
+   ///
+   /// @sa DGMassInverse(FiniteElementSpace&, int) for information about @a
+   /// btype.
+   DGMassInverse(FiniteElementSpace &fes_, Coefficient &coeff,
+                 int btype=BasisType::GaussLegendre);
+   /// @brief Construct the DG inverse mass operator for @a fes_ with
+   /// Coefficient @a coeff and IntegrationRule @a ir.
+   ///
+   /// @sa DGMassInverse(FiniteElementSpace&, int) for information about @a
+   /// btype.
+   DGMassInverse(FiniteElementSpace &fes_, Coefficient &coeff,
+                 const IntegrationRule &ir, int btype=BasisType::GaussLegendre);
+   /// @brief Construct the DG inverse mass operator for @a fes_ with
+   /// IntegrationRule @a ir.
+   ///
+   /// @sa DGMassInverse(FiniteElementSpace&, int) for information about @a
+   /// btype.
+   DGMassInverse(FiniteElementSpace &fes_, const IntegrationRule &ir,
+                 int btype=BasisType::GaussLegendre);
+   /// Solve the system M b = u. Use @a u as an initial guess.
+   void Mult(const Vector &b, Vector &u) const;
+   /// Not implemented. Aborts.
+   void SetOperator(const Operator &op);
+   /// Set the relative tolerance.
+   void SetRelTol(const double rel_tol_);
+   /// Set the absolute tolerance.
+   void SetAbsTol(const double abs_tol_);
+   /// Set the maximum number of iterations.
+   void SetMaxIter(const double max_iter_);
+   /// Sets whether or not to use the second Mult argument as the initial guess.
+   void SetIterativeMode(const bool it_mode_) { iterative_mode = it_mode_; }
+
+   ~DGMassInverse();
+
+   void Setup();
+
+   // Not part of the public interface, must be public because it contains a
+   // kernel
+   template<int DIM, int D1D = 0, int Q1D = 0>
+   void DGMassCGIteration(const Vector &b_, Vector &u_) const;
+};
+
+enum class BatchSolverMode
+{
+   NATIVE,
+   CUSOLVER,
+   CUBLAS
+};
+
+class DGMassInverse_Direct : public Solver
+{
+protected:
+   FiniteElementSpace &fes;
+   const BatchSolverMode mode;
+   Vector blocks;
+   DenseTensor tensor;
+   Array<int> ipiv;
+   class MassIntegrator *m;
+
+   mutable Array<double*> matrix_array;
+   mutable Array<double*> vector_array;
+   mutable Array<int> info_array;
+
+   /// @brief Protected constructor, used internally.
+   ///
+   /// Custom coefficient and integration rule are used in @a coeff and @a ir
+   /// are non-NULL.
+   DGMassInverse_Direct(FiniteElementSpace &fes_, Coefficient *coeff,
+                        const IntegrationRule *ir,
+                        BatchSolverMode mode_ = BatchSolverMode::NATIVE);
+public:
+   /// @brief Construct the DG inverse mass operator for @a fes_.
+   DGMassInverse_Direct(FiniteElementSpace &fes_,
+                        BatchSolverMode mode_ = BatchSolverMode::NATIVE);
+
+   DGMassInverse_Direct(FiniteElementSpace &fes_,
+                        Coefficient &coeff_,
+                        BatchSolverMode mode_ = BatchSolverMode::NATIVE);
+
+   DGMassInverse_Direct(FiniteElementSpace &fes_,
+                        Coefficient &coeff_,
+                        const IntegrationRule &ir_,
+                        BatchSolverMode mode_ = BatchSolverMode::NATIVE);
+   // /// @brief Construct the DG inverse mass operator for @a fes_ with
+   // /// Coefficient @a coeff.
+   // ///
+   // /// @sa DGMassInverse(FiniteElementSpace&, int) for information about @a
+   // /// btype.
+   // DGMassInverse(FiniteElementSpace &fes_, Coefficient &coeff);
+   // /// @brief Construct the DG inverse mass operator for @a fes_ with
+   // /// Coefficient @a coeff and IntegrationRule @a ir.
+   // ///
+   // /// @sa DGMassInverse(FiniteElementSpace&, int) for information about @a
+   // /// btype.
+   // DGMassInverse(FiniteElementSpace &fes_, Coefficient &coeff,
+   //               const IntegrationRule &ir);
+   // /// @brief Construct the DG inverse mass operator for @a fes_ with
+   // /// IntegrationRule @a ir.
+   // ///
+   // /// @sa DGMassInverse(FiniteElementSpace&, int) for information about @a
+   // /// btype.
+   // DGMassInverse(FiniteElementSpace &fes_, const IntegrationRule &ir,
+   //               int btype=BasisType::GaussLegendre);
+
+   void Setup();
+
+   /// Solve the system M b = u.
+   void Mult(const Vector &b, Vector &u) const;
+
+   void Solve(Vector &u) const;
+
+   /// Not implemented. Aborts.
+   void SetOperator(const Operator &op);
+
+   ~DGMassInverse_Direct();
+};
+
+} // namespace mfem
+
+#endif

fem/dgmassinv_kernels.hpp

@@ -0,0 +1,292 @@
+// Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef MFEM_DGMASSINV_KERNELS_HPP
+#define MFEM_DGMASSINV_KERNELS_HPP
+
+#include "bilininteg.hpp"
+#include "../linalg/kernels.hpp"
+#include "../general/forall.hpp"
+#include "kernels.hpp"
+#include "bilininteg_mass_pa.hpp"
+
+namespace mfem
+{
+
+void MakeReciprocal(int n, double *x)
+{
+   MFEM_FORALL(i, n, x[i] = 1.0/x[i]; );
+}
+
+template <int DIM, int D1D, int Q1D>
+MFEM_HOST_DEVICE inline
+void DGMassApply(const int e,
+                 const int NE,
+                 const double *B,
+                 const double *Bt,
+                 const double *pa_data,
+                 const double *x,
+                 double *y,
+                 const int d1d = 0,
+                 const int q1d = 0)
+{
+   constexpr bool use_smem = (D1D > 0 && Q1D > 0);
+   constexpr bool ACCUM = false;
+   constexpr int NBZ = 1;
+   if (use_smem)
+   {
+      // cannot specialize functions below with D1D or Q1D equal to zero
+      // (this branch only runs with D1D and Q1D are both positive)
+      constexpr int TD1D = D1D ? D1D : 1;
+      constexpr int TQ1D = Q1D ? Q1D : 1;
+      if (DIM == 2)
+      {
+         SmemPAMassApply2D_Element<TD1D,TQ1D,NBZ,ACCUM>(e, NE, B, pa_data, x, y);
+      }
+      else if (DIM == 3)
+      {
+         SmemPAMassApply3D_Element<TD1D,TQ1D,ACCUM>(e, NE, B, pa_data, x, y);
+      }
+      else
+      {
+         MFEM_ABORT_KERNEL("Unsupported dimension.");
+      }
+   }
+   else
+   {
+      if (DIM == 2)
+      {
+         PAMassApply2D_Element<ACCUM>(e, NE, B, Bt, pa_data, x, y, d1d, q1d);
+      }
+      else if (DIM == 3)
+      {
+         PAMassApply3D_Element<ACCUM>(e, NE, B, Bt, pa_data, x, y, d1d, q1d);
+      }
+      else
+      {
+         MFEM_ABORT_KERNEL("Unsupported dimension.");
+      }
+   }
+}
+
+MFEM_HOST_DEVICE inline
+void DGMassPreconditioner(const int e,
+                          const int NE,
+                          const int ND,
+                          const double *dinv,
+                          const double *x,
+                          double *y)
+{
+   const auto X = ConstDeviceMatrix(x, ND, NE);
+   const auto D = ConstDeviceMatrix(dinv, ND, NE);
+   auto Y = DeviceMatrix(y, ND, NE);
+
+   const int tid = MFEM_THREAD_ID(x) + MFEM_THREAD_SIZE(x)*MFEM_THREAD_ID(y);
+   const int bxy = MFEM_THREAD_SIZE(x)*MFEM_THREAD_SIZE(y);
+
+   for (int i = tid; i < ND; i += bxy)
+   {
+      Y(i, e) = D(i, e)*X(i, e);
+   }
+   MFEM_SYNC_THREAD;
+}
+
+MFEM_HOST_DEVICE inline
+void DGMassAxpy(const int e,
+                const int NE,
+                const int ND,
+                const double a,
+                const double *x,
+                const double b,
+                const double *y,
+                double *z)
+{
+   const auto X = ConstDeviceMatrix(x, ND, NE);
+   const auto Y = ConstDeviceMatrix(y, ND, NE);
+   auto Z = DeviceMatrix(z, ND, NE);
+
+   const int tid = MFEM_THREAD_ID(x) + MFEM_THREAD_SIZE(x)*MFEM_THREAD_ID(y);
+   const int bxy = MFEM_THREAD_SIZE(x)*MFEM_THREAD_SIZE(y);
+
+   for (int i = tid; i < ND; i += bxy)
+   {
+      Z(i, e) = a*X(i, e) + b*Y(i, e);
+   }
+   MFEM_SYNC_THREAD;
+}
+
+template <int NB>
+MFEM_HOST_DEVICE inline
+double DGMassDot(const int e,
+                 const int NE,
+                 const int ND,
+                 const double *x,
+                 const double *y)
+{
+   const auto X = ConstDeviceMatrix(x, ND, NE);
+   const auto Y = ConstDeviceMatrix(y, ND, NE);
+
+   const int tid = MFEM_THREAD_ID(x) + MFEM_THREAD_SIZE(x)*MFEM_THREAD_ID(y);
+   const int bxy = MFEM_THREAD_SIZE(x)*MFEM_THREAD_SIZE(y);
+
+   MFEM_SHARED double s_dot[NB*NB];
+   s_dot[tid] = 0.0;
+
+   for (int i = tid; i < ND; i += bxy) { s_dot[tid] += X(i,e)*Y(i,e); }
+   MFEM_SYNC_THREAD;
+
+   if (bxy > 512 && tid + 512 < bxy) { s_dot[tid] += s_dot[tid + 512]; }
+   MFEM_SYNC_THREAD;
+
+   if (bxy > 256 && tid < 256 && tid + 256 < bxy) { s_dot[tid] += s_dot[tid + 256]; }
+   MFEM_SYNC_THREAD;
+
+   if (bxy > 128 && tid < 128 && tid + 128 < bxy) { s_dot[tid] += s_dot[tid + 128]; }
+   MFEM_SYNC_THREAD;
+
+   if (bxy > 64 && tid < 64 && tid + 64 < bxy) { s_dot[tid] += s_dot[tid + 64]; }
+   MFEM_SYNC_THREAD;
+
+   if (bxy > 32 && tid < 32 && tid + 32 < bxy) { s_dot[tid] += s_dot[tid + 32]; }
+   MFEM_SYNC_THREAD;
+
+   if (bxy > 16 && tid < 16 && tid + 16 < bxy) { s_dot[tid] += s_dot[tid + 16]; }
+   MFEM_SYNC_THREAD;
+
+   if (bxy > 8 && tid < 8 && tid + 8 < bxy) { s_dot[tid] += s_dot[tid + 8]; }
+   MFEM_SYNC_THREAD;
+
+   if (bxy > 4 && tid < 4 && tid + 4 < bxy) { s_dot[tid] += s_dot[tid + 4]; }
+   MFEM_SYNC_THREAD;
+
+   if (bxy > 2 && tid < 2 && tid + 2 < bxy) { s_dot[tid] += s_dot[tid + 2]; }
+   MFEM_SYNC_THREAD;
+
+   if (bxy > 1 && tid < 1 && tid + 1 < bxy) { s_dot[tid] += s_dot[tid + 1]; }
+   MFEM_SYNC_THREAD;
+
+   return s_dot[0];
+}
+
+template<int T_D1D = 0, int MAX_D1D = 0>
+MFEM_HOST_DEVICE inline
+void DGMassBasis2D(const int e,
+                   const int NE,
+                   const double *b_,
+                   const double *x_,
+                   double *y_,
+                   const int d1d = 0)
+{
+   constexpr int MD1 = T_D1D ? T_D1D : MAX_D1D;
+   const int D1D = T_D1D ? T_D1D : d1d;
+
+   const auto b = Reshape(b_, D1D, D1D);
+   const auto x = Reshape(x_, D1D, D1D, NE);
+   auto y = Reshape(y_, D1D, D1D, NE);
+
+   MFEM_SHARED double sB[MD1*MD1];
+   MFEM_SHARED double sm0[MD1*MD1];
+   MFEM_SHARED double sm1[MD1*MD1];
+
+   kernels::internal::LoadB<MD1,MD1>(D1D,D1D,b,sB);
+
+   ConstDeviceMatrix B(sB, D1D,D1D);
+   DeviceMatrix DD(sm0, MD1, MD1);
+   DeviceMatrix DQ(sm1, MD1, MD1);
+   DeviceMatrix QQ(sm0, MD1, MD1);
+
+   kernels::internal::LoadX(e,D1D,x,DD);
+   kernels::internal::EvalX(D1D,D1D,B,DD,DQ);
+   kernels::internal::EvalY(D1D,D1D,B,DQ,QQ);
+   MFEM_SYNC_THREAD; // sync here to allow in-place evaluations
+   MFEM_FOREACH_THREAD(qy,y,D1D)
+   {
+      MFEM_FOREACH_THREAD(qx,x,D1D)
+      {
+         y(qx,qy,e) = QQ(qx,qy);
+      }
+   }
+   MFEM_SYNC_THREAD;
+}
+
+template<int T_D1D = 0, int MAX_D1D = 0>
+MFEM_HOST_DEVICE inline
+void DGMassBasis3D(const int e,
+                   const int NE,
+                   const double *b_,
+                   const double *x_,
+                   double *y_,
+                   const int d1d = 0)
+{
+   const int D1D = T_D1D ? T_D1D : d1d;
+
+   const auto b = Reshape(b_, D1D, D1D);
+   const auto x = Reshape(x_, D1D, D1D, D1D, NE);
+   auto y = Reshape(y_, D1D, D1D, D1D, NE);
+
+   constexpr int MD1 = T_D1D ? T_D1D : MAX_D1D;
+
+   MFEM_SHARED double sB[MD1*MD1];
+   MFEM_SHARED double sm0[MD1*MD1*MD1];
+   MFEM_SHARED double sm1[MD1*MD1*MD1];
+
+   kernels::internal::LoadB<MD1,MD1>(D1D,D1D,b,sB);
+
+   ConstDeviceMatrix B(sB, D1D,D1D);
+   DeviceCube DDD(sm0, MD1,MD1,MD1);
+   DeviceCube DDQ(sm1, MD1,MD1,MD1);
+   DeviceCube DQQ(sm0, MD1,MD1,MD1);
+   DeviceCube QQQ(sm1, MD1,MD1,MD1);
+
+   kernels::internal::LoadX(e,D1D,x,DDD);
+   kernels::internal::EvalX(D1D,D1D,B,DDD,DDQ);
+   kernels::internal::EvalY(D1D,D1D,B,DDQ,DQQ);
+   kernels::internal::EvalZ(D1D,D1D,B,DQQ,QQQ);
+   MFEM_SYNC_THREAD; // sync here to allow in-place evaluation
+   MFEM_FOREACH_THREAD(qz,z,D1D)
+   {
+      MFEM_FOREACH_THREAD(qy,y,D1D)
+      {
+         for (int qx = 0; qx < D1D; ++qx)
+         {
+            y(qx,qy,qz,e) = QQQ(qz,qy,qx);
+         }
+      }
+   }
+   MFEM_SYNC_THREAD;
+}
+
+template<int DIM, int T_D1D = 0, int MAX_D1D = 0>
+MFEM_HOST_DEVICE inline
+void DGMassBasis(const int e,
+                 const int NE,
+                 const double *b_,
+                 const double *x_,
+                 double *y_,
+                 const int d1d = 0)
+{
+   if (DIM == 2)
+   {
+      DGMassBasis2D<T_D1D, MAX_D1D>(e, NE, b_, x_, y_, d1d);
+   }
+   else if (DIM == 3)
+   {
+      DGMassBasis3D<T_D1D, MAX_D1D>(e, NE, b_, x_, y_, d1d);
+   }
+   else
+   {
+      MFEM_ABORT_KERNEL("Dimension not supported.");
+   }
+}
+
+} // namespace mfem
+
+#endif

fem/fespace.cpp

@@ -466,9 +466,11 @@ void FiniteElementSpace::BuildDofToArrays()
    dof_elem_array.SetSize (ndofs);
    dof_ldof_array.SetSize (ndofs);
    dof_elem_array = -1;
+   elem_dof->HostReadI();
+   elem_dof->HostReadJ();
    for (int i = 0; i < mesh -> GetNE(); i++)
    {
-      const int *dofs = elem_dof -> GetRow(i);
+      const int *dofs = AsConst(*elem_dof).GetRow(i);
       const int n = elem_dof -> RowSize(i);
       for (int j = 0; j < n; j++)
       {

fem/kernels.hpp

@@ -139,6 +139,22 @@ MFEM_HOST_DEVICE inline void LoadX(const int e, const int D1D,
    MFEM_SYNC_THREAD;
 }
 
+MFEM_HOST_DEVICE inline void LoadX(const int e, const int D1D,
+                                   const DeviceTensor<3, const double> &x,
+                                   DeviceMatrix &X)
+{
+   const int tidz = MFEM_THREAD_ID(z);
+
+   MFEM_FOREACH_THREAD(dy,y,D1D)
+   {
+      MFEM_FOREACH_THREAD(dx,x,D1D)
+      {
+         X(dx,dy) = x(dx,dy,e);
+      }
+   }
+   MFEM_SYNC_THREAD;
+}
+
 /// Load 2D input scalar into shared memory, with comp
 MFEM_HOST_DEVICE inline void LoadX(const int e, const int D1D, const int c,
                                    const DeviceTensor<4, const double> &x,

fem/pfespace.cpp

@@ -2630,7 +2630,7 @@ int ParFiniteElementSpace
                                 dof_offs, tdof_offs);
 
 #ifdef HYPRE_USING_GPU
-     hypre_ParCSRMatrixLocalTranspose(**P_);
+      hypre_ParCSRMatrixLocalTranspose(**P_);
 #endif
    }
 

fem/quadinterpolator_face.cpp

@@ -10,6 +10,7 @@
 // CONTRIBUTING.md for details.
 
 #include "quadinterpolator_face.hpp"
+#include "../general/annotation.hpp"
 #include "../general/forall.hpp"
 #include "../linalg/dtensor.hpp"
 #include "../linalg/kernels.hpp"
@@ -217,9 +218,6 @@ void FaceQuadratureInterpolator::Eval3D(
                "Derivatives on the faces are not yet supported.");
    MFEM_FORALL(f, NF,
    {
-      const int ND1D = T_ND1D ? T_ND1D : nd1d;
-      const int NQ1D = T_NQ1D ? T_NQ1D : nq1d;
-      const int VDIM = T_VDIM ? T_VDIM : vdim;
       constexpr int max_ND1D = T_ND1D ? T_ND1D : MAX_ND1D;
       constexpr int max_NQ1D = T_NQ1D ? T_NQ1D : MAX_NQ1D;
       constexpr int max_VDIM = T_VDIM ? T_VDIM : MAX_VDIM3D;
@@ -355,6 +353,184 @@ void FaceQuadratureInterpolator::Eval3D(
    });
 }
 
+template<const int T_VDIM, const int T_ND1D, const int T_NQ1D>
+void FaceQuadratureInterpolator::SmemEval3D(
+   const int NF,
+   const int vdim,
+   const DofToQuad &maps,
+   const Array<bool> &signs,
+   const Vector &e_vec,
+   Vector &q_val,
+   Vector &q_der,
+   Vector &q_det,
+   Vector &q_nor,
+   const int eval_flags)
+{
+   MFEM_PERF_SCOPE("FaceQuadInterpolator::SmemEval3D");
+   const int nd1d = maps.ndof;
+   const int nq1d = maps.nqpt;
+   const int ND1D = T_ND1D ? T_ND1D : nd1d;
+   const int NQ1D = T_NQ1D ? T_NQ1D : nq1d;
+   const int VDIM = T_VDIM ? T_VDIM : vdim;
+   MFEM_VERIFY(ND1D <= MAX_ND1D, "");
+   MFEM_VERIFY(NQ1D <= MAX_NQ1D, "");
+   MFEM_VERIFY(VDIM == 3 || !(eval_flags & DETERMINANTS), "");
+   auto B = Reshape(maps.B.Read(), NQ1D, ND1D);
+   auto G = Reshape(maps.G.Read(), NQ1D, ND1D);
+   auto F = Reshape(e_vec.Read(), ND1D, ND1D, VDIM, NF);
+   auto sign = signs.Read();
+   auto val = Reshape(q_val.Write(), NQ1D, NQ1D, VDIM, NF);
+   // auto der = Reshape(q_der.Write(), NQ1D, VDIM, 3, NF);
+   auto det = Reshape(q_det.Write(), NQ1D, NQ1D, NF);
+   auto nor = Reshape(q_nor.Write(), NQ1D, NQ1D, 3, NF);
+   MFEM_VERIFY(eval_flags | DERIVATIVES,
+               "Derivatives on the faces are not yet supported.");
+
+   MFEM_FORALL_3D(f, NF, NQ1D, NQ1D, VDIM,
+   {
+      constexpr int max_ND1D = T_ND1D ? T_ND1D : MAX_ND1D;
+      constexpr int max_NQ1D = T_NQ1D ? T_NQ1D : MAX_NQ1D;
+      constexpr int max_VDIM = T_VDIM ? T_VDIM : MAX_VDIM3D;
+
+      MFEM_SHARED double sm1[max_NQ1D*max_NQ1D*max_VDIM];
+      MFEM_SHARED double sm2[max_NQ1D*max_ND1D*max_VDIM];
+
+      auto s_F = (double(*)[max_ND1D][max_VDIM])sm1;
+      MFEM_FOREACH_THREAD(d1,x,ND1D)
+      {
+         MFEM_FOREACH_THREAD(d2,y,ND1D)
+         {
+            MFEM_FOREACH_THREAD(c,z,VDIM)
+            {
+               s_F[d1][d2][c] = F(d1,d2,c,f);
+            }
+         }
+      }
+      MFEM_SYNC_THREAD;
+
+      if (eval_flags & VALUES)
+      {
+         auto Bu = (double (*)[max_ND1D][max_VDIM])sm2;
+         MFEM_FOREACH_THREAD(d2,x,ND1D)
+         {
+            MFEM_FOREACH_THREAD(q1,y,NQ1D)
+            {
+               MFEM_FOREACH_THREAD(c,z,VDIM)
+               {
+                  double thrdBu = 0.0;
+                  for (int d1 = 0; d1 < ND1D; ++d1)
+                  {
+                     thrdBu += B(q1,d1)*s_F[d1][d2][c];
+                  }
+                  Bu[q1][d2][c] = thrdBu;
+               }
+            }
+         }
+         MFEM_SYNC_THREAD;
+
+         MFEM_FOREACH_THREAD(q2,x,NQ1D)
+         {
+            MFEM_FOREACH_THREAD(q1,y,NQ1D)
+            {
+               MFEM_FOREACH_THREAD(c,z,VDIM)
+               {
+                  double v = 0.0;
+                  for (int d2 = 0; d2 < ND1D; ++d2)
+                  {
+                     v += B(q2,d2)*Bu[q1][d2][c];
+                  }
+                  val(q1,q2,c,f) = v;
+               }
+            }
+         }
+      }
+
+      if ((eval_flags & DERIVATIVES)
+          || (eval_flags & DETERMINANTS)
+          || (eval_flags & NORMALS))
+      {
+         // We only compute the tangential derivatives
+         auto Gu = (double (*)[max_ND1D][max_VDIM])sm2;
+         MFEM_SHARED double Bu[max_NQ1D][max_ND1D][max_VDIM];
+         MFEM_FOREACH_THREAD(d2,x,ND1D)
+         {
+            MFEM_FOREACH_THREAD(q1,y,NQ1D)
+            {
+               MFEM_FOREACH_THREAD(c,z,VDIM)
+               {
+                  double thrdGu = 0;
+                  double thrdBu = 0;
+                  for (int d1 = 0; d1 < ND1D; ++d1)
+                  {
+                     const double u = s_F[d1][d2][c];
+                     thrdBu += B(q1,d1)*u;
+                     thrdGu += G(q1,d1)*u;
+                  }
+                  Gu[q1][d2][c] = thrdGu;
+                  Bu[q1][d2][c] = thrdBu;
+               }
+            }
+         }
+         MFEM_SYNC_THREAD;
+
+         auto BGu = (double (*)[max_NQ1D][max_VDIM])sm1;
+         MFEM_SHARED double GBu[max_NQ1D][max_NQ1D][max_VDIM];
+         MFEM_FOREACH_THREAD(q2,x,NQ1D)
+         {
+            MFEM_FOREACH_THREAD(q1,y,NQ1D)
+            {
+               MFEM_FOREACH_THREAD(c,z,VDIM)
+               {
+                  double thrdBGu = 0.0;
+                  double thrdGBu = 0.0;
+                  for (int d2 = 0; d2 < ND1D; ++d2)
+                  {
+                     thrdBGu += B(q2,d2)*Gu[q1][d2][c];
+                     thrdGBu += G(q2,d2)*Bu[q1][d2][c];
+                  }
+                  BGu[q2][q1][c] = thrdBGu;
+                  GBu[q2][q1][c] = thrdGBu;
+               }
+            }
+         }
+         MFEM_SYNC_THREAD;
+
+         if (VDIM == 3 && ((eval_flags & NORMALS) ||
+                           (eval_flags & DETERMINANTS)))
+         {
+            double n[3];
+            MFEM_FOREACH_THREAD(q2,x,NQ1D)
+            {
+               MFEM_FOREACH_THREAD(q1,y,NQ1D)
+               {
+                  if (MFEM_THREAD_ID(z) == 0)
+                  {
+                     const double s = sign[f] ? -1.0 : 1.0;
+                     n[0] = s*( BGu[q2][q1][1]*GBu[q2][q1][2]-GBu[q2][q1][1]*
+                                BGu[q2][q1][2] );
+                     n[1] = s*(-BGu[q2][q1][0]*GBu[q2][q1][2]+GBu[q2][q1][0]*
+                               BGu[q2][q1][2] );
+                     n[2] = s*( BGu[q2][q1][0]*GBu[q2][q1][1]-GBu[q2][q1][0]*
+                                BGu[q2][q1][1] );
+
+                     const double norm = sqrt(n[0]*n[0]+n[1]*n[1]+n[2]*n[2]);
+
+                     if (eval_flags & DETERMINANTS) { det(q1,q2,f) = norm; }
+
+                     if (eval_flags & NORMALS)
+                     {
+                        nor(q1,q2,0,f) = n[0]/norm;
+                        nor(q1,q2,1,f) = n[1]/norm;
+                        nor(q1,q2,2,f) = n[2]/norm;
+                     }
+                  }
+               }
+            }
+         }
+      }
+   });
+}
+
 void FaceQuadratureInterpolator::Mult(
    const Vector &e_vec, unsigned eval_flags,
    Vector &q_val, Vector &q_der, Vector &q_det, Vector &q_nor) const
@@ -414,21 +590,22 @@ void FaceQuadratureInterpolator::Mult(
          switch (10*nd1d + nq1d)
          {
             // Q0
-            case 11: eval_func = &Eval3D<1,1,1>; break;
-            case 12: eval_func = &Eval3D<1,1,2>; break;
+            case 11: eval_func = &SmemEval3D<1,1,1>; break;
+            case 12: eval_func = &SmemEval3D<1,1,2>; break;
             // Q1
-            case 22: eval_func = &Eval3D<1,2,2>; break;
-            case 23: eval_func = &Eval3D<1,2,3>; break;
+            case 22: eval_func = &SmemEval3D<1,2,2>; break;
+            case 23: eval_func = &SmemEval3D<1,2,3>; break;
+            case 24: eval_func = &SmemEval3D<1,2,4>; break;
             // Q2
-            case 33: eval_func = &Eval3D<1,3,3>; break;
-            case 34: eval_func = &Eval3D<1,3,4>; break;
+            case 33: eval_func = &SmemEval3D<1,3,3>; break;
+            case 34: eval_func = &SmemEval3D<1,3,4>; break;
             // Q3
-            case 44: eval_func = &Eval3D<1,4,4>; break;
-            case 45: eval_func = &Eval3D<1,4,5>; break;
-            case 46: eval_func = &Eval3D<1,4,6>; break;
+            case 44: eval_func = &SmemEval3D<1,4,4>; break;
+            case 45: eval_func = &SmemEval3D<1,4,5>; break;
+            case 46: eval_func = &SmemEval3D<1,4,6>; break;
             // Q4
-            case 55: eval_func = &Eval3D<1,5,5>; break;
-            case 56: eval_func = &Eval3D<1,5,6>; break;
+            case 55: eval_func = &SmemEval3D<1,5,5>; break;
+            case 56: eval_func = &SmemEval3D<1,5,6>; break;
          }
          if (nq1d >= 10 || !eval_func)
          {
@@ -468,18 +645,19 @@ void FaceQuadratureInterpolator::Mult(
          switch (10*nd1d + nq1d)
          {
             // Q1
-            case 22: eval_func = &Eval3D<3,2,2>; break;
-            case 23: eval_func = &Eval3D<3,2,3>; break;
+            case 22: eval_func = &SmemEval3D<3,2,2>; break;
+            case 23: eval_func = &SmemEval3D<3,2,3>; break;
+            case 24: eval_func = &SmemEval3D<3,2,4>; break;
             // Q2
-            case 33: eval_func = &Eval3D<3,3,3>; break;
-            case 34: eval_func = &Eval3D<3,3,4>; break;
+            case 33: eval_func = &SmemEval3D<3,3,3>; break;
+            case 34: eval_func = &SmemEval3D<3,3,4>; break;
             // Q3
-            case 44: eval_func = &Eval3D<3,4,4>; break;
-            case 45: eval_func = &Eval3D<3,4,5>; break;
-            case 46: eval_func = &Eval3D<3,4,6>; break;
+            case 44: eval_func = &SmemEval3D<3,4,4>; break;
+            case 45: eval_func = &SmemEval3D<3,4,5>; break;
+            case 46: eval_func = &SmemEval3D<3,4,6>; break;
             // Q4
-            case 55: eval_func = &Eval3D<3,5,5>; break;
-            case 56: eval_func = &Eval3D<3,5,6>; break;
+            case 55: eval_func = &SmemEval3D<3,5,5>; break;
+            case 56: eval_func = &SmemEval3D<3,5,6>; break;
          }
          if (nq1d >= 10 || !eval_func)
          {

fem/quadinterpolator_face.hpp

@@ -112,6 +112,18 @@ public:
                       Vector &q_det,
                       Vector &q_nor,
                       const int eval_flags);
+
+   template<const int T_VDIM = 0, const int T_ND = 0, const int T_NQ = 0>
+   static void SmemEval3D(const int NF,
+                          const int vdim,
+                          const DofToQuad &maps,
+                          const Array<bool> &signs,
+                          const Vector &e_vec,
+                          Vector &q_val,
+                          Vector &q_der,
+                          Vector &q_det,
+                          Vector &q_nor,
+                          const int eval_flags);
 };
 
 } // mfem namespace

fem/restriction.hpp

@@ -47,6 +47,7 @@ protected:
    friend class BatchedLOR_ADS;
    friend class BatchedLOR_AMS;
 
+public:
    /// @name Low-level access to the underlying element-dof mappings
    ///@{
    const Array<int> &GatherMap() const { return gather_map; }
@@ -54,7 +55,6 @@ protected:
    const Array<int> &Offsets() const { return offsets; }
    ///@}
 
-public:
    ElementRestriction(const FiniteElementSpace&, ElementDofOrdering);
    void Mult(const Vector &x, Vector &y) const;
    void MultTranspose(const Vector &x, Vector &y) const;

linalg/hypre.cpp

@@ -1834,7 +1834,8 @@ void HypreParMatrix::MultTranspose(double a, const Vector &x,
    }
 
 #ifdef HYPRE_USING_GPU
-   MFEM_VERIFY(A->diagT != NULL, "Transpose action requires hypre_ParCSRMatrixLocalTranspose");
+   MFEM_VERIFY(A->diagT != NULL,
+               "Transpose action requires hypre_ParCSRMatrixLocalTranspose");
 #endif
 
    hypre_ParCSRMatrixMatvecT(a, A, *Y, b, *X);

linalg/hypre_parcsr.cpp

@@ -1560,8 +1560,8 @@ void hypre_ParCSRMatrixAbsMatvecT(hypre_ParCSRMatrix *A,
 
    if (num_cols_offd)
    {
-#if MFEM_HYPRE_VERSION >= 21100
-      if (A->offdT)
+#if MFEM_HYPRE_VERSION >= 21100 && 0
+      if (A->offdT) //need a gpu version
       {
          // offdT is optional. Used only if it's present.
          hypre_CSRMatrixAbsMatvec(A->offdT, alpha, x, 0., y_tmp);
@@ -1575,8 +1575,8 @@ void hypre_ParCSRMatrixAbsMatvecT(hypre_ParCSRMatrix *A,
 
    comm_handle = hypre_ParCSRCommHandleCreate(2, comm_pkg, y_tmp, y_buf);
 
-#if MFEM_HYPRE_VERSION >= 21100
-   if (A->diagT)
+#if MFEM_HYPRE_VERSION >= 21100 && 0
+   if (A->diagT) //need a gpu version
    {
       // diagT is optional. Used only if it's present.
       hypre_CSRMatrixAbsMatvec(A->diagT, alpha, x, beta, y);
@@ -1695,10 +1695,11 @@ void hypre_ParCSRMatrixBooleanMatvecT(hypre_ParCSRMatrix *A,
    y_buf = mfem_hypre_CTAlloc_host(
               HYPRE_Bool, hypre_ParCSRCommPkgSendMapStart(comm_pkg, num_sends));
 
+
    if (num_cols_offd)
    {
-#if MFEM_HYPRE_VERSION >= 21100
-      if (A->offdT)
+#if MFEM_HYPRE_VERSION >= 21100 && 0
+      if (A->offdT) //need a gpu version
       {
          // offdT is optional. Used only if it's present.
          hypre_CSRMatrixBooleanMatvec(A->offdT, alpha, x, 0, y_tmp);
@@ -1712,10 +1713,10 @@ void hypre_ParCSRMatrixBooleanMatvecT(hypre_ParCSRMatrix *A,
 
    comm_handle = hypre_ParCSRCommHandleCreate_bool(2, comm_pkg, y_tmp, y_buf);
 
-#if MFEM_HYPRE_VERSION >= 21100
-   if (A->diagT)
+#if MFEM_HYPRE_VERSION >= 21100 && 0
+   if (A->diagT) //need a gpu version
    {
-      // diagT is optional. Used only if it's present.
+      //diagT is optional. Used only if it's present.
       hypre_CSRMatrixBooleanMatvec(A->diagT, alpha, x, beta, y);
    }
    else