DG mass inverse: support both map type VALUE and INTEGRAL

.gitignore

@@ -320,9 +320,13 @@ tests/unit/debug_device_tests
 
 # Benchmark binaries
 tests/benchmarks/bench_ceed
+tests/benchmarks/bench_dgmassinv
 tests/benchmarks/bench_tmop
 tests/benchmarks/bench_vector
 tests/benchmarks/bench_virtuals
+tests/benchmarks/plot/*.pdf
+tests/benchmarks/plot/*.tex
+tests/benchmarks/plot/*.csv
 
 # Test script output
 tests/scripts/*.err

config/defaults.mk

@@ -453,7 +453,7 @@ GSLIB_LIB = -L$(GSLIB_DIR)/lib -lgs
 
 # CUDA library configuration
 CUDA_OPT =
-CUDA_LIB = -lcusparse
+CUDA_LIB = -lcusparse -lcusolver -lcublas
 
 # HIP library configuration
 HIP_OPT =

fem/CMakeLists.txt

@@ -125,6 +125,7 @@ set(HDRS
   bilinearform.hpp
   bilinearform_ext.hpp
   bilininteg.hpp
+  bilininteg_mass_pa.hpp
   coefficient.hpp
   complex_fem.hpp
   convergence.hpp

fem/bilininteg.hpp

@@ -2170,6 +2170,8 @@ public:
 class MassIntegrator: public BilinearFormIntegrator
 {
 protected:
+   friend class DGMassInverse;
+   friend class DGMassInverse_Direct;
 #ifndef MFEM_THREAD_SAFE
    Vector shape, te_shape;
 #endif

fem/bilininteg_mass_pa.cpp

@@ -13,6 +13,7 @@
 #include "bilininteg.hpp"
 #include "gridfunc.hpp"
 #include "ceed/mass.hpp"
+#include "bilininteg_mass_pa.hpp"
 
 using namespace std;
 
@@ -581,85 +582,18 @@ static void PAMassApply2D(const int NE,
                           const int d1d = 0,
                           const int q1d = 0)
 {
-   const int D1D = T_D1D ? T_D1D : d1d;
-   const int Q1D = T_Q1D ? T_Q1D : q1d;
-   MFEM_VERIFY(D1D <= MAX_D1D, "");
-   MFEM_VERIFY(Q1D <= MAX_Q1D, "");
-   auto B = Reshape(b_.Read(), Q1D, D1D);
-   auto Bt = Reshape(bt_.Read(), D1D, Q1D);
-   auto D = Reshape(d_.Read(), Q1D, Q1D, NE);
-   auto X = Reshape(x_.Read(), D1D, D1D, NE);
-   auto Y = Reshape(y_.ReadWrite(), D1D, D1D, NE);
+   MFEM_VERIFY(T_D1D ? T_D1D : d1d <= MAX_D1D, "");
+   MFEM_VERIFY(T_Q1D ? T_Q1D : q1d <= MAX_Q1D, "");
+
+   const auto B = b_.Read();
+   const auto Bt = bt_.Read();
+   const auto D = d_.Read();
+   const auto X = x_.Read();
+   auto Y = y_.ReadWrite();
+
    MFEM_FORALL(e, NE,
    {
-      const int D1D = T_D1D ? T_D1D : d1d; // nvcc workaround
-      const int Q1D = T_Q1D ? T_Q1D : q1d;
-      // the following variables are evaluated at compile time
-      constexpr int max_D1D = T_D1D ? T_D1D : MAX_D1D;
-      constexpr int max_Q1D = T_Q1D ? T_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];
-            }
-         }
-      }
+      PAMassApply2D_Element(e, NE, B, Bt, D, X, Y, d1d, q1d);
    });
 }
 
@@ -681,108 +615,13 @@ static void SmemPAMassApply2D(const int NE,
    constexpr int MD1 = T_D1D ? T_D1D : MAX_D1D;
    MFEM_VERIFY(D1D <= MD1, "");
    MFEM_VERIFY(Q1D <= MQ1, "");
-   auto b = Reshape(b_.Read(), Q1D, D1D);
-   auto D = Reshape(d_.Read(), Q1D, Q1D, NE);
-   auto x = Reshape(x_.Read(), D1D, D1D, NE);
-   auto Y = Reshape(y_.ReadWrite(), D1D, D1D, NE);
+   const auto b = b_.Read();
+   const auto D = d_.Read();
+   const auto x = x_.Read();
+   auto Y = y_.ReadWrite();
    MFEM_FORALL_2D(e, NE, Q1D, Q1D, NBZ,
    {
-      const int tidz = MFEM_THREAD_ID(z);
-      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;
-      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]);
-            }
-            Y(dx, dy, e) += dd;
-         }
-      }
+      SmemPAMassApply2D_Element<T_D1D,T_Q1D,T_NBZ>(e, NE, b, D, x, Y, d1d, q1d);
    });
 }
 
@@ -796,134 +635,18 @@ static void PAMassApply3D(const int NE,
                           const int d1d = 0,
                           const int q1d = 0)
 {
-   const int D1D = T_D1D ? T_D1D : d1d;
-   const int Q1D = T_Q1D ? T_Q1D : q1d;
-   MFEM_VERIFY(D1D <= MAX_D1D, "");
-   MFEM_VERIFY(Q1D <= MAX_Q1D, "");
-   auto B = Reshape(b_.Read(), Q1D, D1D);
-   auto Bt = Reshape(bt_.Read(), D1D, Q1D);
-   auto D = Reshape(d_.Read(), Q1D, Q1D, Q1D, NE);
-   auto X = Reshape(x_.Read(), D1D, D1D, D1D, NE);
-   auto Y = Reshape(y_.ReadWrite(), D1D, D1D, D1D, NE);
+   MFEM_VERIFY(T_D1D ? T_D1D : d1d <= MAX_D1D, "");
+   MFEM_VERIFY(T_Q1D ? T_Q1D : q1d <= MAX_Q1D, "");
+
+   const auto B = b_.Read();
+   const auto Bt = bt_.Read();
+   const auto D = d_.Read();
+   const auto X = x_.Read();
+   auto Y = y_.ReadWrite();
+
    MFEM_FORALL(e, NE,
    {
-      const int D1D = T_D1D ? T_D1D : d1d;
-      const int Q1D = T_Q1D ? T_Q1D : q1d;
-      constexpr int max_D1D = T_D1D ? T_D1D : MAX_D1D;
-      constexpr int max_Q1D = T_Q1D ? T_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];
-               }
-            }
-         }
-      }
+      PAMassApply3D_Element(e, NE, B, Bt, D, X, Y, d1d, q1d);
    });
 }
 
@@ -944,213 +667,13 @@ static void SmemPAMassApply3D(const int NE,
    constexpr int M1D = T_D1D ? T_D1D : MAX_D1D;
    MFEM_VERIFY(D1D <= M1D, "");
    MFEM_VERIFY(Q1D <= M1Q, "");
-   auto b = Reshape(b_.Read(), Q1D, D1D);
-   auto d = Reshape(d_.Read(), Q1D, Q1D, Q1D, NE);
-   auto x = Reshape(x_.Read(), D1D, D1D, D1D, NE);
-   auto y = Reshape(y_.ReadWrite(), D1D, D1D, D1D, NE);
+   auto b = b_.Read();
+   auto d = d_.Read();
+   auto x = x_.Read();
+   auto y = y_.ReadWrite();
    MFEM_FORALL_3D(e, NE, Q1D, Q1D, 1,
    {
-      const int D1D = T_D1D ? T_D1D : d1d;
-      const 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;
-      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(d,y,D1D)
-      {
-         MFEM_FOREACH_THREAD(q,x,Q1D)
-         {
-            Bt[d][q] = b(q,d);
-         }
-      }
-      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)
-            {
-               y(dx,dy,dz,e) += u[dz];
-            }
-         }
-      }
+      SmemPAMassApply3D_Element<T_D1D,T_Q1D>(e, NE, b, d, x, y, d1d, q1d);
    });
 }
 

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,608 @@
+// 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
+{
+   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 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);
+         // 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;
+
+      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
+
+      // 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 nom0 = DGMassDot<NB>(e, NE, ND, d, r);
+      double nom = nom0;
+      // 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 0x22: return DGMassCGIteration<2,2,2>(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 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 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,
+                                           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
+{
+   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,180 @@
+// 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_);
+   ~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_,
+                        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,290 @@
+// 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_mass_pa.hpp"
+#include "../linalg/kernels.hpp"
+#include "kernels.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/fe/fe_l2.hpp

@@ -55,6 +55,7 @@ public:
                             ElementTransformation &Trans,
                             DenseMatrix &curl) const
    { ProjectCurl_2D(fe, Trans, curl); }
+   using FiniteElement::Project;
    virtual void ProjectDiv(const FiniteElement &fe,
                            ElementTransformation &Trans,
                            DenseMatrix &div) const;
@@ -79,6 +80,7 @@ public:
    virtual void CalcDShape(const IntegrationPoint &ip,
                            DenseMatrix &dshape) const;
    virtual void ProjectDelta(int vertex, Vector &dofs) const;
+   using FiniteElement::Project;
    virtual void ProjectDiv(const FiniteElement &fe,
                            ElementTransformation &Trans,
                            DenseMatrix &div) const;

fem/fem.hpp

@@ -45,6 +45,7 @@
 #include "multigrid.hpp"
 #include "ceed/algebraic.hpp"
 #include "lor.hpp"
+#include "dgmassinv.hpp"
 
 #ifdef MFEM_USE_MPI
 #include "pfespace.hpp"

fem/kernels.hpp

@@ -139,6 +139,20 @@ 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 &DD)
+{
+   MFEM_FOREACH_THREAD(dy,y,D1D)
+   {
+      MFEM_FOREACH_THREAD(dx,x,D1D)
+      {
+         DD(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,

linalg/dtensor.hpp

@@ -127,7 +127,7 @@ public:
 
 /** @brief Wrap a pointer as a DeviceTensor with automatically deduced template
     parameters */
-template <typename T, typename... Dims>
+template <typename T, typename... Dims> MFEM_HOST_DEVICE
 inline DeviceTensor<sizeof...(Dims),T> Reshape(T *ptr, Dims... dims)
 {
    return DeviceTensor<sizeof...(Dims),T>(ptr, dims...);

linalg/solvers.cpp

@@ -861,7 +861,7 @@ void CGSolver::Mult(const Vector &b, Vector &x) const
       }
       nom = betanom;
    }
-   if (print_options.first_and_last)
+   if (print_options.first_and_last && !print_options.iterations)
    {
       mfem::out << "   Iteration : " << setw(3) << final_iter << "  (B r, r) = "
                 << betanom << '\n';

tests/benchmarks/bench_dgmassinv.cpp

@@ -0,0 +1,277 @@
+// 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 "bench.hpp"
+#include "kershaw.hpp"
+#include <memory>
+
+#ifdef MFEM_USE_BENCHMARK
+
+static constexpr double tol = 1e-10;
+
+enum class MassSolverType
+{
+   FULL_CG,
+   LOCAL_CG_LOBATTO,
+   LOCAL_CG_LEGENDRE,
+   DIRECT,
+   DIRECT_CUSOLVER,
+   DIRECT_CUBLAS
+};
+
+Mesh CreateKershawMesh(int N, double eps)
+{
+   Mesh mesh = Mesh::MakeCartesian3D(N, N, N, Element::HEXAHEDRON);
+   KershawTransformation kt(mesh.Dimension(), eps, eps);
+   mesh.Transform(kt);
+   return mesh;
+}
+
+struct DGMassInverse_FullCG : Solver
+{
+   BilinearForm m;
+   OperatorJacobiSmoother jacobi;
+   CGSolver cg;
+
+   DGMassInverse_FullCG(FiniteElementSpace &fes) : m(&fes)
+   {
+      m.AddDomainIntegrator(new MassIntegrator);
+      m.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+      m.Assemble();
+
+      jacobi.SetOperator(m);
+
+      cg.SetAbsTol(tol);
+      cg.SetRelTol(0.0);
+      cg.SetMaxIter(100);
+      cg.SetOperator(m);
+      cg.SetPreconditioner(jacobi);
+   }
+
+   void Mult(const Vector &b, Vector &x) const
+   {
+      cg.Mult(b, x);
+   }
+
+   void Setup()
+   {
+      m.Update();
+      m.Assemble();
+
+      jacobi.SetOperator(m);
+   }
+
+   void SetOperator(const Operator &op) { }
+};
+
+struct DGMassBenchmark
+{
+   MassSolverType solver_type;
+   const int p;
+   const int N;
+   const int dim = 3;
+   Mesh mesh;
+   L2_FECollection fec;
+   FiniteElementSpace fes;
+   const int n;
+
+   std::unique_ptr<Solver> massinv;
+
+   Vector B, X;
+
+   const int dofs;
+   double mdofs;
+
+   DGMassBenchmark(MassSolverType type_, int p_, int N_, double eps_) :
+      solver_type(type_),
+      p(p_),
+      N(N_),
+      mesh(CreateKershawMesh(N,eps_)),
+      fec(p, dim, BasisType::Positive),
+      fes(&mesh, &fec),
+      n(fes.GetTrueVSize()),
+      B(n),
+      X(n),
+      dofs(n),
+      mdofs(0.0)
+   {
+      B.Randomize(1);
+      tic_toc.Clear();
+      NewSolver();
+   }
+
+   void NewFullCG()
+   {
+      if (massinv)
+      {
+         static_cast<DGMassInverse_FullCG*>(massinv.get())->Setup();
+      }
+      else
+      {
+         massinv.reset(new DGMassInverse_FullCG(fes));
+      }
+   }
+
+   void NewLocalCG(int btype)
+   {
+      if (massinv)
+      {
+         static_cast<DGMassInverse*>(massinv.get())->Setup();
+      }
+      else
+      {
+         DGMassInverse *massinv_ = new DGMassInverse(fes, btype);
+         massinv_->SetAbsTol(tol);
+         massinv_->SetRelTol(0.0);
+         massinv.reset(massinv_);
+      }
+   }
+
+   void NewDirect(BatchSolverMode mode)
+   {
+      if (massinv)
+      {
+         static_cast<DGMassInverse_Direct*>(massinv.get())->Setup();
+      }
+      else
+      {
+         DGMassInverse_Direct *massinv_ = new DGMassInverse_Direct(fes, mode);
+         massinv.reset(massinv_);
+      }
+   }
+
+   void NewSolver()
+   {
+      switch (solver_type)
+      {
+         case MassSolverType::FULL_CG:
+            NewFullCG();
+            break;
+         case MassSolverType::LOCAL_CG_LOBATTO:
+            NewLocalCG(BasisType::GaussLobatto);
+            break;
+         case MassSolverType::LOCAL_CG_LEGENDRE:
+            NewLocalCG(BasisType::GaussLegendre);
+            break;
+         case MassSolverType::DIRECT:
+            NewDirect(BatchSolverMode::NATIVE);
+            break;
+         case MassSolverType::DIRECT_CUBLAS:
+            NewDirect(BatchSolverMode::CUBLAS);
+            break;
+         case MassSolverType::DIRECT_CUSOLVER:
+            NewDirect(BatchSolverMode::CUSOLVER);
+            break;
+      }
+   }
+
+   void Setup()
+   {
+      MFEM_DEVICE_SYNC;
+      tic_toc.Start();
+      NewSolver();
+      MFEM_DEVICE_SYNC;
+      tic_toc.Stop();
+      mdofs += 1e-6 * dofs;
+   }
+
+   void Solve()
+   {
+      if (!massinv) { NewSolver(); }
+      X = 0.0;
+      MFEM_DEVICE_SYNC;
+      tic_toc.Start();
+      massinv->Mult(B, X);
+      MFEM_DEVICE_SYNC;
+      tic_toc.Stop();
+      mdofs += 1e-6 * dofs;
+   }
+
+   void SetupAndSolve()
+   {
+      X = 0.0;
+      MFEM_DEVICE_SYNC;
+      tic_toc.Start();
+      NewSolver();
+      massinv->Mult(B, X);
+      MFEM_DEVICE_SYNC;
+      tic_toc.Stop();
+      mdofs += 1e-6 * dofs;
+   }
+
+   double Mdofs() const { return mdofs / tic_toc.RealTime(); }
+};
+
+// The different orders the tests can run
+#define P_ORDERS bm::CreateDenseRange(1,5,1)
+// The different sides of the mesh
+#define N_SIDES bm::CreateDenseRange(2,30,2)
+
+#define MAX_NDOFS 2*1024*1024
+
+/// Kernels definitions and registrations
+#define Benchmark(solver_type, op_name, suffix, eps)\
+static void solver_type##_##op_name##_##suffix(bm::State &state){\
+   const int side = state.range(0);\
+   const int p = state.range(1);\
+   DGMassBenchmark mb(MassSolverType:: solver_type, p, side, eps);\
+   if (mb.dofs > MAX_NDOFS) { state.SkipWithError("MAX_NDOFS"); }\
+   while (state.KeepRunning()) { mb.op_name(); }\
+   state.counters["MDof/s"] = bm::Counter(mb.Mdofs());\
+   state.counters["dofs"] = bm::Counter(mb.dofs);\
+   state.counters["p"] = bm::Counter(p);\
+}\
+BENCHMARK(solver_type##_##op_name##_##suffix)\
+            -> ArgsProduct({N_SIDES,P_ORDERS})\
+            -> Unit(bm::kMillisecond);
+
+#define MassBenchmarks(solver_type, suffix, eps) \
+   Benchmark(solver_type, Setup, suffix, eps) \
+   Benchmark(solver_type, Solve, suffix, eps) \
+   Benchmark(solver_type, SetupAndSolve, suffix, eps)
+
+#define AllMassBenchmarks(suffix, eps) \
+   MassBenchmarks(FULL_CG, suffix, eps) \
+   MassBenchmarks(LOCAL_CG_LOBATTO, suffix, eps) \
+   MassBenchmarks(LOCAL_CG_LEGENDRE, suffix, eps) \
+   MassBenchmarks(DIRECT, suffix, eps) \
+   MassBenchmarks(DIRECT_CUBLAS, suffix, eps) \
+   MassBenchmarks(DIRECT_CUSOLVER, suffix, eps)
+
+AllMassBenchmarks(1_0, 1.0)
+AllMassBenchmarks(0_5, 0.5)
+AllMassBenchmarks(0_3, 0.3)
+
+int main(int argc, char *argv[])
+{
+   bm::ConsoleReporter CR;
+   bm::Initialize(&argc, argv);
+
+   // Device setup, cpu by default
+   std::string device_config = "cpu";
+   if (bmi::global_context != nullptr)
+   {
+      const auto device = bmi::global_context->find("device");
+      if (device != bmi::global_context->end())
+      {
+         mfem::out << device->first << " : " << device->second << std::endl;
+         device_config = device->second;
+      }
+   }
+   Device device(device_config.c_str());
+   device.Print();
+
+   if (bm::ReportUnrecognizedArguments(argc, argv)) { return 1; }
+   bm::RunSpecifiedBenchmarks(&CR);
+   return 0;
+}
+
+#endif // MFEM_USE_BENCHMARK

tests/benchmarks/kershaw.hpp

@@ -0,0 +1,113 @@
+#ifndef __KERSHAW_HPP__
+#define __KERSHAW_HPP__
+
+#include "mfem.hpp"
+
+// 1D transformation at the right boundary.
+double right(const double eps, const double x)
+{
+   return (x <= 0.5) ? (2-eps) * x : 1 + eps*(x-1);
+}
+
+// 1D transformation at the left boundary
+double left(const double eps, const double x)
+{
+   return 1-right(eps,1-x);
+}
+
+// Transition from a value of "a" for x=0, to a value of "b" for x=1. Smoothness
+// is controlled by the parameter "s", taking values 0, 1, or 2.
+double step(const double a, const double b, double x, int s)
+{
+   if (x <= 0) { return a; }
+   if (x >= 1) { return b; }
+   switch (s)
+   {
+      case 0:
+      default:
+         return a + (b-a) * (x);
+      case 1: return a + (b-a) * (x*x*(3-2*x));
+      case 2: return a + (b-a) * (x*x*x*(x*(6*x-15)+10));
+   }
+}
+
+// 3D version of a generalized Kershaw mesh transformation, see D. Kershaw,
+// "Differencing of the diffusion equation in Lagrangian hydrodynamic codes",
+// JCP, 39:375–395, 1981.
+//
+// The input mesh should be Cartesian nx x ny x nz with nx divisible by 6 and
+// ny, nz divisible by 2.
+//
+// The eps parameters are in (0, 1]. Uniform mesh is recovered for epsy=epsz=1.
+void kershaw(const double epsy, const double epsz, const int smoothness,
+             const double x, const double y, const double z,
+             double &X, double &Y, double &Z)
+{
+   X = x;
+
+   int layer = x*6.0;
+   double lambda = (x-layer/6.0)*6;
+
+   // The x-range is split in 6 layers going from left-to-left, left-to-right,
+   // right-to-left (2 layers), left-to-right and right-to-right yz-faces.
+   switch (layer)
+   {
+      case 0:
+         Y = left(epsy, y);
+         Z = left(epsz, z);
+         break;
+      case 1:
+      case 4:
+         Y = step(left(epsy, y), right(epsy, y), lambda, smoothness);
+         Z = step(left(epsz, z), right(epsz, z), lambda, smoothness);
+         break;
+      case 2:
+         Y = step(right(epsy, y), left(epsy, y), lambda/2, smoothness);
+         Z = step(right(epsz, z), left(epsz, z), lambda/2, smoothness);
+         break;
+      case 3:
+         Y = step(right(epsy, y), left(epsy, y), (1+lambda)/2, smoothness);
+         Z = step(right(epsz, z), left(epsz, z), (1+lambda)/2, smoothness);
+         break;
+      default:
+         Y = right(epsy, y);
+         Z = right(epsz, z);
+         break;
+   }
+}
+
+namespace mfem
+{
+
+struct KershawTransformation : VectorCoefficient
+{
+   double epsy, epsz;
+   int dim, s;
+   KershawTransformation(int dim_, double epsy_, double epsz_, int s_=0)
+      : VectorCoefficient(dim_), epsy(epsy_), epsz(epsz_), dim(dim_), s(s_) { }
+   using VectorCoefficient::Eval;
+   void Eval(Vector &V, ElementTransformation &T,
+             const IntegrationPoint &ip) override
+   {
+      double xyz[3];
+      Vector transip(xyz, 3);
+      T.Transform(ip, transip);
+      if (dim == 1)
+      {
+         V[0] = xyz[0]; // no transformation in 1D
+      }
+      else if (dim == 2)
+      {
+         double z=0, zt;
+         kershaw(epsy, epsz, s, xyz[0], xyz[1], z, V[0], V[1], zt);
+      }
+      else // dim == 3
+      {
+         kershaw(epsy, epsz, s, xyz[0], xyz[1], xyz[2], V[0], V[1], V[2]);
+      }
+   }
+};
+
+}
+
+#endif

tests/benchmarks/makefile

@@ -18,8 +18,8 @@ CONFIG_MK = $(MFEM_BUILD_DIR)/config/config.mk
 MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 
-SEQ_TESTS = bench_ceed bench_tmop bench_vector bench_virtuals
-PAR_TESTS = 
+SEQ_TESTS = bench_ceed bench_tmop bench_vector bench_virtuals bench_dgmassinv
+PAR_TESTS =
 ifeq ($(MFEM_USE_MPI),NO)
    TESTS = $(SEQ_TESTS)
 else
@@ -59,4 +59,3 @@ clean-build:
 
 clean-exec:
 	@rm -f refined.mesh sol.gf
-

tests/benchmarks/plot/plot_benchmarks.sh

@@ -0,0 +1,26 @@
+#!/usr/bin/env bash
+
+solver_types="Full_CG Local_CG_Lobatto Local_CG_Legendre Direct Direct_cuSolver Direct_cuBLAS"
+op_type="Setup Solve Setup_and_Solve"
+
+for stype in $solver_types
+do
+   for op in $op_type
+   do
+      bench_name="${stype}_${op//_/}"
+      csv_name="${bench_name,,}.csv"
+      tex_name="${bench_name,,}.tex"
+
+      title_name="${stype} ${op}"
+      title_name=${title_name//_/ }
+
+      cat << EOF > ${tex_name}
+\def\RUNTITLE{${title_name}}
+\def\DATA{${csv_name}}
+\input{template.tex}
+EOF
+
+      latexmk -pdflatex="pdflatex -halt-on-error" -synctex=1 -outdir=build ${tex_name}
+      cp build/${bench_name,,}.pdf ./
+   done
+done

tests/benchmarks/run_all_mass.sh

@@ -0,0 +1,20 @@
+#!/usr/bin/env bash
+
+solver_types="FULL_CG LOCAL_CG_LOBATTO LOCAL_CG_LEGENDRE DIRECT DIRECT_CUSOLVER DIRECT_CUBLAS"
+op_type="Setup Solve SetupAndSolve"
+
+GREEN='\033[0;32m'
+NC='\033[0m'
+
+for stype in $solver_types
+do
+   for op in $op_type
+   do
+      bench_name="${stype}_${op}"
+      bench_filter="${bench_name}_0_5/.*/.*"
+      output_name="${bench_name,,}.csv"
+      cmd="./bench_dgmassinv --benchmark_filter=${bench_filter} --benchmark_context=device=cuda --benchmark_out_format=csv --benchmark_out=${output_name}"
+      echo -e "${GREEN}${cmd}${NC}"
+      $cmd
+   done
+done

tests/unit/fem/test_dgmassinv.cpp

@@ -0,0 +1,101 @@
+// 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 "mfem.hpp"
+#include "unit_tests.hpp"
+
+using namespace mfem;
+
+TEST_CASE("DG Mass Inverse", "[CUDA]")
+{
+   auto mesh_filename = GENERATE(
+                           "../../data/star.mesh",
+                           "../../data/star-q3.mesh",
+                           "../../data/fichera.mesh"
+                        );
+   auto order = GENERATE(2, 3, 4, 5);
+   auto btype1 = GENERATE(BasisType::GaussLobatto, BasisType::GaussLegendre);
+
+   CAPTURE(mesh_filename, order, btype1);
+
+   Mesh mesh = Mesh::LoadFromFile(mesh_filename);
+   DG_FECollection fec(order, mesh.Dimension(), btype1);
+   FiniteElementSpace fes(&mesh, &fec);
+
+   BilinearForm m(&fes);
+   m.AddDomainIntegrator(new MassIntegrator);
+   m.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+   m.Assemble();
+
+   Array<int> empty;
+   OperatorJacobiSmoother jacobi(m, empty);
+
+   int n = fes.GetTrueVSize();
+   Vector B(n), X1(n), X2(n), X3(n);
+   B.Randomize(1);
+   X1.Randomize(2);
+   X2 = X1;
+   X3 = X1;
+
+   const double tol = 1e-8;
+
+   CGSolver cg;
+   cg.SetAbsTol(tol);
+   cg.SetRelTol(0.0);
+   cg.SetMaxIter(100);
+   cg.SetPrintLevel(IterativeSolver::PrintLevel().None());
+   cg.SetOperator(m);
+   cg.SetPreconditioner(jacobi);
+   cg.Mult(B, X1);
+
+   SECTION("Local CG")
+   {
+      auto btype2 = GENERATE(BasisType::GaussLobatto, BasisType::GaussLegendre);
+      CAPTURE(btype2);
+
+      DGMassInverse m_inv(fes, btype2);
+      m_inv.SetAbsTol(tol);
+      m_inv.SetRelTol(0.0);
+      m_inv.Mult(B, X2);
+
+      X2 -= X1;
+      REQUIRE(X2.Normlinf() == MFEM_Approx(0.0, 1e2*tol, 1e2*tol));
+   }
+
+   if (fes.GetFE(0)->GetDof() <= 64)
+   {
+      // Test direct solvers when the blocks are relatively small
+      SECTION("Direct")
+      {
+         DGMassInverse_Direct m_inv_direct(fes, BatchSolverMode::NATIVE);
+         m_inv_direct.Mult(B, X3);
+         X3 -= X1;
+         REQUIRE(X3.Normlinf() == MFEM_Approx(0.0, 1e2*tol, 1e2*tol));
+      }
+      if (Device::Allows(Backend::CUDA))
+      {
+         SECTION("Direct CuSolver")
+         {
+            DGMassInverse_Direct m_inv_direct(fes, BatchSolverMode::CUSOLVER);
+            m_inv_direct.Mult(B, X3);
+            X3 -= X1;
+            REQUIRE(X3.Normlinf() == MFEM_Approx(0.0, 1e2*tol, 1e2*tol));
+         }
+         SECTION("Direct CuBLAS")
+         {
+            DGMassInverse_Direct m_inv_direct(fes, BatchSolverMode::CUBLAS);
+            m_inv_direct.Mult(B, X3);
+            X3 -= X1;
+            REQUIRE(X3.Normlinf() == MFEM_Approx(0.0, 1e2*tol, 1e2*tol));
+         }
+      }
+   }
+}