Merge master in lor-gpu

Co-authored-by: Yohann Dudouit <dudouit1@llnl.gov>

.gitignore

@@ -18,6 +18,11 @@ CMakeFiles/
 # Backup files
 *~
 
+*.sqlite
+*.nsys-rep
+*.qdstrm
+*.csv
+
 # Default install location
 /mfem/
 
@@ -329,6 +334,7 @@ tests/benchmarks/bench_ceed
 tests/benchmarks/bench_tmop
 tests/benchmarks/bench_vector
 tests/benchmarks/bench_virtuals
+tests/benchmarks/bench_lor
 
 # Test script output
 tests/scripts/*.err

fem/bilinearform_ext.cpp

@@ -18,6 +18,8 @@
 #include "pgridfunc.hpp"
 #include "ceed/interface/util.hpp"
 
+#include "../general/nvtx.hpp"
+
 namespace mfem
 {
 
@@ -289,6 +291,10 @@ void PABilinearFormExtension::SetupRestrictionOperators(const L2FaceValues m)
 
 void PABilinearFormExtension::Assemble()
 {
+#undef MFEM_NVTX_COLOR
+#define MFEM_NVTX_COLOR NavyBlue
+   NVTX("HO Assemble");
+
    SetupRestrictionOperators(L2FaceValues::DoubleValued);
 
    Array<BilinearFormIntegrator*> &integrators = *a->GetDBFI();
@@ -383,6 +389,10 @@ void PABilinearFormExtension::FormLinearSystem(const Array<int> &ess_tdof_list,
 
 void PABilinearFormExtension::Mult(const Vector &x, Vector &y) const
 {
+#undef MFEM_NVTX_COLOR
+#define MFEM_NVTX_COLOR MediumSpringGreen
+   NVTX("HO Apply");
+
    Array<BilinearFormIntegrator*> &integrators = *a->GetDBFI();
 
    const int iSz = integrators.Size();

fem/bilininteg_diffusion_pa.cpp

@@ -1686,7 +1686,7 @@ static void PADiffusionApply(const int dim,
          case 0x77: return SmemPADiffusionApply2D<7,7,4>(NE,symm,B,G,D,X,Y);
          case 0x88: return SmemPADiffusionApply2D<8,8,2>(NE,symm,B,G,D,X,Y);
          case 0x99: return SmemPADiffusionApply2D<9,9,2>(NE,symm,B,G,D,X,Y);
-         default:   return PADiffusionApply2D(NE,symm,B,G,Bt,Gt,D,X,Y,D1D,Q1D);
+            // default:   return PADiffusionApply2D(NE,symm,B,G,Bt,Gt,D,X,Y,D1D,Q1D);
       }
    }
 
@@ -1704,7 +1704,14 @@ static void PADiffusionApply(const int dim,
          case 0x67: return SmemPADiffusionApply3D<6,7>(NE,symm,B,G,D,X,Y);
          case 0x78: return SmemPADiffusionApply3D<7,8>(NE,symm,B,G,D,X,Y);
          case 0x89: return SmemPADiffusionApply3D<8,9>(NE,symm,B,G,D,X,Y);
-         default:   return PADiffusionApply3D(NE,symm,B,G,Bt,Gt,D,X,Y,D1D,Q1D);
+         case 0x33: return SmemPADiffusionApply3D<3,3>(NE,symm,B,G,D,X,Y);
+         case 0x44: return SmemPADiffusionApply3D<4,4>(NE,symm,B,G,D,X,Y);
+         case 0x55: return SmemPADiffusionApply3D<5,5>(NE,symm,B,G,D,X,Y);
+         case 0x66: return SmemPADiffusionApply3D<6,6>(NE,symm,B,G,D,X,Y);
+         case 0x77: return SmemPADiffusionApply3D<7,7>(NE,symm,B,G,D,X,Y);
+         case 0x88: return SmemPADiffusionApply3D<8,8>(NE,symm,B,G,D,X,Y);
+         case 0x99: return SmemPADiffusionApply3D<9,9>(NE,symm,B,G,D,X,Y);
+            // default:   return PADiffusionApply3D(NE,symm,B,G,Bt,Gt,D,X,Y,D1D,Q1D);
       }
    }
    MFEM_ABORT("Unknown kernel: 0x"<<std::hex << id << std::dec);

fem/lor/lor.cpp

@@ -365,8 +365,11 @@ FiniteElementSpace &LORBase::GetFESpace() const
 
 void LORBase::AssembleSystem(BilinearForm &a_ho, const Array<int> &ess_dofs)
 {
-   A.Clear();
-   delete a;
+   if (a)
+   {
+      A.Clear();
+      delete a;
+   }
    if (BatchedLORAssembly::FormIsSupported(a_ho))
    {
       // Skip forming the space

fem/lor/lor_ads.cpp

@@ -13,6 +13,9 @@
 #include "../../general/forall.hpp"
 #include "../../fem/pbilinearform.hpp"
 
+#define MFEM_NVTX_COLOR DeepSkyBlue
+#include "../../general/nvtx.hpp"
+
 namespace mfem
 {
 
@@ -78,8 +81,10 @@ void BatchedLOR_ADS::Form3DFaceToEdge(Array<int> &face2edge)
    }
 }
 
-void BatchedLOR_ADS::FormCurlMatrix()
+void BatchedLOR_ADS::FormCurlMatrixLocal()
 {
+   NVTX("Discrete Curl");
+
    // The curl matrix maps from LOR edges to LOR faces. Given a quadrilateral
    // face (defined by its four edges) f_i = (e_j1, e_j2, e_j3, e_j4), the
    // matrix has nonzeros A(i, jk), so there are always exactly four nonzeros
@@ -87,10 +92,9 @@ void BatchedLOR_ADS::FormCurlMatrix()
    const int nface_dof = face_fes.GetNDofs();
    const int nedge_dof = edge_fes.GetNDofs();
 
-   SparseMatrix C_local;
    C_local.OverrideSize(nface_dof, nedge_dof);
-
-   C_local.GetMemoryI().New(nedge_dof+1, Device::GetDeviceMemoryType());
+   EnsureCapacity(C_local.GetMemoryI(), nedge_dof+1,
+                  Device::GetDeviceMemoryType());
    // Each row always has four nonzeros
    const int nnz = 4*nedge_dof;
    auto I = C_local.WriteI();
@@ -120,8 +124,8 @@ void BatchedLOR_ADS::FormCurlMatrix()
    const auto f2e = Reshape(face2edge.Read(), 4, nface_per_el);
 
    // Fill J and data
-   C_local.GetMemoryJ().New(nnz, Device::GetDeviceMemoryType());
-   C_local.GetMemoryData().New(nnz, Device::GetDeviceMemoryType());
+   EnsureCapacity(C_local.GetMemoryJ(), nnz, Device::GetDeviceMemoryType());
+   EnsureCapacity(C_local.GetMemoryData(), nnz, Device::GetDeviceMemoryType());
 
    auto J = C_local.WriteJ();
    auto V = C_local.WriteData();
@@ -146,6 +150,11 @@ void BatchedLOR_ADS::FormCurlMatrix()
          V[i*4 + k] = sgn*sgn_f*sgn_e;
       }
    });
+}
+
+void BatchedLOR_ADS::FormCurlMatrix()
+{
+   FormCurlMatrixLocal();
 
    // Create a block diagonal parallel matrix
    OperatorHandle C_diag(Operator::Hypre_ParCSR);
@@ -179,6 +188,8 @@ void BatchedLOR_ADS::FormCurlMatrix()
    }
    C->CopyRowStarts();
    C->CopyColStarts();
+
+   C_local.Clear();
 }
 
 HypreParMatrix *BatchedLOR_ADS::StealCurlMatrix()

fem/lor/lor_ads.hpp

@@ -36,7 +36,9 @@ protected:
    ND_FECollection edge_fec; ///< The associated Nedelec collection.
    ParFiniteElementSpace edge_fes; ///< The associated Nedelec space.
    BatchedLOR_AMS ams; ///< The associated AMS object.
-   HypreParMatrix *C; ///< The discrete curl matrix.
+   HypreParMatrix *C = nullptr; ///< The discrete curl matrix.
+
+   SparseMatrix C_local;
 
    /// Form the local elementwise discrete curl matrix.
    void Form3DFaceToEdge(Array<int> &face2edge);
@@ -64,6 +66,9 @@ public:
 
    /// Form the discrete curl matrix (not part of the public API).
    void FormCurlMatrix();
+
+   void FormCurlMatrixLocal();
+
    ~BatchedLOR_ADS();
 };
 

fem/lor/lor_ams.cpp

@@ -13,6 +13,9 @@
 #include "../../general/forall.hpp"
 #include "../../fem/pbilinearform.hpp"
 
+#define MFEM_NVTX_COLOR DeepSkyBlue
+#include "../../general/nvtx.hpp"
+
 namespace mfem
 {
 
@@ -136,8 +139,9 @@ void BatchedLOR_AMS::Form3DEdgeToVertex(Array<int> &edge2vert)
    }
 }
 
-void BatchedLOR_AMS::FormGradientMatrix()
+void BatchedLOR_AMS::FormGradientMatrixLocal()
 {
+   NVTX("Discrete Gradient");
    // The gradient matrix maps from LOR vertices to LOR edges. Given an edge
    // (defined by its two vertices) e_i = (v_j1, v_j2), the matrix has nonzeros
    // A(i, j1) = -1 and A(i, j2) = 1, so there are always exactly two nonzeros
@@ -145,10 +149,10 @@ void BatchedLOR_AMS::FormGradientMatrix()
    const int nedge_dof = edge_fes.GetNDofs();
    const int nvert_dof = vert_fes.GetNDofs();
 
-   SparseMatrix G_local;
    G_local.OverrideSize(nedge_dof, nvert_dof);
+   EnsureCapacity(G_local.GetMemoryI(), nedge_dof+1,
+                  Device::GetDeviceMemoryType());
 
-   G_local.GetMemoryI().New(nedge_dof+1, Device::GetDeviceMemoryType());
    // Each row always has two nonzeros
    const int nnz = 2*nedge_dof;
    auto I = G_local.WriteI();
@@ -180,8 +184,8 @@ void BatchedLOR_AMS::FormGradientMatrix()
    const auto e2v = Reshape(edge2vertex.Read(), 2, nedge_per_el);
 
    // Fill J and data
-   G_local.GetMemoryJ().New(nnz, Device::GetDeviceMemoryType());
-   G_local.GetMemoryData().New(nnz, Device::GetDeviceMemoryType());
+   EnsureCapacity(G_local.GetMemoryJ(), nnz, Device::GetDeviceMemoryType());
+   EnsureCapacity(G_local.GetMemoryData(), nnz, Device::GetDeviceMemoryType());
 
    auto J = G_local.WriteJ();
    auto V = G_local.WriteData();
@@ -204,6 +208,11 @@ void BatchedLOR_AMS::FormGradientMatrix()
       V[i*2 + 0] = -sgn;
       V[i*2 + 1] = sgn;
    });
+}
+
+void BatchedLOR_AMS::FormGradientMatrix()
+{
+   FormGradientMatrixLocal();
 
    // Create a block diagonal parallel matrix
    OperatorHandle G_diag(Operator::Hypre_ParCSR);
@@ -237,6 +246,8 @@ void BatchedLOR_AMS::FormGradientMatrix()
    }
    G->CopyRowStarts();
    G->CopyColStarts();
+
+   G_local.Clear();
 }
 
 template <typename T>
@@ -281,7 +292,7 @@ void BatchedLOR_AMS::FormCoordinateVectors(const Vector &X_vert)
    const MemoryClass mc = GetHypreMemoryClass();
    bool dev = (mc == MemoryClass::DEVICE);
 
-   xyz_tvec = new Vector(ntdofs*dim);
+   if (xyz_tvec == nullptr) { xyz_tvec = new Vector(ntdofs*dim); }
 
    auto xyz_tv = Reshape(HypreWrite(xyz_tvec->GetMemory()), ntdofs, dim);
    const auto xyz_e =
@@ -302,21 +313,24 @@ void BatchedLOR_AMS::FormCoordinateVectors(const Vector &X_vert)
    });
 
    // Make x, y, z HypreParVectors point to T-vector data
-   HYPRE_BigInt glob_size = vert_fes.GlobalTrueVSize();
-   HYPRE_BigInt *cols = vert_fes.GetTrueDofOffsets();
+   if (x == nullptr)
+   {
+      HYPRE_BigInt glob_size = vert_fes.GlobalTrueVSize();
+      HYPRE_BigInt *cols = vert_fes.GetTrueDofOffsets();
 
-   double *d_x_ptr = xyz_tv + 0*ntdofs;
-   x = new HypreParVector(vert_fes.GetComm(), glob_size, d_x_ptr, cols, dev);
-   double *d_y_ptr = xyz_tv + 1*ntdofs;
-   y = new HypreParVector(vert_fes.GetComm(), glob_size, d_y_ptr, cols, dev);
-   if (dim == 3)
-   {
-      double *d_z_ptr = xyz_tv + 2*ntdofs;
-      z = new HypreParVector(vert_fes.GetComm(), glob_size, d_z_ptr, cols, dev);
-   }
-   else
-   {
-      z = NULL;
+      double *d_x_ptr = xyz_tv + 0*ntdofs;
+      x = new HypreParVector(vert_fes.GetComm(), glob_size, d_x_ptr, cols, dev);
+      double *d_y_ptr = xyz_tv + 1*ntdofs;
+      y = new HypreParVector(vert_fes.GetComm(), glob_size, d_y_ptr, cols, dev);
+      if (dim == 3)
+      {
+         double *d_z_ptr = xyz_tv + 2*ntdofs;
+         z = new HypreParVector(vert_fes.GetComm(), glob_size, d_z_ptr, cols, dev);
+      }
+      else
+      {
+         z = NULL;
+      }
    }
 }
 

fem/lor/lor_ams.hpp

@@ -33,12 +33,14 @@ protected:
    const int order; ///< Polynomial degree.
    H1_FECollection vert_fec; ///< The corresponding H1 collection.
    ParFiniteElementSpace vert_fes; ///< The corresponding H1 space.
-   Vector *xyz_tvec; ///< Mesh vertex coordinates in true-vector format.
-   HypreParMatrix *G; ///< Discrete gradient matrix.
+   Vector *xyz_tvec = nullptr; ///< Mesh vertex coordinates in true-vector format.
+   HypreParMatrix *G = nullptr; ///< Discrete gradient matrix.
+
+   SparseMatrix G_local;
 
    /// @name Mesh coordinate vectors in HypreParVector format
    ///@{
-   HypreParVector *x, *y, *z;
+   HypreParVector *x = nullptr, *y = nullptr, *z = nullptr;
    ///@}
 
    /// @name Construct the local (elementwise) discrete gradient
@@ -89,6 +91,9 @@ public:
 
    /// Construct the discrete gradient matrix (not part of the public API).
    void FormGradientMatrix();
+
+   void FormGradientMatrixLocal();
+
    ~BatchedLOR_AMS();
 };
 

fem/lor/lor_batched.cpp

@@ -15,6 +15,8 @@
 #include <climits>
 #include "../pbilinearform.hpp"
 
+#include "../../general/nvtx.hpp"
+
 // Specializations
 #include "lor_h1.hpp"
 #include "lor_nd.hpp"
@@ -70,8 +72,13 @@ bool BatchedLORAssembly::FormIsSupported(BilinearForm &a)
 }
 
 void BatchedLORAssembly::FormLORVertexCoordinates(FiniteElementSpace &fes_ho,
-                                                  Vector &X_vert)
+                                                  Vector &X_vert,
+                                                  Vector *evec)
 {
+#undef MFEM_NVTX_COLOR
+#define MFEM_NVTX_COLOR DeepSkyBlue
+   NVTX("LOR Coordinates");
+
    Mesh &mesh_ho = *fes_ho.GetMesh();
    mesh_ho.EnsureNodes();
 
@@ -87,9 +94,22 @@ void BatchedLORAssembly::FormLORVertexCoordinates(FiniteElementSpace &fes_ho,
    const Operator *nodal_restriction =
       nodal_fes->GetElementRestriction(ElementDofOrdering::LEXICOGRAPHIC);
 
+   Vector *tmp_evec = nullptr;
+   Vector *nodal_evec;
+
+   if (evec)
+   {
+      nodal_evec = evec;
+      nodal_evec->SetSize(nodal_restriction->Height());
+   }
+   else
+   {
+      tmp_evec = new Vector(nodal_restriction->Height());
+      nodal_evec = tmp_evec;
+   }
+
    // Map from nodal L-vector to E-vector
-   Vector nodal_evec(nodal_restriction->Height());
-   nodal_restriction->Mult(*nodal_gf, nodal_evec);
+   nodal_restriction->Mult(*nodal_gf, *nodal_evec);
 
    IntegrationRule ir = GetCollocatedIntRule(fes_ho);
 
@@ -98,7 +118,9 @@ void BatchedLORAssembly::FormLORVertexCoordinates(FiniteElementSpace &fes_ho,
    const QuadratureInterpolator *quad_interp =
       nodal_fes->GetQuadratureInterpolator(ir);
    quad_interp->SetOutputLayout(QVectorLayout::byVDIM);
-   quad_interp->Values(nodal_evec, X_vert);
+   quad_interp->Values(*nodal_evec, X_vert);
+
+   delete tmp_evec;
 }
 
 // The following two functions (GetMinElt and GetAndIncrementNnzIndex) are
@@ -372,11 +394,11 @@ void BatchedLORAssembly::SparseIJToCSR(OperatorHandle &A) const
 
    A_mat->OverrideSize(nvdof, nvdof);
 
-   A_mat->GetMemoryI().New(nvdof+1, Device::GetDeviceMemoryType());
+   EnsureCapacity(A_mat->GetMemoryI(), nvdof+1, Device::GetDeviceMemoryType());
    int nnz = FillI(*A_mat);
 
-   A_mat->GetMemoryJ().New(nnz, Device::GetDeviceMemoryType());
-   A_mat->GetMemoryData().New(nnz, Device::GetDeviceMemoryType());
+   EnsureCapacity(A_mat->GetMemoryJ(), nnz, Device::GetDeviceMemoryType());
+   EnsureCapacity(A_mat->GetMemoryData(), nnz, Device::GetDeviceMemoryType());
    FillJAndData(*A_mat);
 }
 
@@ -455,7 +477,6 @@ void BatchedLORAssembly::ParAssemble(
    BilinearForm &a, const Array<int> &ess_dofs, OperatorHandle &A)
 {
    // Assemble the system matrix local to this partition
-   OperatorHandle A_local;
    AssembleWithoutBC(a, A_local);
 
    ParBilinearForm *pa =
@@ -471,6 +492,9 @@ void BatchedLORAssembly::ParAssemble(
 void BatchedLORAssembly::Assemble(
    BilinearForm &a, const Array<int> ess_dofs, OperatorHandle &A)
 {
+#undef MFEM_NVTX_COLOR
+#define MFEM_NVTX_COLOR NavyBlue
+   NVTX("LOR Assemble");
 #ifdef MFEM_USE_MPI
    if (dynamic_cast<ParFiniteElementSpace*>(&fes_ho))
    {

fem/lor/lor_batched.hpp

@@ -53,6 +53,8 @@ protected:
    /// nonzero).
    Array<int> sparse_mapping;
 
+   OperatorHandle A_local; // Cache this
+
 public:
    /// Construct the batched assembly object corresponding to @a fes_ho_.
    BatchedLORAssembly(FiniteElementSpace &fes_ho_);
@@ -68,12 +70,12 @@ public:
 
    /// Compute the vertices of the LOR mesh and place the result in @a X_vert.
    static void FormLORVertexCoordinates(FiniteElementSpace &fes_ho,
-                                        Vector &X_vert);
+                                        Vector &X_vert,
+                                        Vector *evec = nullptr);
 
    /// Return the vertices of the LOR mesh in E-vector format
    const Vector &GetLORVertexCoordinates() { return X_vert; }
 
-protected:
    /// After assembling the "sparse IJ" format, convert it to CSR.
    void SparseIJToCSR(OperatorHandle &A) const;
 
@@ -117,12 +119,12 @@ public:
 /// If the capacity of @a mem is not large enough, delete it and allocate new
 /// memory with size @a capacity.
 template <typename T>
-void EnsureCapacity(Memory<T> &mem, int capacity)
+void EnsureCapacity(Memory<T> &mem, int capacity, MemoryType mt)
 {
    if (mem.Capacity() < capacity)
    {
       mem.Delete();
-      mem.New(capacity, mem.GetMemoryType());
+      mem.New(capacity, mt);
    }
 }
 

fem/lor/lor_h1.cpp

@@ -205,6 +205,49 @@ void BatchedLOR_H1::Assemble2D()
    }
 }
 
+template<int ORDER>
+static void SparseMapping3D(Array<int> &sparse_mapping)
+{
+   static constexpr int nnz_per_row = 27;
+   static constexpr int nd1d = ORDER + 1;
+   static constexpr int ndof_per_el = nd1d*nd1d*nd1d;
+
+   sparse_mapping.SetSize(nnz_per_row*ndof_per_el);
+   sparse_mapping = -1;
+   auto map = Reshape(sparse_mapping.HostReadWrite(), nnz_per_row, ndof_per_el);
+
+   for (int iz=0; iz<nd1d; ++iz)
+   {
+      const int jz_begin = (iz > 0) ? iz - 1 : 0;
+      const int jz_end = (iz < ORDER) ? iz + 1 : ORDER;
+      for (int iy=0; iy<nd1d; ++iy)
+      {
+         const int jy_begin = (iy > 0) ? iy - 1 : 0;
+         const int jy_end = (iy < ORDER) ? iy + 1 : ORDER;
+         for (int ix=0; ix<nd1d; ++ix)
+         {
+            const int jx_begin = (ix > 0) ? ix - 1 : 0;
+            const int jx_end = (ix < ORDER) ? ix + 1 : ORDER;
+
+            const int ii_el = ix + nd1d*(iy + nd1d*iz);
+
+            for (int jz=jz_begin; jz<=jz_end; ++jz)
+            {
+               for (int jy=jy_begin; jy<=jy_end; ++jy)
+               {
+                  for (int jx=jx_begin; jx<=jx_end; ++jx)
+                  {
+                     const int jj_off = (jx-ix+1) + 3*(jy-iy+1) + 9*(jz-iz+1);
+                     const int jj_el = jx + nd1d*(jy + nd1d*jz);
+                     map(jj_off, ii_el) = jj_el;
+                  }
+               }
+            }
+         }
+      }
+   }
+}
+
 template <int ORDER>
 void BatchedLOR_H1::Assemble3D()
 {
@@ -330,21 +373,16 @@ void BatchedLOR_H1::Assemble3D()
                   }
                }
 
-               //MFEM_UNROLL(2)
                for (int iqx=0; iqx<2; ++iqx)
                {
-                  //MFEM_UNROLL(2)
                   for (int jz=0; jz<2; ++jz)
                   {
                      // Note loop starts at iz=jz here, taking advantage of
                      // symmetries.
-                     //MFEM_UNROLL(2)
                      for (int iz=jz; iz<2; ++iz)
                      {
-                        //MFEM_UNROLL(2)
                         for (int iqy=0; iqy<2; ++iqy)
                         {
-                           //MFEM_UNROLL(2)
                            for (int iqz=0; iqz<2; ++iqz)
                            {
                               const double mq = const_mq ? MQ(0,0,0,0) : MQ(kx+iqx, ky+iqy, kz+iqz, iel_ho);
@@ -379,10 +417,8 @@ void BatchedLOR_H1::Assemble3D()
                               double wdetJ = Q(6,iqz,iqy,iqx);
                               mass_A(iqy,iz,jz,iqx) += mq*wdetJ*biz*bjz;
                            }
-                           //MFEM_UNROLL(2)
                            for (int jy=0; jy<2; ++jy)
                            {
-                              //MFEM_UNROLL(2)
                               for (int iy=0; iy<2; ++iy)
                               {
                                  const double biy = (iy == iqy) ? 1.0 : 0.0;
@@ -405,16 +441,12 @@ void BatchedLOR_H1::Assemble3D()
                               }
                            }
                         }
-                        //MFEM_UNROLL(2)
                         for (int jy=0; jy<2; ++jy)
                         {
-                           //MFEM_UNROLL(2)
                            for (int jx=0; jx<2; ++jx)
                            {
-                              //MFEM_UNROLL(2)
                               for (int iy=0; iy<2; ++iy)
                               {
-                                 //MFEM_UNROLL(2)
                                  for (int ix=0; ix<2; ++ix)
                                  {
                                     const double bix = (ix == iqx) ? 1.0 : 0.0;
@@ -482,40 +514,176 @@ void BatchedLOR_H1::Assemble3D()
          }
       }
    });
+   SparseMapping3D<ORDER>(sparse_mapping);
+}
 
-   sparse_mapping.SetSize(nnz_per_row*ndof_per_el);
-   sparse_mapping = -1;
-   auto map = Reshape(sparse_mapping.HostReadWrite(), nnz_per_row, ndof_per_el);
-   for (int iz=0; iz<nd1d; ++iz)
+template <>
+void BatchedLOR_H1::Assemble3D<1>()
+{
+   static constexpr int nv = 8;
+   static constexpr int nd1d = 2;
+   static constexpr int ndof_per_el = 8;
+   static constexpr int nnz_per_row = 27;
+   static constexpr int sz_local_mat = nv*nv;
+
+   const int nel_ho = fes_ho.GetNE();
+
+   const bool const_mq = c1.Size() == 1;
+   const auto MQ = const_mq
+                   ? Reshape(c1.Read(), 1, 1, 1, 1)
+                   : Reshape(c1.Read(), nd1d, nd1d, nd1d, nel_ho);
+   const bool const_dq = c2.Size() == 1;
+   const auto DQ = const_dq
+                   ? Reshape(c2.Read(), 1, 1, 1, 1)
+                   : Reshape(c2.Read(), nd1d, nd1d, nd1d, nel_ho);
+
+   sparse_ij.SetSize(nel_ho*ndof_per_el*nnz_per_row);
+   auto V = Reshape(sparse_ij.Write(), nnz_per_row, nd1d, nd1d, nd1d, nel_ho);
+
+   const auto X = X_vert.Read();
+
+   MFEM_FORALL_3D(iel_ho, nel_ho, 8, 4, 1,
    {
-      const int jz_begin = (iz > 0) ? iz - 1 : 0;
-      const int jz_end = (iz < ORDER) ? iz + 1 : ORDER;
-      for (int iy=0; iy<nd1d; ++iy)
+      static constexpr int e[8] = {0,1,3,2,4,5,7,6};
+      MFEM_SHARED double vx[8], vy[8], vz[8];
+      const int tidz = MFEM_THREAD_ID(z);
+
+      MFEM_SHARED double local_mat_[sz_local_mat];
+      DeviceTensor<4> local_mat(local_mat_, 2,2,2, nv);
+
+      if (tidz == 0)
       {
-         const int jy_begin = (iy > 0) ? iy - 1 : 0;
-         const int jy_end = (iy < ORDER) ? iy + 1 : ORDER;
-         for (int ix=0; ix<nd1d; ++ix)
+         MFEM_FOREACH_THREAD(xyz,x,8)
          {
-            const int jx_begin = (ix > 0) ? ix - 1 : 0;
-            const int jx_end = (ix < ORDER) ? ix + 1 : ORDER;
-
-            const int ii_el = ix + nd1d*(iy + nd1d*iz);
-
-            for (int jz=jz_begin; jz<=jz_end; ++jz)
+            const int z = xyz%2, y = (xyz/2)%2, x = xyz/2/2;
+            MFEM_FOREACH_THREAD(j,y,nnz_per_row)
             {
-               for (int jy=jy_begin; jy<=jy_end; ++jy)
+               if (j < 8) { local_mat(z,y,x,j) = 0.0; }
+               V(j,x,y,z,iel_ho) = 0.0;
+               if (j == 0)
                {
-                  for (int jx=jx_begin; jx<=jx_end; ++jx)
-                  {
-                     const int jj_off = (jx-ix+1) + 3*(jy-iy+1) + 9*(jz-iz+1);
-                     const int jj_el = jx + nd1d*(jy + nd1d*jz);
-                     map(jj_off, ii_el) = jj_el;
-                  }
+                  const int i = x + 2*y + 4*z;
+                  const int ei = 3*(e[i] + 8*iel_ho);
+                  vx[i] = X[ei + 0];
+                  vy[i] = X[ei + 1];
+                  vz[i] = X[ei + 2];
                }
             }
          }
       }
-   }
+      MFEM_SYNC_THREAD;
+
+      MFEM_FOREACH_THREAD(xyz,x,8)
+      {
+         const int qz = xyz%2, qy = (xyz/2)%2, qx = xyz/2/2;
+         static constexpr double w = 1.0/8.0;
+
+         double J_[3*3];
+         DeviceTensor<2> J(J_, 3,3);
+         Jacobian3D(qx,qy,qz, vx,vy,vz, J);
+         const double detJ = Det3D(J);
+         const double w_detJ = w/detJ;
+
+         // adj(J)
+         double A_[3*3];
+         DeviceTensor<2> A(A_, 3, 3);
+         Adjugate3D(J, A);
+
+         const double J11 = w_detJ*(A(0,0)*A(0,0)+A(0,1)*A(0,1)+A(0,2)*A(0,2)); // 1,1
+         const double J21 = w_detJ*(A(0,0)*A(1,0)+A(0,1)*A(1,1)+A(0,2)*A(1,2)); // 2,1
+         const double J31 = w_detJ*(A(0,0)*A(2,0)+A(0,1)*A(2,1)+A(0,2)*A(2,2)); // 3,1
+         const double J12 = J21;
+         const double J22 = w_detJ*(A(1,0)*A(1,0)+A(1,1)*A(1,1)+A(1,2)*A(1,2)); // 2,2
+         const double J32 = w_detJ*(A(1,0)*A(2,0)+A(1,1)*A(2,1)+A(1,2)*A(2,2)); // 3,2
+         const double J13 = J31;
+         const double J23 = J32;
+         const double J33 = w_detJ*(A(2,0)*A(2,0)+A(2,1)*A(2,1)+A(2,2)*A(2,2)); // 3,3
+         const double wdetJ =  w*detJ;
+
+         const double mq = const_mq ? MQ(0,0,0,0) : MQ(qx,qy,qz, iel_ho);
+         const double dq = const_dq ? DQ(0,0,0,0) : DQ(qx,qy,qz, iel_ho);
+
+         MFEM_FOREACH_THREAD(xyz,y,8)
+         {
+            const int jz = xyz%2, jy = (xyz/2)%2, jx = xyz/2/2;
+            const double bjz = (jz == qz) ? 1.0 : 0.0;
+            const double gjz = (jz == 0) ? -1.0 : 1.0;
+            const double bjy = (jy == qy) ? 1.0 : 0.0;
+            const double gjy = (jy == 0) ? -1.0 : 1.0;
+            const double bjx = (jx == qx) ? 1.0 : 0.0;
+            const double gjx = (jx == 0) ? -1.0 : 1.0;
+            const double djx = gjx*bjy*bjz;
+            const double djy = bjx*gjy*bjz;
+            const double djz = bjx*bjy*gjz;
+            const int jj_loc = jx + 2*jy + 4*jz;
+
+            MFEM_FOREACH_THREAD(xyz,z,8)
+            {
+               const int iz = xyz%2, iy = (xyz/2)%2, ix = xyz/2/2;
+               const double biz = (iz == qz) ? 1.0 : 0.0;
+               const double giz = (iz == 0) ? -1.0 : 1.0;
+               const double biy = (iy == qy) ? 1.0 : 0.0;
+               const double giy = (iy == 0) ? -1.0 : 1.0;
+               const double bix = (ix == qx) ? 1.0 : 0.0;
+               const double gix = (ix == 0) ? -1.0 : 1.0;
+               const double dix = gix*biy*biz;
+               const double diy = bix*giy*biz;
+               const double diz = bix*biy*giz;
+               const int ii_loc = ix + 2*iy + 4*iz;
+
+               // Only store the lower-triangular part of
+               // the matrix (by symmetry).
+               if (jj_loc > ii_loc) { continue; }
+
+               double grad_grad = 0.0;
+               grad_grad += dix*djx*J11;
+               grad_grad += diy*djx*J12;
+               grad_grad += diz*djx*J13;
+
+               grad_grad += dix*djy*J21;
+               grad_grad += diy*djy*J22;
+               grad_grad += diz*djy*J23;
+
+               grad_grad += dix*djz*J31;
+               grad_grad += diy*djz*J32;
+               grad_grad += diz*djz*J33;
+
+               const double basis_basis = wdetJ*bix*biy*biz*bjx*bjy*bjz;
+               const double value = dq*grad_grad + mq*basis_basis;
+               AtomicAdd(local_mat(iz,iy,ix, jj_loc), value);
+            } // i
+         } // j
+      } // q
+      MFEM_SYNC_THREAD;
+
+      // Assemble the local matrix into the macro-element sparse matrix
+      // in a format similar to coordinate format.
+      // The (I,J) arrays are implicit (not stored explicitly).
+      if (tidz == 0)
+      {
+         MFEM_FOREACH_THREAD(xyz,x,8)
+         {
+            const int iz = xyz%2, iy = (xyz/2)%2, ix = xyz/2/2;
+            const int ii_loc = ix + 2*iy + 4*iz;
+            MFEM_FOREACH_THREAD(jj_loc,y,8)
+            {
+               const int jx = jj_loc%2, jy = (jj_loc/2)%2, jz = jj_loc/2/2;
+               const int jj_off = (jx-ix+1) + 3*(jy-iy+1) + 9*(jz-iz+1);
+
+               if (jj_loc <= ii_loc)
+               {
+                  AtomicAdd(V(jj_off, ix,iy,iz, iel_ho), local_mat(iz,iy,ix, jj_loc));
+               }
+               else
+               {
+                  AtomicAdd(V(jj_off, ix,iy,iz, iel_ho), local_mat(jz,jy,jx, ii_loc));
+               }
+            }
+         }
+      }
+      MFEM_SYNC_THREAD;
+   });
+   SparseMapping3D<1>(sparse_mapping);
 }
 
 // Explicit template instantiations
@@ -528,7 +696,7 @@ template void BatchedLOR_H1::Assemble2D<6>();
 template void BatchedLOR_H1::Assemble2D<7>();
 template void BatchedLOR_H1::Assemble2D<8>();
 
-template void BatchedLOR_H1::Assemble3D<1>();
+//template void BatchedLOR_H1::Assemble3D<1>(); // explicitly specialized
 template void BatchedLOR_H1::Assemble3D<2>();
 template void BatchedLOR_H1::Assemble3D<3>();
 template void BatchedLOR_H1::Assemble3D<4>();

fem/lor/lor_nd.cpp

@@ -216,6 +216,90 @@ void BatchedLOR_ND::Assemble2D()
    }
 }
 
+
+template<int ORDER>
+static void SparseMapping3D(Array<int> &sparse_mapping)
+{
+   static constexpr int nnz_per_row = 33;
+   static constexpr int dim = 3;
+   static constexpr int o = ORDER;
+   static constexpr int op1 = ORDER + 1;
+   static constexpr int ndof_per_el = dim*o*op1*op1;
+
+   sparse_mapping.SetSize(nnz_per_row*ndof_per_el);
+   sparse_mapping = -1;
+   auto map = Reshape(sparse_mapping.HostReadWrite(), nnz_per_row, ndof_per_el);
+
+   for (int ci=0; ci<dim; ++ci)
+   {
+      const int i_off = ci*o*op1*op1;
+      const int id0 = ci;
+      const int id1 = (ci+1)%3;
+      const int id2 = (ci+2)%3;
+
+      const int nxi = (ci == 0) ? o : op1;
+      const int nyi = (ci == 1) ? o : op1;
+
+      for (int i0=0; i0<o; ++i0)
+      {
+         for (int i1=0; i1<op1; ++i1)
+         {
+            for (int i2=0; i2<op1; ++i2)
+            {
+               int ii_lex[3];
+               ii_lex[id0] = i0;
+               ii_lex[id1] = i1;
+               ii_lex[id2] = i2;
+               const int ii_el = i_off + ii_lex[0] + ii_lex[1]*nxi + ii_lex[2]*nxi*nyi;
+
+               for (int cj_rel=0; cj_rel<dim; ++cj_rel)
+               {
+                  const int cj = (ci + cj_rel) % 3;
+                  const int j_off = cj*o*op1*op1;
+
+                  const int nxj = (cj == 0) ? o : op1;
+                  const int nyj = (cj == 1) ? o : op1;
+
+                  const int j0_begin = i0;
+                  const int j0_end = (cj_rel == 0) ? i0 : i0 + 1;
+                  const int j1_begin = (i1 > 0) ? i1-1 : i1;
+                  const int j1_end = (cj_rel == 1)
+                                     ? ((i1 < o) ? i1 : i1-1)
+                                     : ((i1 < o) ? i1+1 : i1);
+                  const int j2_begin = (i2 > 0) ? i2-1 : i2;
+                  const int j2_end = (cj_rel == 2)
+                                     ? ((i2 < o) ? i2 : i2-1)
+                                     : ((i2 < o) ? i2+1 : i2);
+
+                  for (int j0=j0_begin; j0<=j0_end; ++j0)
+                  {
+                     const int d0 = j0 - i0;
+                     for (int j1=j1_begin; j1<=j1_end; ++j1)
+                     {
+                        const int d1 = j1 - i1 + 1;
+                        for (int j2=j2_begin; j2<=j2_end; ++j2)
+                        {
+                           const int d2 = j2 - i2 + 1;
+                           int jj_lex[3];
+                           jj_lex[id0] = j0;
+                           jj_lex[id1] = j1;
+                           jj_lex[id2] = j2;
+                           const int jj_el = j_off + jj_lex[0] + jj_lex[1]*nxj + jj_lex[2]*nxj*nyj;
+                           int jj_off;
+                           if (cj_rel == 0) { jj_off = d1 + 3*d2; }
+                           else if (cj_rel == 1) { jj_off = 9 + d0 + 2*d1 + 4*d2; }
+                           else /* if (cj_rel == 2) */ { jj_off = 21 + d0 + 2*d1 + 6*d2; }
+                           map(jj_off, ii_el) = jj_el;
+                        }
+                     }
+                  }
+               }
+            }
+         }
+      }
+   }
+}
+
 template <int ORDER>
 void BatchedLOR_ND::Assemble3D()
 {
@@ -244,7 +328,7 @@ void BatchedLOR_ND::Assemble3D()
    sparse_ij.SetSize(nnz_per_row*ndof_per_el*nel_ho);
    auto V = Reshape(sparse_ij.Write(), nnz_per_row, o*op1*op1, dim, nel_ho);
 
-   auto X = X_vert.Read();
+   const auto X = X_vert.Read();
 
    // Last thread dimension is lowered to avoid "too many resources" error
    MFEM_FORALL_3D(iel_ho, nel_ho, ORDER, ORDER, (ORDER>6)?4:ORDER,
@@ -488,78 +572,231 @@ void BatchedLOR_ND::Assemble3D()
          }
       }
    });
+   SparseMapping3D<ORDER>(sparse_mapping);
+}
 
-   sparse_mapping.SetSize(nnz_per_row*ndof_per_el);
-   sparse_mapping = -1;
-   auto map = Reshape(sparse_mapping.HostReadWrite(), nnz_per_row, ndof_per_el);
-   for (int ci=0; ci<dim; ++ci)
+template<>
+void BatchedLOR_ND::Assemble3D<1>()
+{
+   static constexpr int ne = 12; // number of edges in hexahedron
+   static constexpr int dim = 3;
+   static constexpr int o = 1;
+   static constexpr int op1 = 2;
+   static constexpr int ndof_per_el = dim*o*op1*op1;
+   static constexpr int nnz_per_row = 33;
+   static constexpr int sz_local_mat = ne*ne;
+
+   const int nel_ho = fes_ho.GetNE();
+   const bool const_mq = c1.Size() == 1;
+   const auto MQ = const_mq
+                   ? Reshape(c1.Read(), 1, 1, 1, 1)
+                   : Reshape(c1.Read(), op1, op1, op1, nel_ho);
+   const bool const_dq = c2.Size() == 1;
+   const auto DQ = const_dq
+                   ? Reshape(c2.Read(), 1, 1, 1, 1)
+                   : Reshape(c2.Read(), op1, op1, op1, nel_ho);
+
+   sparse_ij.SetSize(nnz_per_row*ndof_per_el*nel_ho);
+   auto V = Reshape(sparse_ij.Write(), nnz_per_row, o*op1*op1, dim, nel_ho);
+
+   const auto X = X_vert.Read();
+
+   MFEM_FORALL_3D(iel_ho, nel_ho, 8, 1, 4,
    {
-      const int i_off = ci*o*op1*op1;
-      const int id0 = ci;
-      const int id1 = (ci+1)%3;
-      const int id2 = (ci+2)%3;
-
-      const int nxi = (ci == 0) ? o : op1;
-      const int nyi = (ci == 1) ? o : op1;
-
-      for (int i0=0; i0<o; ++i0)
+      MFEM_FOREACH_THREAD(iz,z,o) // 1
       {
-         for (int i1=0; i1<op1; ++i1)
+         MFEM_FOREACH_THREAD(iy,y,op1) // 2
          {
-            for (int i2=0; i2<op1; ++i2)
+            MFEM_FOREACH_THREAD(ix,x,op1) // 2
             {
-               int ii_lex[3];
-               ii_lex[id0] = i0;
-               ii_lex[id1] = i1;
-               ii_lex[id2] = i2;
-               const int ii_el = i_off + ii_lex[0] + ii_lex[1]*nxi + ii_lex[2]*nxi*nyi;
-
-               for (int cj_rel=0; cj_rel<dim; ++cj_rel)
+               for (int c=0; c<dim; ++c)
                {
-                  const int cj = (ci + cj_rel) % 3;
-                  const int j_off = cj*o*op1*op1;
-
-                  const int nxj = (cj == 0) ? o : op1;
-                  const int nyj = (cj == 1) ? o : op1;
-
-                  const int j0_begin = i0;
-                  const int j0_end = (cj_rel == 0) ? i0 : i0 + 1;
-                  const int j1_begin = (i1 > 0) ? i1-1 : i1;
-                  const int j1_end = (cj_rel == 1)
-                                     ? ((i1 < o) ? i1 : i1-1)
-                                     : ((i1 < o) ? i1+1 : i1);
-                  const int j2_begin = (i2 > 0) ? i2-1 : i2;
-                  const int j2_end = (cj_rel == 2)
-                                     ? ((i2 < o) ? i2 : i2-1)
-                                     : ((i2 < o) ? i2+1 : i2);
-
-                  for (int j0=j0_begin; j0<=j0_end; ++j0)
+                  for (int j=0; j<nnz_per_row; ++j)
                   {
-                     const int d0 = j0 - i0;
-                     for (int j1=j1_begin; j1<=j1_end; ++j1)
-                     {
-                        const int d1 = j1 - i1 + 1;
-                        for (int j2=j2_begin; j2<=j2_end; ++j2)
-                        {
-                           const int d2 = j2 - i2 + 1;
-                           int jj_lex[3];
-                           jj_lex[id0] = j0;
-                           jj_lex[id1] = j1;
-                           jj_lex[id2] = j2;
-                           const int jj_el = j_off + jj_lex[0] + jj_lex[1]*nxj + jj_lex[2]*nxj*nyj;
-                           int jj_off;
-                           if (cj_rel == 0) { jj_off = d1 + 3*d2; }
-                           else if (cj_rel == 1) { jj_off = 9 + d0 + 2*d1 + 4*d2; }
-                           else /* if (cj_rel == 2) */ { jj_off = 21 + d0 + 2*d1 + 6*d2; }
-                           map(jj_off, ii_el) = jj_el;
-                        }
-                     }
+                     V(j,ix+iy*op1+iz*op1*op1,c,iel_ho) = 0.0;
                   }
                }
             }
          }
       }
-   }
+      MFEM_SYNC_THREAD;
+
+      MFEM_SHARED double local_mat_[sz_local_mat];
+      DeviceTensor<4> local_mat(local_mat_, 3,4, 3,4);
+      /// should be optimized
+      for (int i=0; i<sz_local_mat; ++i) { local_mat[i] = 0.0; }
+
+      MFEM_SHARED double vx[8], vy[8], vz[8];
+      /// should be optimized
+      LORVertexCoordinates3D<1>(X, iel_ho, 0,0,0, vx,vy,vz);
+      MFEM_SYNC_THREAD;
+
+      MFEM_FOREACH_THREAD(xyz,x,8)
+      {
+         const int qz = xyz%2, qy = (xyz/2)%2, qx = xyz/2/2;
+         static constexpr double w = 1.0/8.0;
+
+         double J_[3*3];
+         DeviceTensor<2> J(J_, 3,3);
+         Jacobian3D(qx,qy,qz, vx,vy,vz, J);
+         const double detJ = Det3D(J);
+         const double w_detJ = w/detJ;
+
+         // adj(J)
+         double A_[3*3];
+         DeviceTensor<2> A(A_, 3,3);
+         Adjugate3D(J, A);
+
+         const double Q0 = w_detJ*(A(0,0)*A(0,0)+A(0,1)*A(0,1)+A(0,2)*A(0,2)); // 1,1
+         const double Q1 = w_detJ*(A(0,0)*A(1,0)+A(0,1)*A(1,1)+A(0,2)*A(1,2)); // 2,1
+         const double Q2 = w_detJ*(A(0,0)*A(2,0)+A(0,1)*A(2,1)+A(0,2)*A(2,2)); // 3,1
+         const double Q3 = w_detJ*(A(1,0)*A(1,0)+A(1,1)*A(1,1)+A(1,2)*A(1,2)); // 2,2
+         const double Q4 = w_detJ*(A(1,0)*A(2,0)+A(1,1)*A(2,1)+A(1,2)*A(2,2)); // 3,2
+         const double Q5 = w_detJ*(A(2,0)*A(2,0)+A(2,1)*A(2,1)+A(2,2)*A(2,2)); // 3,3
+
+         // w J^T J / det(J)
+         const double Q6  = w_detJ*(J(0,0)*J(0,0)+J(1,0)*J(1,0)+J(2,0)*J(2,0)); // 1,1
+         const double Q7  = w_detJ*(J(0,0)*J(0,1)+J(1,0)*J(1,1)+J(2,0)*J(2,1)); // 2,1
+         const double Q8  = w_detJ*(J(0,0)*J(0,2)+J(1,0)*J(1,2)+J(2,0)*J(2,2)); // 3,1
+         const double Q9  = w_detJ*(J(0,1)*J(0,1)+J(1,1)*J(1,1)+J(2,1)*J(2,1)); // 2,2
+         const double Q10 = w_detJ*(J(0,1)*J(0,2)+J(1,1)*J(1,2)+J(2,1)*J(2,2)); // 3,2
+         const double Q11 = w_detJ*(J(0,2)*J(0,2)+J(1,2)*J(1,2)+J(2,2)*J(2,2)); // 3,3
+
+         const double mq = const_mq ? MQ(0,0,0,0) : MQ(qx,qy,qz, iel_ho);
+         const double dq = const_dq ? DQ(0,0,0,0) : DQ(qx,qy,qz, iel_ho);
+
+         // Loop over x,y,z components. 0 => x, 1 => y, 2 => z
+         MFEM_FOREACH_THREAD(cj,y,dim)
+         {
+            const double jq1 = (cj == 0) ? qy : ((cj == 1) ? qz : qx);
+            const double jq2 = (cj == 0) ? qz : ((cj == 1) ? qx : qy);
+
+            const int jd_0 = cj;
+            const int jd_1 = (cj + 1)%3;
+            const int jd_2 = (cj + 2)%3;
+
+            MFEM_FOREACH_THREAD(bj,z,4) // 4 edges in each dim
+            {
+               const int bj1 = bj%2;
+               const int bj2 = bj/2;
+
+               double curl_j[3];
+               curl_j[jd_0] = 0.0;
+               curl_j[jd_1] = ((bj1 == 0) ? jq1 - 1 : -jq1)*((bj2 == 0) ? 1 : -1);
+               curl_j[jd_2] = ((bj2 == 0) ? 1 - jq2 : jq2)*((bj1 == 0) ? 1 : -1);
+
+               double basis_j[3];
+               basis_j[jd_0] = ((bj1 == 0) ? 1 - jq1 : jq1)*((bj2 == 0) ? 1 - jq2 : jq2);
+               basis_j[jd_1] = 0.0;
+               basis_j[jd_2] = 0.0;
+
+               const int jj_loc = bj + 4*cj;
+
+               for (int ci=0; ci<dim; ++ci)
+               {
+                  const double iq1 = (ci == 0) ? qy : ((ci == 1) ? qz : qx);
+                  const double iq2 = (ci == 0) ? qz : ((ci == 1) ? qx : qy);
+                  const int id_0 = ci, id_1 = (ci + 1)%3, id_2 = (ci + 2)%3;
+
+                  for (int bi=0; bi<4; ++bi)
+                  {
+                     const int bi1 = bi%2, bi2 = bi/2;
+
+                     double curl_i[3];
+                     curl_i[id_0] = 0.0;
+                     curl_i[id_1] = ((bi1 == 0) ? iq1 - 1 : -iq1)*((bi2 == 0) ? 1 : -1);
+                     curl_i[id_2] = ((bi2 == 0) ? 1 - iq2 : iq2)*((bi1 == 0) ? 1 : -1);
+
+                     double basis_i[3];
+                     basis_i[id_0] = ((bi1 == 0) ? 1 - iq1 : iq1)*((bi2 == 0) ? 1 - iq2 : iq2);
+                     basis_i[id_1] = 0.0;
+                     basis_i[id_2] = 0.0;
+
+                     const int ii_loc = bi + 4*ci;
+
+                     // Only store the lower-triangular part of
+                     // the matrix (by symmetry).
+                     if (jj_loc > ii_loc) { continue; }
+
+                     double curl_curl = 0.0;
+                     curl_curl += Q6*curl_i[0]*curl_j[0];
+                     curl_curl += Q7*(curl_i[0]*curl_j[1] + curl_i[1]*curl_j[0]);
+                     curl_curl += Q8*(curl_i[0]*curl_j[2] + curl_i[2]*curl_j[0]);
+                     curl_curl += Q9*curl_i[1]*curl_j[1];
+                     curl_curl += Q10*(curl_i[1]*curl_j[2] + curl_i[2]*curl_j[1]);
+                     curl_curl += Q11*curl_i[2]*curl_j[2];
+
+                     double basis_basis = 0.0;
+                     basis_basis += Q0*basis_i[0]*basis_j[0];
+                     basis_basis += Q1*(basis_i[0]*basis_j[1] + basis_i[1]*basis_j[0]);
+                     basis_basis += Q2*(basis_i[0]*basis_j[2] + basis_i[2]*basis_j[0]);
+                     basis_basis += Q3*basis_i[1]*basis_j[1];
+                     basis_basis += Q4*(basis_i[1]*basis_j[2] + basis_i[2]*basis_j[1]);
+                     basis_basis += Q5*basis_i[2]*basis_j[2];
+
+                     const double val = dq*curl_curl + mq*basis_basis;
+                     AtomicAdd(local_mat(ci,bi, cj,bj), val);
+                  } // bi
+               } // ci
+            } // bj
+         } // cj
+      } // q
+      MFEM_SYNC_THREAD;
+
+      // Assemble the local matrix into the macro-element sparse matrix
+      // The nonzeros of the macro-element sparse matrix are ordered as
+      // follows:
+      //
+      // The axes are ordered relative to the direction of the basis
+      // vector, e.g. for x-vectors, the axes are (x,y,z), for
+      // y-vectors the axes are (y,z,x), and for z-vectors the axes are
+      // (z,x,y).
+      //
+      // The nonzeros are then given in "rotated lexicographic"
+      // ordering, according to these axes.
+      if (MFEM_THREAD_ID(x) == 0 && MFEM_THREAD_ID(y) == 0)
+      {
+         MFEM_FOREACH_THREAD(ii_loc,z,ne)
+         {
+            const int ci = ii_loc/4, bi = ii_loc%4;
+            const int i0 = 0, i1 = bi%2, i2 = bi/2;
+            const int id0 = ci, id1 = (ci+1)%3, id2 = (ci+2)%3;
+
+            int ii_lex[3];
+            ii_lex[id0] = i0, ii_lex[id1] = i1, ii_lex[id2] = i2;
+
+            const int nx = (ci == 0) ? o : op1, ny = (ci == 1) ? o : op1;
+            const int ii = ii_lex[0] + (ii_lex[1])*nx + (ii_lex[2])*nx*ny;
+
+            for (int jj_loc=0; jj_loc<ne; ++jj_loc)
+            {
+               const int cj = jj_loc/4, bj = jj_loc%4;
+               // add 3 to take modulus (rather than remainder) when
+               // (cj - ci) is negative
+               const int cj_rel = (3 + cj - ci)%3;
+               const int jd0 = cj_rel, jd1 = (cj_rel+1)%3, jd2 = (cj_rel+2)%3;
+
+               int jj_rel[3];
+               jj_rel[jd0] = 0, jj_rel[jd1] = bj%2, jj_rel[jd2] = bj/2;
+
+               const int d0 = jj_rel[0] - i0;
+               const int d1 = 1 + jj_rel[1] - i1;
+               const int d2 = 1 + jj_rel[2] - i2;
+
+               const int jj_off = (cj_rel == 0) ? d1 + 3*d2 :
+                                  (cj_rel == 1) ? 9 + d0 + 2*d1 + 4*d2 :
+                                  (cj_rel == 2) ? 21 + d0 + 2*d1 + 6*d2 : -1;
+
+               // Symmetry
+               const double val = (jj_loc <= ii_loc)
+                                  ? local_mat(ci,bi, cj,bj)
+                                  : local_mat(cj,bj, ci,bi);
+               AtomicAdd(V(jj_off, ii, ci, iel_ho), val);
+            }
+         }
+      }
+   });
+   SparseMapping3D<1>(sparse_mapping);
 }
 
 // Explicit template instantiations
@@ -572,7 +809,7 @@ template void BatchedLOR_ND::Assemble2D<6>();
 template void BatchedLOR_ND::Assemble2D<7>();
 template void BatchedLOR_ND::Assemble2D<8>();
 
-template void BatchedLOR_ND::Assemble3D<1>();
+//template void BatchedLOR_ND::Assemble3D<1>();  // explicitly specialized
 template void BatchedLOR_ND::Assemble3D<2>();
 template void BatchedLOR_ND::Assemble3D<3>();
 template void BatchedLOR_ND::Assemble3D<4>();

fem/lor/lor_rt.cpp

@@ -233,6 +233,87 @@ void BatchedLOR_RT::Assemble2D()
    }
 }
 
+template<int ORDER>
+static void SparseMapping3D(Array<int> &sparse_mapping)
+{
+   static constexpr int nnz_per_row = 11;
+   static constexpr int dim = 3;
+   static constexpr int o = ORDER;
+   static constexpr int op1 = ORDER + 1;
+   static constexpr int ndof_per_el = dim*o*o*op1;
+
+   sparse_mapping.SetSize(nnz_per_row*ndof_per_el);
+   sparse_mapping = -1;
+   auto map = Reshape(sparse_mapping.HostReadWrite(), nnz_per_row, ndof_per_el);
+
+   for (int ci=0; ci<dim; ++ci)
+   {
+      const int i_off = ci*o*o*op1;
+      const int id0 = ci;
+      const int id1 = (ci+1)%3;
+      const int id2 = (ci+2)%3;
+
+      const int nxi = (ci == 0) ? op1 : o;
+      const int nyi = (ci == 1) ? op1 : o;
+
+      for (int i0=0; i0<op1; ++i0)
+      {
+         for (int i1=0; i1<o; ++i1)
+         {
+            for (int i2=0; i2<o; ++i2)
+            {
+               int ii_lex[3];
+               ii_lex[id0] = i0;
+               ii_lex[id1] = i1;
+               ii_lex[id2] = i2;
+               const int ii_el = i_off + ii_lex[0] + ii_lex[1]*nxi + ii_lex[2]*nxi*nyi;
+
+               for (int cj_rel=0; cj_rel<dim; ++cj_rel)
+               {
+                  const int cj = (ci + cj_rel) % 3;
+                  const int j_off = cj*o*o*op1;
+
+                  const int nxj = (cj == 0) ? op1 : o;
+                  const int nyj = (cj == 1) ? op1 : o;
+
+                  const int j0_begin = (i0 > 0) ? i0-1 : i0;
+                  const int j0_end = (cj_rel == 0)
+                                     ? ((i0 < o) ? i0+1 : i0)
+                                     : ((i0 < o) ? i0 : i0-1);
+                  const int j1_begin = i1;
+                  const int j1_end = (cj_rel == 1) ? i1+1 : i1;
+                  const int j2_begin = i2;
+                  const int j2_end = (cj_rel == 2) ? i2+1 : i2;
+
+                  for (int j0=j0_begin; j0<=j0_end; ++j0)
+                  {
+                     const int d0 = 1 + j0 - i0;
+                     for (int j1=j1_begin; j1<=j1_end; ++j1)
+                     {
+                        const int d1 = j1 - i1;
+                        for (int j2=j2_begin; j2<=j2_end; ++j2)
+                        {
+                           const int d2 = j2 - i2;
+                           int jj_lex[3];
+                           jj_lex[id0] = j0;
+                           jj_lex[id1] = j1;
+                           jj_lex[id2] = j2;
+                           const int jj_el = j_off + jj_lex[0] + jj_lex[1]*nxj + jj_lex[2]*nxj*nyj;
+                           int jj_off;
+                           if (cj_rel == 0) { jj_off = d0; }
+                           else if (cj_rel == 1) { jj_off = 3 + d0 + 2*d1; }
+                           else /* if (cj_rel == 2) */ { jj_off = 7 + d0 + 2*d2; }
+                           map(jj_off, ii_el) = jj_el;
+                        }
+                     }
+                  }
+               }
+            }
+         }
+      }
+   }
+}
+
 template <int ORDER>
 void BatchedLOR_RT::Assemble3D()
 {
@@ -261,7 +342,7 @@ void BatchedLOR_RT::Assemble3D()
    sparse_ij.SetSize(nnz_per_row*ndof_per_el*nel_ho);
    auto V = Reshape(sparse_ij.Write(), nnz_per_row, o*o*op1, dim, nel_ho);
 
-   auto X = X_vert.Read();
+   const auto X = X_vert.Read();
 
    // Last thread dimension is lowered to avoid "too many resources" error
    MFEM_FORALL_3D(iel_ho, nel_ho, ORDER, ORDER, (ORDER>6)?4:ORDER,
@@ -474,76 +555,185 @@ void BatchedLOR_RT::Assemble3D()
          }
       }
    });
+   SparseMapping3D<ORDER>(sparse_mapping);
+}
 
-   sparse_mapping.SetSize(nnz_per_row*ndof_per_el);
-   sparse_mapping = -1;
-   auto map = Reshape(sparse_mapping.HostReadWrite(), nnz_per_row, ndof_per_el);
-   for (int ci=0; ci<dim; ++ci)
+template <>
+void BatchedLOR_RT::Assemble3D<1>()
+{
+
+   static constexpr int o = 1;
+   static constexpr int op1 = 2;
+   static constexpr int nf = 6; // number of faces in hexahedron
+   static constexpr int dim = 3;
+   static constexpr int ndof_per_el = dim*o*o*op1;
+   static constexpr int nnz_per_row = 11;
+   static constexpr int sz_local_mat = nf*nf;
+
+   const int nel_ho = fes_ho.GetNE();
+   const bool const_mq = c1.Size() == 1;
+   const auto MQ = const_mq
+                   ? Reshape(c1.Read(), 1, 1, 1, 1)
+                   : Reshape(c1.Read(), op1, op1, op1, nel_ho);
+   const bool const_dq = c2.Size() == 1;
+   const auto DQ = const_dq
+                   ? Reshape(c2.Read(), 1, 1, 1, 1)
+                   : Reshape(c2.Read(), op1, op1, op1, nel_ho);
+
+   sparse_ij.SetSize(nnz_per_row*ndof_per_el*nel_ho);
+   auto V = Reshape(sparse_ij.Write(), nnz_per_row, o*o*op1, dim, nel_ho);
+
+   const auto X = X_vert.Read();
+
+   MFEM_FORALL_3D(iel_ho, nel_ho, 8, 1, 1,
    {
-      const int i_off = ci*o*o*op1;
-      const int id0 = ci;
-      const int id1 = (ci+1)%3;
-      const int id2 = (ci+2)%3;
-
-      const int nxi = (ci == 0) ? op1 : o;
-      const int nyi = (ci == 1) ? op1 : o;
-
-      for (int i0=0; i0<op1; ++i0)
+      MFEM_FOREACH_THREAD(j,x,nnz_per_row)
       {
-         for (int i1=0; i1<o; ++i1)
+         for (int ix = 0; ix < op1; ++ix)
          {
-            for (int i2=0; i2<o; ++i2)
+            for (int c = 0; c < dim; ++c)
             {
-               int ii_lex[3];
-               ii_lex[id0] = i0;
-               ii_lex[id1] = i1;
-               ii_lex[id2] = i2;
-               const int ii_el = i_off + ii_lex[0] + ii_lex[1]*nxi + ii_lex[2]*nxi*nyi;
-
-               for (int cj_rel=0; cj_rel<dim; ++cj_rel)
-               {
-                  const int cj = (ci + cj_rel) % 3;
-                  const int j_off = cj*o*o*op1;
-
-                  const int nxj = (cj == 0) ? op1 : o;
-                  const int nyj = (cj == 1) ? op1 : o;
-
-                  const int j0_begin = (i0 > 0) ? i0-1 : i0;
-                  const int j0_end = (cj_rel == 0)
-                                     ? ((i0 < o) ? i0+1 : i0)
-                                     : ((i0 < o) ? i0 : i0-1);
-                  const int j1_begin = i1;
-                  const int j1_end = (cj_rel == 1) ? i1+1 : i1;
-                  const int j2_begin = i2;
-                  const int j2_end = (cj_rel == 2) ? i2+1 : i2;
-
-                  for (int j0=j0_begin; j0<=j0_end; ++j0)
-                  {
-                     const int d0 = 1 + j0 - i0;
-                     for (int j1=j1_begin; j1<=j1_end; ++j1)
-                     {
-                        const int d1 = j1 - i1;
-                        for (int j2=j2_begin; j2<=j2_end; ++j2)
-                        {
-                           const int d2 = j2 - i2;
-                           int jj_lex[3];
-                           jj_lex[id0] = j0;
-                           jj_lex[id1] = j1;
-                           jj_lex[id2] = j2;
-                           const int jj_el = j_off + jj_lex[0] + jj_lex[1]*nxj + jj_lex[2]*nxj*nyj;
-                           int jj_off;
-                           if (cj_rel == 0) { jj_off = d0; }
-                           else if (cj_rel == 1) { jj_off = 3 + d0 + 2*d1; }
-                           else /* if (cj_rel == 2) */ { jj_off = 7 + d0 + 2*d2; }
-                           map(jj_off, ii_el) = jj_el;
-                        }
-                     }
-                  }
-               }
+               V(j,ix,c,iel_ho) = 0.0;
             }
          }
       }
-   }
+      MFEM_SYNC_THREAD;
+
+      MFEM_FOREACH_THREAD(xyz,x,8)
+      {
+         const int qz = xyz%2, qy = (xyz/2)%2, qx = xyz/2/2;
+         static constexpr double w = 1.0/8.0;
+
+         double local_mat_[sz_local_mat];
+         DeviceTensor<4> local_mat(local_mat_, 3,2, 3,2);
+         for (int i=0; i<sz_local_mat; ++i) { local_mat[i] = 0.0; }
+
+         double vx[8], vy[8], vz[8];
+         LORVertexCoordinates3D<1>(X, iel_ho, 0,0,0, vx, vy, vz);
+
+         double J_[3*3];
+         DeviceTensor<2> J(J_, 3,3);
+         Jacobian3D(qx,qy,qz, vx,vy,vz, J);
+         const double detJ = Det3D(J);
+         const double w_detJ = w/detJ;
+
+         const double Q0 = w_detJ*(J(0,0)*J(0,0)+J(1,0)*J(1,0)+J(2,0)*J(2,0)); // 1,1
+         const double Q1 = w_detJ*(J(0,1)*J(0,0)+J(1,1)*J(1,0)+J(2,1)*J(2,0)); // 2,1
+         const double Q2 = w_detJ*(J(0,2)*J(0,0)+J(1,2)*J(1,0)+J(2,2)*J(2,0)); // 3,1
+         const double Q3 = w_detJ*(J(0,1)*J(0,1)+J(1,1)*J(1,1)+J(2,1)*J(2,1)); // 2,2
+         const double Q4 = w_detJ*(J(0,2)*J(0,1)+J(1,2)*J(1,1)+J(2,2)*J(2,1)); // 3,2
+         const double Q5 = w_detJ*(J(0,2)*J(0,2)+J(1,2)*J(1,2)+J(2,2)*J(2,2)); // 3,3
+         const double Q6 = w_detJ;
+
+         const double mq = const_mq ? MQ(0,0,0,0) : MQ(qx,qy,qz, iel_ho);
+         const double dq = const_dq ? DQ(0,0,0,0) : DQ(qx,qy,qz, iel_ho);
+
+         // Loop over x,y,z components. 0 => x, 1 => y, 2 => z
+         for (int cj=0; cj<dim; ++cj)//MFEM_FOREACH_THREAD(cj,y,dim)
+         {
+            const int jq0 = (cj == 0) ? qx : (cj == 1) ? qy : qz;
+            const int jd0 = cj, jd1 = (cj + 1)%3, jd2 = (cj + 2)%3;
+
+            for (int bj=0; bj<2; ++bj)//MFEM_FOREACH_THREAD(bj,z,2)  // 2 faces in each dim
+            {
+               const double div_j = (bj == 0) ? -1.0 : 1.0;
+
+               double basis_j[3];
+               basis_j[jd0] = (bj == jq0) ? 1.0 : 0.0;
+               basis_j[jd1] = 0.0;
+               basis_j[jd2] = 0.0;
+
+               const int jj_loc = bj + 2*cj;
+
+               for (int ci=0; ci<dim; ++ci)
+               {
+                  const double iq0 = (ci == 0) ? qx : ((ci == 1) ? qy : qz);
+                  const int id0 = ci, id1 = (ci + 1)%3, id2 = (ci + 2)%3;
+
+                  for (int bi=0; bi<2; ++bi)
+                  {
+                     const double div_i = (bi == 0) ? -1.0 : 1.0;
+                     double basis_i[3];
+                     basis_i[id0] = (bi == iq0) ? 1.0 : 0.0;
+                     basis_i[id1] = 0.0;
+                     basis_i[id2] = 0.0;
+
+                     const int ii_loc = bi + 2*ci;
+
+                     // Only store the lower-triangular part of
+                     // the matrix (by symmetry).
+                     if (jj_loc > ii_loc) { continue; }
+
+                     const double div_div = Q6*div_i*div_j;
+
+                     double basis_basis = 0.0;
+                     basis_basis += Q0*basis_i[0]*basis_j[0];
+                     basis_basis += Q1*(basis_i[0]*basis_j[1] + basis_i[1]*basis_j[0]);
+                     basis_basis += Q2*(basis_i[0]*basis_j[2] + basis_i[2]*basis_j[0]);
+                     basis_basis += Q3*basis_i[1]*basis_j[1];
+                     basis_basis += Q4*(basis_i[1]*basis_j[2] + basis_i[2]*basis_j[1]);
+                     basis_basis += Q5*basis_i[2]*basis_j[2];
+
+                     const double val = dq*div_div + mq*basis_basis;
+                     local_mat(ci,bi, cj,bj) += val;
+                  } // bi
+               } // ci
+            } // bj
+         } // cj
+
+         // Assemble the local matrix into the macro-element sparse matrix
+         // The nonzeros of the macro-element sparse matrix are ordered as
+         // follows:
+         //
+         // The axes are ordered relative to the direction of the basis
+         // vector, e.g. for x-vectors, the axes are (x,y,z), for
+         // y-vectors the axes are (y,z,x), and for z-vectors the axes are
+         // (z,x,y).
+         //
+         // The nonzeros are then given in "rotated lexicographic"
+         // ordering, according to these axes.
+         for (int ii_loc=0; ii_loc<nf; ++ii_loc)
+         {
+            const int ci = ii_loc/2, bi = ii_loc%2;
+            const int id0 = ci, id1 = (ci+1)%3, id2 = (ci+2)%3;
+            const int i0 = bi, i1 = 0, i2 = 0;
+
+            int ii_lex[3];
+            ii_lex[id0] = i0, ii_lex[id1] = i1, ii_lex[id2] = i2;
+
+            const int nx = (ci == 0) ? op1 : o;
+            const int ny = (ci == 1) ? op1 : o;
+            const int ii = ii_lex[0] + ii_lex[1]*nx + ii_lex[2]*nx*ny;
+
+            for (int jj_loc=0; jj_loc<nf; ++jj_loc)
+            {
+               const int cj = jj_loc/2, bj = jj_loc%2;
+               // add 3 to take modulus (rather than remainder) when
+               // (cj - ci) is negative
+               const int cj_rel = (3 + cj - ci)%3;
+               const int jd0 = cj_rel, jd1 = (cj_rel+1)%3, jd2 = (cj_rel+2)%3;
+
+               int jj_rel[3];
+               jj_rel[jd0] = bj, jj_rel[jd1] = 0, jj_rel[jd2] = 0;
+
+               const int d0 = jj_rel[0] - i0 + 1;
+               const int d1 = jj_rel[1] - i1;
+               const int d2 = jj_rel[2] - i2;
+
+               const int jj_off = (cj_rel == 0) ?d0 :
+                                  (cj_rel == 1) ?  3 + d0 + 2*d1 :
+                                  (cj_rel == 2) ?  7 + d0 + 2*d2 : -1;
+
+               // Symmetry
+               const double val = (jj_loc <= ii_loc)
+                                  ? local_mat(ci,bi, cj,bj)
+                                  : local_mat(cj,bj, ci,bi);
+               AtomicAdd(V(jj_off, ii, ci, iel_ho), val);
+            } // jj_loc
+         } // ii_loc
+      } // q
+   });
+   SparseMapping3D<1>(sparse_mapping);
 }
 
 // Explicit template instantiations
@@ -556,7 +746,7 @@ template void BatchedLOR_RT::Assemble2D<6>();
 template void BatchedLOR_RT::Assemble2D<7>();
 template void BatchedLOR_RT::Assemble2D<8>();
 
-template void BatchedLOR_RT::Assemble3D<1>();
+//template void BatchedLOR_RT::Assemble3D<1>(); // explicitly specialized
 template void BatchedLOR_RT::Assemble3D<2>();
 template void BatchedLOR_RT::Assemble3D<3>();
 template void BatchedLOR_RT::Assemble3D<4>();

fem/pbilinearform.cpp

@@ -16,6 +16,9 @@
 #include "fem.hpp"
 #include "../general/sort_pairs.hpp"
 
+#define MFEM_NVTX_COLOR Crimson
+#include "../general/nvtx.hpp"
+
 namespace mfem
 {
 
@@ -124,6 +127,7 @@ void ParBilinearForm::pAllocMat()
 void ParBilinearForm::ParallelRAP(SparseMatrix &loc_A, OperatorHandle &A,
                                   bool steal_loc_A)
 {
+   NVTX("RAP");
    ParFiniteElementSpace &pfespace = *ParFESpace();
 
    // Create a block diagonal parallel matrix

fem/qinterp/eval_by_nodes.cpp

@@ -128,8 +128,8 @@ void TensorValues<QVectorLayout::byNODES>(const int NE,
 
          default:
          {
-            constexpr int MD = 8;
-            constexpr int MQ = 8;
+            constexpr int MD = 9;
+            constexpr int MQ = 9;
             MFEM_VERIFY(D1D <= MD, "Orders higher than " << MD-1
                         << " are not supported!");
             MFEM_VERIFY(Q1D <= MQ, "Quadrature rules with more than "

fem/qinterp/eval_by_vdim.cpp

@@ -102,8 +102,8 @@ void TensorValues<QVectorLayout::byVDIM>(const int NE,
 
          default:
          {
-            constexpr int MD = 8;
-            constexpr int MQ = 8;
+            constexpr int MD = 9;
+            constexpr int MQ = 9;
             MFEM_VERIFY(D1D <= MD, "Orders higher than " << MD-1
                         << " are not supported!");
             MFEM_VERIFY(Q1D <= MQ, "Quadrature rules with more than "

fem/qinterp/grad_by_nodes.cpp

@@ -110,8 +110,8 @@ void TensorDerivatives<QVectorLayout::byNODES>(const int NE,
          case 0x348: return Derivatives3D<L,P,3,4,8>(NE,B,G,J,X,Y);
          default:
          {
-            constexpr int MD = 8;
-            constexpr int MQ = 8;
+            constexpr int MD = 9;
+            constexpr int MQ = 9;
             MFEM_VERIFY(D1D <= MD, "Orders higher than " << MD-1
                         << " are not supported!");
             MFEM_VERIFY(Q1D <= MQ, "Quadrature rules with more than "

general/debug.hpp

@@ -0,0 +1,178 @@
+// 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_DEBUG_HPP
+#define MFEM_DEBUG_HPP
+
+#include <string>
+#include <cstring>
+#include <iomanip>
+#include <iostream>
+
+#include "globals.hpp"
+
+#ifdef MFEM_USE_MPI
+#include <mpi.h>
+#endif
+
+namespace mfem
+{
+
+class Debug
+{
+   const bool debug = false;
+public:
+   inline Debug() {}
+
+   inline Debug(const int mpi_rank,
+                const char *FILE, const int LINE,
+                const char *FUNC, int COLOR): debug(true)
+   {
+      if (!debug) { return; }
+      const char *base = Strrnchr(FILE,'/', 2);
+      const char *file = base ? base + 1 : FILE;
+      const uint8_t color = COLOR ? COLOR : 20 + Checksum8(FILE) % 210;
+      mfem::out << "\033[38;5;" << std::to_string(color) << "m";
+      mfem::out << mpi_rank << std::setw(30) << file << ":";
+      mfem::out << "\033[2m" << std::setw(4) << LINE << "\033[22m: ";
+      if (FUNC) { mfem::out << "[" << FUNC << "] "; }
+      mfem::out << "\033[1m";
+   }
+
+   ~Debug()
+   {
+      if (!debug) { return; }
+      mfem::out << "\033[m";
+      mfem::out << std::endl;
+   }
+
+   template <typename T>
+   inline void operator<<(const T &arg) const noexcept { mfem::out << arg; }
+
+   template<typename T, typename... Args>
+   inline void operator()(const char *fmt, const T &arg,
+                          Args... args) const noexcept
+   {
+      if (!debug) { return; }
+      for (; *fmt != '\0'; fmt++ )
+      {
+         if (*fmt == '%')
+         {
+            fmt++;
+            const char c = *fmt;
+            if (c == 'p') { operator<<(arg); }
+            if (c == 's' || c == 'd' || c == 'f') { operator<<(arg); }
+            if (c == 'x' || c == 'X')
+            {
+               mfem::out << std::hex;
+               if (c == 'X') { mfem::out << std::uppercase; }
+               operator<<(arg);
+               mfem::out << std::nouppercase << std::dec;
+            }
+            if (c == '.')
+            {
+               fmt++;
+               const char c2 = *fmt;
+               char num[8] = { 0 };
+               for (int k = 0; *fmt != '\0'; fmt++, k++)
+               {
+                  if (*fmt == 'e' || *fmt == 'f') { break; }
+                  if (*fmt < 0x30 || *fmt > 0x39) { break; }
+                  num[k] = *fmt;
+               }
+               const int fx = std::atoi(num);
+               if (c2 == 'e') { mfem::out << std::scientific; }
+               if (c2 == 'f') { mfem::out << std::fixed; }
+               mfem::out << std::setprecision(fx);
+               operator<<(arg);
+               mfem::out << std::setprecision(6);
+            }
+            return operator()(fmt + 1, args...);
+         }
+         operator<<(*fmt);
+      }
+   }
+
+   template<typename T>
+   inline void operator()(const T &arg) const noexcept
+   {
+      if (!debug) { return; }
+      operator<<(arg);
+   }
+
+   inline void operator()() const noexcept { }
+
+public:
+   static const Debug Set(const char *FILE, const int LINE, const char *FUNC,
+                          int COLOR = 0)
+   {
+      static int mpi_dbg = 0, mpi_rank = 0;
+      static bool env_mpi = false, env_dbg = false;
+      static bool ini_dbg = false;
+      if (!ini_dbg)
+      {
+         const char *DBG = getenv("MFEM_DEBUG");
+         const char *MPI = getenv("MFEM_DEBUG_MPI");
+         env_dbg = DBG != nullptr;
+         env_mpi = MPI != nullptr;
+#ifdef MFEM_USE_MPI
+         int mpi_ini = false;
+         MPI_Initialized(&mpi_ini);
+         if (mpi_ini) { MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank); }
+         mpi_dbg = atoi(env_mpi ? MPI : "0");
+#endif
+         ini_dbg = true;
+      }
+      const bool debug = (env_dbg && (!env_mpi || mpi_rank == mpi_dbg));
+      return debug ? Debug(mpi_rank, FILE, LINE, FUNC, COLOR) : Debug();
+   }
+
+private:
+   inline uint8_t Checksum8(const char *bfr)
+   {
+      unsigned int chk = 0;
+      size_t len = strlen(bfr);
+      for (; len; len--,bfr++) { chk += static_cast<unsigned int>(*bfr); }
+      return (uint8_t) chk;
+   }
+
+   inline const char *Strrnchr(const char *s, const unsigned char c, int n)
+   {
+      size_t len = strlen(s);
+      char *p = const_cast<char*>(s) + len - 1;
+      for (; n; n--,p--,len--)
+      {
+         for (; len; p--,len--)
+            if (*p == c) { break; }
+         if (!len) { return nullptr; }
+         if (n == 1) { return p; }
+      }
+      return nullptr;
+   }
+
+};
+
+#ifndef MFEM_DEBUG_COLOR
+#define MFEM_DEBUG_COLOR 0
+#endif
+
+#define dbg(...) \
+    mfem::Debug::Set(__FILE__,__LINE__,__FUNCTION__,MFEM_DEBUG_COLOR).\
+    operator()(__VA_ARGS__)
+
+} // mfem namespace
+
+#define DBG(...) { printf("\033[33m");  \
+                   printf(__VA_ARGS__); \
+                   printf(" \n\033[m"); \
+                   fflush(0); }
+
+#endif // MFEM_DEBUG_HPP

general/nvtx.hpp

@@ -0,0 +1,312 @@
+// Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights
+// reserved. See files LICENSE and NOTICE for details.
+//
+// This file is part of CEED, a collection of benchmarks, miniapps, software
+// libraries and APIs for efficient high-order finite element and spectral
+// element discretizations for exascale applications. For more information and
+// source code availability see http://github.com/ceed.
+//
+// The CEED research is supported by the Exascale Computing Project 17-SC-20-SC,
+// a collaborative effort of two U.S. Department of Energy organizations (Office
+// of Science and the National Nuclear Security Administration) responsible for
+// the planning and preparation of a capable exascale ecosystem, including
+// software, applications, hardware, advanced system engineering and early
+// testbed platforms, in support of the nation's exascale computing imperative.
+
+#ifndef MFEM_NVVP_HPP
+#define MFEM_NVVP_HPP
+
+#include <string>
+#include <cassert>
+
+#include "globals.hpp"
+#include "backends.hpp"
+
+#if defined(MFEM_USE_CUDA) //||1
+#include <cuda.h>
+#include <nvToolsExt.h>
+#include <cudaProfiler.h>
+#include <cuda_runtime.h>
+#else // MFEM_USE_CUDA
+typedef struct
+{
+   uint16_t version;
+   uint16_t size;
+   uint32_t category;
+   int32_t colorType;
+   uint32_t color;
+   int32_t payloadType;
+   int32_t reserved0;
+   int32_t messageType;
+   struct { const char* ascii; } message;
+} nvtxEventAttributes_t;
+#define NVTX_VERSION 2
+#define NVTX_COLOR_ARGB 1
+#define NVTX_MESSAGE_TYPE_ASCII 1
+#define NVTX_EVENT_ATTRIB_STRUCT_SIZE \
+    ((uint16_t)(sizeof(nvtxEventAttributes_t)))
+inline int nvtxRangePushEx(const nvtxEventAttributes_t*) { return 0; }
+inline int nvtxRangePop(void) { return 0; }
+#endif // MFEM_USE_CUDA
+
+namespace mfem
+{
+
+// https://en.wikipedia.org/wiki/Web_colors#Extended_colors
+typedef enum
+{
+   Black, NavyBlue, DarkBlue, MediumBlue, Blue, DarkGreen, WebGreen, Teal,
+   DarkCyan, DeepSkyBlue, DarkTurquoise, MediumSpringGreen, Green, Lime,
+   SpringGreen, Aqua, Cyan, MidnightBlue, DodgerBlue, LightSeaGreen,
+   ForestGreen, SeaGreen, DarkSlateGray, LimeGreen, MediumSeaGreen,
+   Turquoise, RoyalBlue, SteelBlue, DarkSlateBlue, MediumTurquoise, Indigo,
+   DarkOliveGreen, CadetBlue, Cornflower, RebeccaPurple, MediumAquamarine,
+   DimGray, SlateBlue, OliveDrab, SlateGray, LightSlateGray,
+   MediumSlateBlue, LawnGreen, WebMaroon, WebPurple, Chartreuse,
+   Aquamarine, Olive, WebGray, SkyBlue, LightSkyBlue, BlueViolet, DarkRed,
+   DarkMagenta, SaddleBrown, DarkSeaGreen, LightGreen, MediumPurple,
+   DarkViolet, PaleGreen, DarkOrchid, YellowGreen, Purple, Sienna, Brown,
+   DarkGray, LightBlue, GreenYellow, PaleTurquoise, Maroon,
+   LightSteelBlue, PowderBlue, Firebrick, DarkGoldenrod, MediumOrchid,
+   RosyBrown, DarkKhaki, Gray, Silver, MediumVioletRed, IndianRed, Peru,
+   Chocolate, Tan, LightGray, Thistle, Orchid, Goldenrod, PaleVioletRed,
+   Crimson, Gainsboro, Plum, Burlywood, LightCyan, Lavender, DarkSalmon,
+   Violet, PaleGoldenrod, LightCoral, Khaki, AliceBlue, Honeydew, Azure,
+   SandyBrown, Wheat, Beige, WhiteSmoke, MintCream, GhostWhite, Salmon,
+   AntiqueWhite, Linen, LightGoldenrod, OldLace, Red, Fuchsia, Magenta,
+   DeepPink, OrangeRed, Tomato, HotPink, Coral, DarkOrange, LightSalmon,
+   Orange, LightPink, Pink, Gold, PeachPuff, NavajoWhite, Moccasin,
+   Bisque, MistyRose, BlanchedAlmond, PapayaWhip, LavenderBlush, Seashell,
+   Cornsilk, LemonChiffon, FloralWhite, Snow, Yellow, LightYellow, Ivory,
+   White
+} COLOR_NAMES;
+
+static constexpr uint32_t HEX_COLORS[] =
+{
+   0x000000, 0x000080, 0x00008B, 0x0000CD, 0x0000FF, 0x006400, 0x008000,
+   0x008080, 0x008B8B, 0x00BFFF, 0x00CED1, 0x00FA9A, 0x00FF00, 0x00FF00,
+   0x00FF7F, 0x00FFFF, 0x00FFFF, 0x191970, 0x1E90FF, 0x20B2AA, 0x228B22,
+   0x2E8B57, 0x2F4F4F, 0x32CD32, 0x3CB371, 0x40E0D0, 0x4169E1, 0x4682B4,
+   0x483D8B, 0x48D1CC, 0x4B0082, 0x556B2F, 0x5F9EA0, 0x6495ED, 0x663399,
+   0x66CDAA, 0x696969, 0x6A5ACD, 0x6B8E23, 0x708090, 0x778899, 0x7B68EE,
+   0x7CFC00, 0x7F0000, 0x7F007F, 0x7FFF00, 0x7FFFD4, 0x808000, 0x808080,
+   0x87CEEB, 0x87CEFA, 0x8A2BE2, 0x8B0000, 0x8B008B, 0x8B4513, 0x8FBC8F,
+   0x90EE90, 0x9370DB, 0x9400D3, 0x98FB98, 0x9932CC, 0x9ACD32, 0xA020F0,
+   0xA0522D, 0xA52A2A, 0xA9A9A9, 0xADD8E6, 0xADFF2F, 0xAFEEEE, 0xB03060,
+   0xB0C4DE, 0xB0E0E6, 0xB22222, 0xB8860B, 0xBA55D3, 0xBC8F8F, 0xBDB76B,
+   0xBEBEBE, 0xC0C0C0, 0xC71585, 0xCD5C5C, 0xCD853F, 0xD2691E, 0xD2B48C,
+   0xD3D3D3, 0xD8BFD8, 0xDA70D6, 0xDAA520, 0xDB7093, 0xDC143C, 0xDCDCDC,
+   0xDDA0DD, 0xDEB887, 0xE0FFFF, 0xE6E6FA, 0xE9967A, 0xEE82EE, 0xEEE8AA,
+   0xF08080, 0xF0E68C, 0xF0F8FF, 0xF0FFF0, 0xF0FFFF, 0xF4A460, 0xF5DEB3,
+   0xF5F5DC, 0xF5F5F5, 0xF5FFFA, 0xF8F8FF, 0xFA8072, 0xFAEBD7, 0xFAF0E6,
+   0xFAFAD2, 0xFDF5E6, 0xFF0000, 0xFF00FF, 0xFF00FF, 0xFF1493, 0xFF4500,
+   0xFF6347, 0xFF69B4, 0xFF7F50, 0xFF8C00, 0xFFA07A, 0xFFA500, 0xFFB6C1,
+   0xFFC0CB, 0xFFD700, 0xFFDAB9, 0xFFDEAD, 0xFFE4B5, 0xFFE4C4, 0xFFE4E1,
+   0xFFEBCD, 0xFFEFD5, 0xFFF0F5, 0xFFF5EE, 0xFFF8DC, 0xFFFACD, 0xFFFAF0,
+   0xFFFAFA, 0xFFFF00, 0xFFFFE0, 0xFFFFF0, 0xFFFFFF
+};
+
+static const int NUM_HEX_COLORS = sizeof(HEX_COLORS)/sizeof(uint32_t);
+
+class Nvtx
+{
+   const bool nvtx = false;
+   const bool enforce_kernel_sync = false;
+   const char *base, *file;
+   const uint32_t color = Black;
+   mutable std::string ascii;
+   mutable nvtxEventAttributes_t event;
+
+public:
+   Nvtx() { }
+
+   Nvtx(bool enforce_kernel_sync,
+        const char *FILE, const int LINE, const char *FUNC, uint32_t COLOR):
+      nvtx(true),
+      enforce_kernel_sync(enforce_kernel_sync),
+      base(Strrnchr(FILE,'/', 2)),
+      file(base ? base + 1 : FILE),
+      color(COLOR),
+      ascii(file),
+      event({})
+   {
+      event.version = NVTX_VERSION;
+      event.size = NVTX_EVENT_ATTRIB_STRUCT_SIZE;
+      event.colorType = NVTX_COLOR_ARGB;
+      event.color = HEX_COLORS[color % NUM_HEX_COLORS];
+      event.messageType = NVTX_MESSAGE_TYPE_ASCII;
+
+      ascii += ":";
+      ascii += std::to_string(LINE);
+      ascii += ":[";
+      ascii += FUNC;
+      ascii += "] ";
+   }
+
+   Nvtx(const char *title, int color = Wheat, bool enforce_kernel_sync = true):
+      nvtx(true),
+      enforce_kernel_sync(enforce_kernel_sync),
+      color(color),
+      ascii(title),
+      event({})
+   {
+      event.version = NVTX_VERSION;
+      event.size = NVTX_EVENT_ATTRIB_STRUCT_SIZE;
+      event.colorType = NVTX_COLOR_ARGB;
+      event.color = HEX_COLORS[color % NUM_HEX_COLORS];
+      event.messageType = NVTX_MESSAGE_TYPE_ASCII;
+      event.message.ascii = ascii.c_str();
+      nvtxRangePushEx(&event); // push
+   }
+
+   ~Nvtx()
+   {
+      if (!nvtx) { return; }
+      if (enforce_kernel_sync)
+      {
+         nvtxEventAttributes_t eks = {};
+         eks.version = NVTX_VERSION;
+         eks.size = NVTX_EVENT_ATTRIB_STRUCT_SIZE;
+         eks.category = 0; // user value
+         eks.colorType = NVTX_COLOR_ARGB;
+         eks.messageType = NVTX_MESSAGE_TYPE_ASCII;
+         eks.message.ascii = "Sync";
+         eks.color = HEX_COLORS[Yellow];
+         nvtxRangePushEx(&eks);
+         MFEM_STREAM_SYNC;
+         nvtxRangePop();
+      }
+      nvtxRangePop(); // pop
+   }
+
+   // used through MFEM_NVTX
+   inline void operator()() const noexcept
+   {
+      event.message.ascii = ascii.c_str();
+      nvtxRangePushEx(&event); // push
+   }
+
+   template<typename T>
+   inline void operator()(const T &arg) const noexcept
+   {
+      if (!nvtx) { return; }
+      operator<<(arg);
+      event.message.ascii = ascii.c_str();
+      nvtxRangePushEx(&event); // push
+   }
+
+   template<typename T, typename... Args>
+   inline void operator()(const char *fmt, const T &arg,
+                          Args... args) const noexcept
+   {
+      if (!nvtx) { return; }
+
+      for (; *fmt != '\0'; fmt++ )
+      {
+         if (*fmt == '%')
+         {
+            fmt++;
+            const char c = *fmt;
+            if (c == 'p') { operator<<(arg); }
+            if (c == 's' || c == 'd' || c == 'f') { operator<<(arg); }
+            if (c == 'x' || c == 'X')
+            {
+               mfem::out << std::hex;
+               if (c == 'X') { mfem::out << std::uppercase; }
+               operator<<(arg);
+               mfem::out << std::nouppercase << std::dec;
+            }
+            if (c == '.')
+            {
+               fmt++;
+               const char c2 = *fmt;
+               char num[8] = { 0 };
+               for (int k = 0; *fmt != '\0'; fmt++, k++)
+               {
+                  if (*fmt == 'e' || *fmt == 'f') { break; }
+                  if (*fmt < 0x30 || *fmt > 0x39) { break; }
+                  num[k] = *fmt;
+               }
+               const int fx = std::atoi(num);
+               if (c2 == 'e') { mfem::out << std::scientific; }
+               if (c2 == 'f') { mfem::out << std::fixed; }
+               mfem::out << std::setprecision(fx);
+               operator<<(arg);
+               mfem::out << std::setprecision(6);
+            }
+            return operator()(fmt + 1, args...);
+         }
+         operator<<(*fmt);
+      }
+      // should never be here
+      assert(false);
+   }
+
+   template <typename T>
+   inline void operator<<(const T &arg) const noexcept
+   {
+      if (!nvtx) { return; }
+      ascii += arg;
+   }
+
+   inline void operator<<(const int &arg) const noexcept
+   {
+      if (!nvtx) { return; }
+      ascii += std::to_string(arg);
+   }
+
+
+public:
+   static const Nvtx Set(const char *FILE, const int LINE, const char *FUNC,
+                         uint32_t COLOR)
+   {
+      static bool env_nvtx = false;
+      static bool env_eks = false;
+      static bool ini_nvtx = false;
+      if (!ini_nvtx)
+      {
+         env_nvtx = getenv("MFEM_NVTX") != nullptr;
+         env_eks = getenv("MFEM_EKS") != nullptr;
+         ini_nvtx = true;
+      }
+      return env_nvtx ? Nvtx(env_eks, FILE, LINE, FUNC, COLOR) : Nvtx();
+   }
+
+private:
+   inline const char *Strrnchr(const char *s, const unsigned char c, int n)
+   {
+      size_t len = strlen(s);
+      char *p = const_cast<char*>(s) + len - 1;
+      for (; n; n--,p--,len--)
+      {
+         for (; len; p--,len--)
+            if (*p == c) { break; }
+         if (!len) { return nullptr; }
+         if (n == 1) { return p; }
+      }
+      return nullptr;
+   }
+};
+
+#ifndef MFEM_NVTX_COLOR
+#define MFEM_NVTX_COLOR SeaGreen
+#endif
+
+// Helpers for generating unique variable names
+#define NVTX_PRIVATE_NAME(name) NVTX_PRIVATE_CONCAT(name, __LINE__)
+#define NVTX_PRIVATE_CONCAT(a, b) NVTX_PRIVATE_CONCAT2(a, b)
+#define NVTX_PRIVATE_CONCAT2(a, b) a##b
+
+// temporary object which is only alive in the expression
+// __PRETTY_FUNCTION__
+#define NVTX(...) \
+    mfem::Nvtx NVTX_PRIVATE_NAME(nvtx) = \
+    mfem::Nvtx::Set(__FILE__,__LINE__,__FUNCTION__,MFEM_NVTX_COLOR);\
+    NVTX_PRIVATE_NAME(nvtx).operator()(__VA_ARGS__)
+
+#define MFEM_NVTX NVTX()
+
+} // namespace mfem
+
+#endif // MFEM_NVVP_HPP

linalg/hypre.cpp

@@ -22,6 +22,8 @@
 #include <cmath>
 #include <cstdlib>
 
+#include "../general/nvtx.hpp"
+
 using namespace std;
 
 namespace mfem
@@ -1732,6 +1734,10 @@ void HypreParMatrix::ResetTranspose() const
 HYPRE_Int HypreParMatrix::Mult(HypreParVector &x, HypreParVector &y,
                                double a, double b) const
 {
+#undef MFEM_NVTX_COLOR
+#define MFEM_NVTX_COLOR Crimson
+   NVTX("SpMV");
+
    x.HypreRead();
    (b == 0.0) ? y.HypreWrite() : y.HypreReadWrite();
    return hypre_ParCSRMatrixMatvec(a, A, x, b, y);
@@ -1739,6 +1745,10 @@ HYPRE_Int HypreParMatrix::Mult(HypreParVector &x, HypreParVector &y,
 
 void HypreParMatrix::Mult(double a, const Vector &x, double b, Vector &y) const
 {
+#undef MFEM_NVTX_COLOR
+#define MFEM_NVTX_COLOR Crimson
+   NVTX("SpMV");
+
    MFEM_ASSERT(x.Size() == Width(), "invalid x.Size() = " << x.Size()
                << ", expected size = " << Width());
    MFEM_ASSERT(y.Size() == Height(), "invalid y.Size() = " << y.Size()
@@ -1797,6 +1807,10 @@ void HypreParMatrix::Mult(double a, const Vector &x, double b, Vector &y) const
 void HypreParMatrix::MultTranspose(double a, const Vector &x,
                                    double b, Vector &y) const
 {
+#undef MFEM_NVTX_COLOR
+#define MFEM_NVTX_COLOR Crimson
+   NVTX("SpMV_T");
+
    MFEM_ASSERT(x.Size() == Height(), "invalid x.Size() = " << x.Size()
                << ", expected size = " << Height());
    MFEM_ASSERT(y.Size() == Width(), "invalid y.Size() = " << y.Size()
@@ -3880,6 +3894,10 @@ void HypreSolver::Setup(const HypreParVector &b, HypreParVector &x) const
 {
    if (setup_called) { return; }
 
+#undef MFEM_NVTX_COLOR
+#define MFEM_NVTX_COLOR Indigo
+   NVTX("AMG Setup");
+
    MFEM_VERIFY(A != NULL, "HypreParMatrix A is missing");
 
    HYPRE_Int err_flag = SetupFcn()(*this, *A, b, x);
@@ -3905,6 +3923,10 @@ void HypreSolver::Setup(const Vector &b, Vector &x) const
 
 void HypreSolver::Mult(const HypreParVector &b, HypreParVector &x) const
 {
+#undef MFEM_NVTX_COLOR
+#define MFEM_NVTX_COLOR MediumSpringGreen
+   NVTX("AMG V-cycle");
+
    HYPRE_Int err_flag;
    if (A == NULL)
    {

linalg/solvers.cpp

@@ -20,6 +20,9 @@
 #include <cmath>
 #include <set>
 
+#define MFEM_NVTX_COLOR DarkGreen
+#include "../general/nvtx.hpp"
+
 namespace mfem
 {
 
@@ -710,6 +713,8 @@ void CGSolver::UpdateVectors()
 
 void CGSolver::Mult(const Vector &b, Vector &x) const
 {
+   NVTX("CG");
+
    int i;
    double r0, den, nom, nom0, betanom, alpha, beta;
 

miniapps/solvers/plor_solvers.cpp

@@ -84,7 +84,7 @@ int main(int argc, char *argv[])
    int order = 3;
    const char *fe = "h";
    const char *device_config = "cpu";
-   bool visualization = true;
+   bool visualization = false;
 
    OptionsParser args(argc, argv);
    args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
@@ -123,6 +123,9 @@ int main(int argc, char *argv[])
    for (int l = 0; l < par_ref_levels; l++) { mesh.UniformRefinement(); }
    serial_mesh.Clear();
 
+   HYPRE_Int global_ne = mesh.GetGlobalNE();
+   if (Mpi::Root()) { cout << "Number of elements: " << global_ne << endl; }
+
    if (mesh.ncmesh && (RT || ND))
    { MFEM_ABORT("LOR AMS and ADS solvers are not supported with AMR meshes."); }
 
@@ -188,15 +191,23 @@ int main(int argc, char *argv[])
    unique_ptr<Solver> solv_lor;
    if (H1 || L2)
    {
-      solv_lor.reset(new LORSolver<HypreBoomerAMG>(a, ess_dofs));
+      auto *lor_solver = new LORSolver<HypreBoomerAMG>(a, ess_dofs);
+      HYPRE_BigInt nnz =
+         lor_solver->GetLOR().GetAssembledSystem().As<HypreParMatrix>()->NNZ();
+      if (Mpi::Root()) { cout << "Number of NNZ:  " << nnz << endl; }
+
+      solv_lor.reset(lor_solver);
+      ((LORSolver<HypreBoomerAMG>&)*solv_lor).GetSolver().Setup(B, X);
    }
    else if (RT && dim == 3)
    {
       solv_lor.reset(new LORSolver<HypreADS>(a, ess_dofs));
+      ((LORSolver<HypreADS>&)*solv_lor).GetSolver().Setup(B, X);
    }
    else
    {
       solv_lor.reset(new LORSolver<HypreAMS>(a, ess_dofs));
+      ((LORSolver<HypreAMS>&)*solv_lor).GetSolver().Setup(B, X);
    }
 
    CGSolver cg(MPI_COMM_WORLD);
@@ -208,14 +219,16 @@ int main(int argc, char *argv[])
    cg.SetPreconditioner(*solv_lor);
    cg.Mult(B, X);
 
-   a.RecoverFEMSolution(X, b, x);
+   return 0;
 
-   double er =
-      (H1 || L2) ? x.ComputeL2Error(u_coeff) : x.ComputeL2Error(u_vec_coeff);
-   if (Mpi::Root()) { cout << "L2 error: " << er << endl; }
+   a.RecoverFEMSolution(X, b, x);
 
    if (visualization)
    {
+      double er =
+         (H1 || L2) ? x.ComputeL2Error(u_coeff) : x.ComputeL2Error(u_vec_coeff);
+      if (Mpi::Root()) { cout << "L2 error: " << er << endl; }
+
       // Save the solution and mesh to disk. The output can be viewed using
       // GLVis as follows: "glvis -np <np> -m mesh -g sol"
       x.Save("sol");

miniapps/solvers/run_amr.sh

@@ -0,0 +1,15 @@
+#!/bin/bash
+
+ORDERS="1 2 3 4 5 6 7"
+
+for p in $ORDERS
+do
+   MFEM_EKS=1 MFEM_NVTX=1 lrun -n 12 nsys profile -t nvtx,cuda -s none -o amr/prof/o${p}_%q{OMPI_COMM_WORLD_RANK} -f true ./plor_solvers -m ./amr_mesh.mesh -o $p -rs 0 -rp 0 -d cuda | tee amr/out/o$p.txt
+done
+
+# After creating the NVTX files, convert to CSV using the following:
+
+# for p in 1 2 3 4 5 6 7
+# do
+#    nsys stats -r nvtxsum -f csv o${p}_0.nsys-rep > o${p}.csv
+# done

tests/benchmarks/bench_lor.cpp

@@ -0,0 +1,428 @@
+// 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 <cmath>
+#include <cassert>
+
+#include "bench.hpp"
+
+#include "fem/lor/lor_ads.hpp"
+#include "fem/lor/lor_ams.hpp"
+
+#ifdef MFEM_USE_BENCHMARK
+
+#include "fem/lor/lor.hpp"
+
+#define MFEM_DEBUG_COLOR 119
+#include "general/debug.hpp"
+
+#define MFEM_NVTX_COLOR Lime
+#include "general/nvtx.hpp"
+
+Mesh MakeCartesianMesh(int p, int requested_ndof, int dim)
+{
+   const int ne = std::max(1, (int)std::ceil(requested_ndof / pow(p, dim)));
+   if (dim == 2)
+   {
+      const int nx = sqrt(ne);
+      const int ny = ne / nx;
+      return Mesh::MakeCartesian2D(nx, ny, Element::QUADRILATERAL);
+   }
+   else
+   {
+      const int nx = cbrt(ne);
+      const int ny = sqrt(ne / nx);
+      const int nz = ne / nx / ny;
+      return Mesh::MakeCartesian3D(nx, ny, nz, Element::HEXAHEDRON);
+   }
+}
+
+ParMesh MakeParCartesianMesh(int p, int requested_ndof, int dim)
+{
+   Mesh mesh = MakeCartesianMesh(p, requested_ndof, dim);
+   return ParMesh(MPI_COMM_WORLD, mesh);
+}
+
+struct RT_LORBench
+{
+   ParMesh mesh;
+   RT_FECollection fec_ho;
+   ParFiniteElementSpace fes_ho;
+
+   BatchedLORAssembly lor;
+   BatchedLOR_ADS ads;
+   OperatorHandle A_lor;
+
+   LORSolver<HypreADS> *amg = nullptr;
+
+   ParBilinearForm a_ho;
+   Array<int> ess_dofs;
+
+   const int ndofs;
+   double mdof;
+
+   Vector x, y;
+
+   RT_LORBench(int p, int requested_ndof, int dim, const std::string &name) :
+      mesh(MakeParCartesianMesh(p, requested_ndof, dim)),
+      fec_ho(p - 1, dim, BasisType::GaussLobatto, BasisType::IntegratedGLL),
+      fes_ho(&mesh, &fec_ho),
+      lor(fes_ho),
+      ads(fes_ho, lor.GetLORVertexCoordinates()),
+      a_ho(&fes_ho),
+      ndofs(fes_ho.GetTrueVSize()),
+      mdof(0.0)
+   {
+      fes_ho.GetBoundaryTrueDofs(ess_dofs);
+
+      a_ho.AddDomainIntegrator(new VectorFEMassIntegrator);
+      a_ho.AddDomainIntegrator(new DivDivIntegrator);
+
+      if (name == "RTAssembleBatched") { RTAssembleBatched(); }
+      if (name == "DiscreteCurl") { DiscreteCurl(); }
+      if (name == "ADSApply")
+      {
+         amg = new LORSolver<HypreADS>(a_ho, ess_dofs);
+         x.SetSize(ndofs);
+         y.SetSize(ndofs);
+         ADSApply();
+      }
+   }
+
+   void RTAssembleBatched()
+   {
+      NVTX("RTAssembleBatched");
+      MFEM_DEVICE_SYNC;
+
+      lor.AssembleWithoutBC(a_ho, A_lor);
+      A_lor.As<SparseMatrix>()->EliminateBC(ess_dofs,
+                        Operator::DiagonalPolicy::DIAG_KEEP);
+   }
+
+   void DiscreteCurl()
+   {
+      NVTX("DiscreteCurl");
+      MFEM_DEVICE_SYNC;
+
+      ads.FormCurlMatrix();
+   }
+
+   void ADSApply()
+   {
+      NVTX("ADSApply");
+      MFEM_DEVICE_SYNC;
+      amg->Mult(x, y);
+   }
+
+   ~RT_LORBench() { delete amg; }
+};
+
+struct ND_LORBench
+{
+   ParMesh mesh;
+   ND_FECollection fec_ho;
+   ParFiniteElementSpace fes_ho;
+
+   BatchedLORAssembly lor;
+   BatchedLOR_AMS ams;
+   OperatorHandle A_lor;
+
+   LORSolver<HypreAMS> *amg = nullptr;
+
+   ParBilinearForm a_ho;
+   Array<int> ess_dofs;
+
+   Vector evec, xvert;
+
+   const int ndofs;
+   double mdof;
+
+   Vector x, y;
+
+   ND_LORBench(int p, int requested_ndof, int dim, const std::string &name) :
+      mesh(MakeParCartesianMesh(p, requested_ndof, dim)),
+      fec_ho(p, dim, BasisType::GaussLobatto, BasisType::IntegratedGLL),
+      fes_ho(&mesh, &fec_ho),
+      lor(fes_ho),
+      ams(fes_ho, lor.GetLORVertexCoordinates()),
+      a_ho(&fes_ho),
+      ndofs(fes_ho.GetTrueVSize()),
+      mdof(0.0)
+   {
+      fes_ho.GetBoundaryTrueDofs(ess_dofs);
+
+      a_ho.AddDomainIntegrator(new VectorFEMassIntegrator);
+      a_ho.AddDomainIntegrator(new CurlCurlIntegrator);
+
+      if (name == "NDAssembleBatched") { NDAssembleBatched(); }
+      if (name == "DiscreteGradient") { DiscreteGradient(); }
+      if (name == "CoordinateVectors") { CoordinateVectors(); }
+      if (name == "VertexCoordinates") { VertexCoordinates(); }
+      if (name == "AMSApply")
+      {
+         amg = new LORSolver<HypreAMS>(a_ho, ess_dofs);
+         x.SetSize(ndofs);
+         y.SetSize(ndofs);
+         AMSApply();
+      }
+   }
+
+   void NDAssembleBatched()
+   {
+      NVTX("NDAssembleBatched");
+      MFEM_DEVICE_SYNC;
+
+      lor.AssembleWithoutBC(a_ho, A_lor);
+      A_lor.As<SparseMatrix>()->EliminateBC(ess_dofs,
+                        Operator::DiagonalPolicy::DIAG_KEEP);
+   }
+
+   void DiscreteGradient()
+   {
+      NVTX("DiscreteGradient");
+      MFEM_DEVICE_SYNC;
+
+      ams.FormGradientMatrixLocal();
+   }
+
+   void CoordinateVectors()
+   {
+      NVTX("CoordinateVectors");
+      MFEM_DEVICE_SYNC;
+
+      ams.FormCoordinateVectors(lor.GetLORVertexCoordinates());
+   }
+
+   void VertexCoordinates()
+   {
+      NVTX("VertexCoordinates");
+      MFEM_DEVICE_SYNC;
+
+      BatchedLORAssembly::FormLORVertexCoordinates(fes_ho, xvert, &evec);
+   }
+
+   void AMSApply()
+   {
+      NVTX("AMSApply");
+      MFEM_DEVICE_SYNC;
+      amg->Mult(x, y);
+   }
+
+   ~ND_LORBench() { delete amg; }
+};
+
+struct LORBench
+{
+   Mesh mesh;
+
+   H1_FECollection fec_ho;
+   FiniteElementSpace fes_ho;
+
+   IntegrationRules irs;
+   const IntegrationRule &ir;
+
+   LORDiscretization lor;
+   BatchedLORAssembly lor_b;
+   BilinearForm a_ho, a_lor;
+
+   OperatorHandle A_ho, A_lor;
+
+   HYPRE_Int row_starts[2];
+   HypreParMatrix *A = nullptr;
+   HypreBoomerAMG amg;
+
+   Array<int> ess_dofs;
+
+   const int ndofs;
+   double mdof;
+   Vector x, y;
+
+   LORBench(int p, int requested_ndof, int dim, const std::string &name) :
+      mesh(MakeCartesianMesh(p, requested_ndof, dim)),
+      fec_ho(p, dim),
+      fes_ho(&mesh, &fec_ho),
+      irs(0, Quadrature1D::GaussLobatto),
+      ir(irs.Get(mesh.GetElementGeometry(0), 1)),
+      lor(fes_ho),
+      lor_b(fes_ho),
+      a_ho(&fes_ho),
+      a_lor(&lor.GetFESpace()),
+      ndofs(fes_ho.GetTrueVSize()),
+      mdof(0.0)
+   {
+      // std::cout << "Requested ndof: "
+      //           << std::setw(10) << requested_ndof
+      //           << " Actual: "
+      //           << std::setw(10) << ndofs << '\n';
+      fes_ho.GetBoundaryTrueDofs(ess_dofs);
+
+      a_ho.AddDomainIntegrator(new DiffusionIntegrator);
+      a_ho.AddDomainIntegrator(new MassIntegrator);
+      a_ho.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+
+      a_lor.AddDomainIntegrator(new DiffusionIntegrator(&ir));
+      a_lor.AddDomainIntegrator(new MassIntegrator(&ir));
+      a_lor.SetAssemblyLevel(AssemblyLevel::FULL);
+
+      if (name == "ApplyHO" || name =="Vcycle")
+      {
+         x.SetSize(ndofs);
+         y.SetSize(ndofs);
+         x.Randomize(1);
+         y.Randomize(2);
+      }
+
+      // warm up
+      if (name == "AssembleHO" || name == "ApplyHO") { AssembleHO(); }
+      if (name == "ApplyHO") { ApplyHO(); }
+      if (name == "AssembleBatched") { AssembleBatched(); }
+      if (name == "AssembleFull") { AssembleFull(); }
+      if (name == "AMGSetup" || name == "Vcycle")
+      {
+         AssembleBatched();
+         SparseMatrix &A_serial = lor.GetAssembledMatrix();
+         row_starts[0] = 0;
+         row_starts[1] = A_serial.Height();
+         A = new HypreParMatrix(MPI_COMM_WORLD, A_serial.Height(), row_starts, &A_serial);
+         amg.SetOperator(*A);
+         amg.SetPrintLevel(0);
+      }
+      if (name == "Vcycle") { amg.Setup(x,y); }
+   }
+
+   void AssembleHO()
+   {
+      NVTX("AssembleHO");
+      MFEM_DEVICE_SYNC;
+      a_ho.Assemble();
+      a_ho.FormSystemMatrix(ess_dofs, A_ho);
+   }
+
+   void AssembleFull()
+   {
+      NVTX("AssembleFull");
+      MFEM_DEVICE_SYNC;
+      a_lor.Assemble();
+      a_lor.FormSystemMatrix(ess_dofs, A_lor);
+   }
+
+   void AssembleBatched()
+   {
+      NVTX("AssembleBatched");
+      MFEM_DEVICE_SYNC;
+      // lor.AssembleSystem(a_ho, ess_dofs);
+      lor_b.Assemble(a_ho, ess_dofs, A_lor);
+   }
+
+   void ApplyHO()
+   {
+      NVTX("ApplyHO");
+      MFEM_DEVICE_SYNC;
+      A_ho->Mult(x, y);
+   }
+
+   void AMGSetup()
+   {
+      NVTX("AMG Setup");
+      MFEM_DEVICE_SYNC;
+      amg.SetOperator(*A);
+      amg.Setup(x, y);
+   }
+
+   void Vcycle()
+   {
+      NVTX("Vcycle");
+      MFEM_DEVICE_SYNC;
+      amg.Mult(x, y);
+   }
+
+   ~LORBench() { delete A; }
+};
+
+// The different orders the tests can run
+#define P_ORDERS bm::CreateDenseRange(1,7,1)
+
+// The different sides of the mesh
+#define LOG_NDOFS bm::CreateDenseRange(7,23,1)
+
+// Dimensions: 2 or 3
+#define DIMS bm::CreateDenseRange(2,3,1)
+
+/// Kernels definitions and registrations
+#define Benchmark(Class, Name)\
+static void Name(bm::State &state){\
+   const int p = state.range(0);\
+   const int log_ndof = state.range(1);\
+   const int requested_ndof = pow(2, log_ndof);\
+   const int dim = state.range(2);\
+   const std::string name = #Name;\
+   if (p == 1 && log_ndof >= 21) { state.SkipWithError("Problem size"); return; }\
+   if (p == 2 && log_ndof >= 23) { state.SkipWithError("Problem size"); return; }\
+   if (p == 3 && log_ndof >= 23) { state.SkipWithError("Problem size"); return; }\
+   if ((name == "ADSApply" || name == "AMSApply") && log_ndof >= 21)\
+   { state.SkipWithError("Problem size"); return; }\
+   Class lor(p, requested_ndof, dim, name);\
+   while (state.KeepRunning()) { lor.Name(); }\
+   bm::Counter::Flags flags = bm::Counter::kIsIterationInvariantRate;\
+   state.counters["MDof/s"] = bm::Counter(1e-6*lor.ndofs, flags);\
+   state.counters["dofs"] = bm::Counter(lor.ndofs);\
+   state.counters["p"] = bm::Counter(p);\
+}\
+BENCHMARK(Name)\
+            -> ArgsProduct({P_ORDERS, LOG_NDOFS, DIMS})\
+            -> Unit(bm::kMillisecond)\
+            -> Iterations(10);
+
+Benchmark(LORBench, AssembleHO)
+Benchmark(LORBench, AssembleFull)
+Benchmark(LORBench, AssembleBatched)
+
+Benchmark(LORBench, ApplyHO)
+Benchmark(LORBench, AMGSetup)
+Benchmark(LORBench, Vcycle)
+
+Benchmark(ND_LORBench, NDAssembleBatched)
+Benchmark(ND_LORBench, DiscreteGradient)
+Benchmark(ND_LORBench, CoordinateVectors)
+Benchmark(ND_LORBench, VertexCoordinates)
+Benchmark(ND_LORBench, AMSApply)
+
+Benchmark(RT_LORBench, RTAssembleBatched)
+Benchmark(RT_LORBench, DiscreteCurl)
+Benchmark(RT_LORBench, ADSApply)
+
+int main(int argc, char *argv[])
+{
+   Mpi::Init();
+
+   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/makefile

@@ -19,8 +19,9 @@ MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 
 SEQ_TESTS = bench_assembly_levels bench_ceed bench_dg_amr bench_elasticity \
-            bench_tmop bench_vector bench_virtuals
+            bench_tmop bench_vector bench_virtuals bench_lor
 PAR_TESTS = 
+
 ifeq ($(MFEM_USE_MPI),NO)
    TESTS = $(SEQ_TESTS)
 else
@@ -60,4 +61,3 @@ clean-build:
 
 clean-exec:
 	@rm -f refined.mesh sol.gf
-

tests/unit/fem/test_lor_batched.cpp

@@ -89,7 +89,8 @@ RT_FECollection *NewLOR_FE_Collection<RT_FECollection>(int order, int dim)
 template <typename FE_COLL, typename INTEG_1, typename INTEG_2>
 void TestBatchedLOR()
 {
-   const int order = 5;
+   const bool all_tests = launch_all_non_regression_tests;
+   const int order = !all_tests ? 5 : GENERATE(1,2,3,5);
    const auto mesh_fname = GENERATE(
                               "../../data/star-q3.mesh",
                               "../../data/fichera-q3.mesh"