Finished DeviceMult for H1Space

examples/CMakeLists.txt

@@ -10,6 +10,9 @@
 # CONTRIBUTING.md for details.
 
 list(APPEND ALL_EXE_SRCS
+  blas.cpp
+  blas-ex.cpp
+  blas-ex2.cpp
   ex0.cpp
   ex1.cpp
   ex2.cpp
@@ -49,6 +52,9 @@ list(APPEND ALL_EXE_SRCS
 
 if (MFEM_USE_MPI)
   list(APPEND ALL_EXE_SRCS
+    blas.cpp
+    blas-ex.cpp
+    blas-ex.cpp
     ex0p.cpp
     ex1p.cpp
     ex2p.cpp

examples/blas.cpp

@@ -0,0 +1,225 @@
+//                      MFEM cuBLAS and hipBLAS gemvBatched() Example
+//
+// Compile with: make blas
+//
+// Device sample runs:
+//               blas -d cuda
+//               blas -d cuda -n 2 -ne 2 
+//               blas -pa -d cuda
+//               blas -fa -d cuda
+//               blas -d hip 
+//               blas -d hip -pa 
+//               blas -d hip -fa
+//               blas -d hip
+
+// Description:  This example code demonstrates the use of cu or hipBLAS on MFEM 
+//               objects to multiply batched square matrices with batched vectors. 
+//               It utilizes the BLAS functions 
+//                  cublasDgemvBatched()  or  hipblasDgemvBatched().
+//               Note that version cuda/11.7.0 or newer is needed. 
+// 
+//               User can specify the number of rows (cols) n for the (nxn) square 
+//               square matrices stored in DenseTensor A and (nx1) vectors stored in 
+//               Vector X, as well as scalars alpha and beta to transform the operations: 
+//                  alpha * A[i] * X[i] + beta = Y[i]. 
+//               Output also solves the linear system Ab = x as needed: 
+//                  Ainv[i] * X[i] = B[i]. 
+//               This function is specifically included to demonstrate the Batch Linear
+//               Algebra functions ported with cu/hipBLAS (see batchlinalg.cpp).
+
+#include "mfem.hpp"
+#include "../fem/bilinearform.cpp" // TODO: causes Seg Fault when removed (Note: not included in batchlinalg.cpp so Seg Fault occurs when running inverse portion)
+#include <fstream>
+#include <iostream>
+#define IDXT(i,j,k,ld) (((ld)*(ld)*(k))+((j)*(ld))+(i))
+#define IDXM(i,j,ld) ((ld*j)+i)
+
+using namespace std;
+using namespace mfem;
+
+int main (int argc, char *argv[]){
+    // 1. Parse command-line options.
+    int ndof = 4;
+    int num_elem = 3;
+    double alpha = 1.;
+    double beta = 0.;
+    bool inverse = true;
+    bool pa = false;
+    bool fa = false;
+    const char *device_config = "cuda";
+
+    OptionsParser args(argc, argv);
+    args.AddOption(&ndof, "-n", "--ndof",
+                    "Number of Rows/Columns of batched square matrices.");
+    args.AddOption(&num_elem, "-ne", "--num_elem",
+                    "Number of matrices in tensor (batch size).");
+    args.AddOption(&alpha, "-a", "--alpha",
+                    "Scalar for batched matrices.");
+    args.AddOption(&beta, "-b", "--beta",
+                    "Scalar to add to multiplication.");
+    args.AddOption(&inverse, "-inv", "--inverse", "-no-inv", 
+                    "--no-inverse", "Disable inverse of batched matrices.");
+    args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
+                    "--no-partial-assembly", "Enable Partial Assembly.");
+    args.AddOption(&fa, "-fa", "--full-assembly", "-no-fa",
+                    "--no-full-assembly", "Enable Full Assembly.");
+    args.AddOption(&device_config, "-d", "--device",
+                    "Device configuration string, see Device::Configure().");
+    args.Parse();
+    if (!args.Good())
+    {
+        args.PrintUsage(cout);
+        return 1;
+    }
+    args.PrintOptions(cout);
+
+    // 2. Specify hardware devices such as GPUs.
+    Device device(device_config);
+    device.Print();
+    printf ("\n");
+
+    // 3. Create a DenseTensor A and Vector X to batch multiply. Create a Vector Y container for the result. 
+    DenseTensor A(ndof,ndof,num_elem);
+    Vector X(ndof*num_elem);
+    Vector Y(ndof*num_elem);
+    printf ("A transpose is \n"); // col-major, so prints vectors of columns of A out together in each row
+    for (int k = 0; k < num_elem; k++) {
+        for (int j = 0; j < ndof; j++) {
+            X.HostReadWrite()[IDXM(j,k,ndof)] = (double) (IDXM(j,k,ndof));
+            for (int i = 0; i < ndof; i++) {
+                if (i==j) {A.HostReadWrite()[IDXT(i,j,k,ndof)] = 1;}
+                else {A.HostReadWrite()[IDXT(i,j,k,ndof)] = (double) (IDXT(j,i,k,ndof));}
+                printf ("%9.3f", A.HostRead()[IDXT(i,j,k,ndof)]);
+            }
+            printf ("\n");  
+        }
+        printf ("\n");
+    }
+    printf ("X transpose is \n");
+    for (int k = 0; k < num_elem; k++) {
+        for (int j = 0; j < ndof; j++) {
+            printf ("%9.3f", X.HostRead()[IDXM(j,k,ndof)]);
+        }
+        printf ("\n");
+    }
+    printf ("\n");
+
+    // 4. Run CUDA or HIP
+    //    Full commands not written out; using MFEM_USE_CUDA_OR_HIP implementation.
+    //    See commented section below for example in cuBLAS and hipBLAS.
+    
+    // gemvBatched() requires array of pointers for each A, X, Y 
+    Array<double *>devPtrA(num_elem);  
+    Array<double *>devPtrX(num_elem);  
+    Array<double *>devPtrY(num_elem);  
+
+    for (int k = 0; k < num_elem; k++) {
+        devPtrA[k] = &A.ReadWrite()[ndof*ndof*k];
+        devPtrX[k] = &X.ReadWrite()[ndof*k];
+        devPtrY[k] = &Y.ReadWrite()[ndof*k];
+    }
+    
+    MFEM_cu_or_hip(blasStatus_t) stat;  // collects generation of cublasStatus_t
+    MFEM_cu_or_hip(blasHandle_t) handle;  // tracks handle into API; can be specified futher but NULL works too
+
+    stat = MFEM_cu_or_hip(blasCreate)(&handle);  // create handle to start CUBLAS work on the device; i.e. initialize CUBLAS
+    if (stat != MFEM_CU_or_HIP(BLAS_STATUS_SUCCESS)) {
+        printf ("BLAS initialization failed\n");
+        return EXIT_FAILURE;
+    }
+
+    stat = MFEM_cu_or_hip(blasDgemvBatched) (handle, MFEM_CU_or_HIP(BLAS_OP_N), ndof, ndof,
+                                &alpha, devPtrA.Read(), ndof, devPtrX.Read(), 1,
+                                &beta, devPtrY.ReadWrite(), 1, num_elem);  // version 11.9.0 needs batchCount = W as the last parameter
+    if (stat != MFEM_CU_or_HIP(BLAS_STATUS_SUCCESS)) {
+        printf ("BLAS gemvBatched() failed\n");
+        printf("\n");
+        return EXIT_FAILURE;
+    }
+
+    MFEM_cu_or_hip(blasDestroy)(handle);  // end API stream
+    
+
+    // 5. Done! We can now output.
+    printf ("Y tranpose is \n");
+    for (int k = 0; k < num_elem; k++) {
+        for (int j = 0; j < ndof; j++) {
+            printf ("%9.3f", Y.HostRead()[IDXM(j,k,ndof)]);  
+        }
+        printf ("\n");
+    }
+    printf ("\n");
+
+
+    // 6. If wanted, solve the linear systems
+    if (inverse) {
+        Vector B(ndof*num_elem);
+        DenseTensor Ainv(ndof,ndof,num_elem);
+        BatchSolver batchSolver(BatchSolver::SolveMode::INVERSE);
+        batchSolver.AssignMatrices(A);
+        batchSolver.GetInverse(Ainv);
+        batchSolver.Mult(X, B);
+
+        printf ("A inverse transpose is \n");
+        for (int k = 0; k < num_elem; k++) {
+            for (int j = 0; j < ndof; j++) {
+                for (int i = 0; i < ndof; i++) {
+                    printf ("%9.3f", Ainv.HostRead()[IDXT(i,j,k,ndof)]);
+                }
+                printf ("\n");  
+            }
+            printf ("\n");
+        }
+
+        printf ("B tranpose is \n");
+        for (int k = 0; k < num_elem; k++) {
+            for (int j = 0; j < ndof; j++) {
+                printf ("%9.3f", B.HostRead()[IDXM(j,k,ndof)]);  // col-major, so prints vectors of columns out together in each row
+            }
+            printf ("\n");
+        }
+        printf ("\n");
+    }
+
+    return EXIT_SUCCESS;
+
+
+    // *****CUBLAS IMPLEMENTATIONS*****
+    // cublasStatus_t stat;  // collects generation of cublasStatus_t
+    // cublasHandle_t handle;  // tracks handle into API; can be specified futher but NULL works too
+    // stat = cublasCreate(&handle);  // create handle to start CUBLAS work on the device; i.e. initialize CUBLAS
+    // if (stat != CUBLAS_STATUS_SUCCESS) {
+    //     printf ("CUBLAS initialization failed\n");
+    //     return EXIT_FAILURE;
+    // }
+    // stat = cublasDgemvBatched (handle, CUBLAS_OP_N, M, N,
+    //                             &alpha, devPtrA.READ(), M, devPtrX.Read(), 1,
+    //                             &beta, devPtrY.ReadWrite(), 1, W);  // version 11.9.0 needs batchCount = W as the last parameter
+    // if (stat != CUBLAS_STATUS_SUCCESS) {
+    //     printf ("CUBLAS gemvBatched() failed\n");
+    //     printf("%s", cublasGetStatusString(stat));
+    //     printf("\n");
+    //     return EXIT_FAILURE;
+    // }
+    // cublasDestroy(handle);  // end API stream
+
+
+    // *****HIPBLAS IMPLEMENTATIONS*****
+    // hipblasStatus_t stat;  // collects generation of hipblasStatus_t
+    // hipblasHandle_t handle = nullptr;  // tracks handle into API; can be specified futher but NULL works too
+    // stat = hipblasCreate(&handle);  // create handle to start HIPBLAS work on the device; i.e. initialize HIPBLAS
+    // if (stat != HIPBLAS_STATUS_SUCCESS) {
+    //     printf ("HIPBLAS initialization failed\n");
+    //     return EXIT_FAILURE;
+    // }
+    // stat = hipblasDgemvBatched (handle, HIPBLAS_OP_N, M, N,
+    //                             &alpha, devPtrA.Read(), M, devPtrX.Read(), 1,
+    //                             &beta, devPtrY.ReadWrite(), 1, W);  
+    // if (stat != HIPBLAS_STATUS_SUCCESS) {
+    //     printf ("HIPBLAS gemvBatched() failed\n");
+    //     printf("%s", hipblasStatusToString(stat));
+    //     printf("\n");
+    //     return EXIT_FAILURE;
+    // }
+    // hipblasDestroy(handle);  // end API stream
+}

examples/ex0.cpp

@@ -13,12 +13,19 @@
 //              can be specified by command line options.
 
 #include "mfem.hpp"
+#include "../linalg/dtensor.hpp"
 #include <fstream>
 #include <iostream>
+#define IDX2C(i,j,k,inc) ((i)+(j*inc)+(k*inc*inc))
+#define IDXV(i,j,inc) ((i)+(j*inc))
 
 using namespace std;
 using namespace mfem;
 
+/// @brief 
+/// @param argc 
+/// @param argv 
+/// @return 
 int main(int argc, char *argv[])
 {
    // 1. Parse command line options.
@@ -77,5 +84,111 @@ int main(int argc, char *argv[])
    x.Save("sol.gf");
    mesh.Save("mesh.mesh");
 
+   cout << "Let's see if this works \n";
+
+   ConstantCoefficient q(111.0);
+   Coefficient *Q(&q); 
+   Q->SetTime(222.);
+   real_t time = Q->GetTime();
+
+
+   cout << "Pointer Q is " << Q << "\n";
+   cout << "This matches with the address of q, which is " << &q << "\n";
+   // cout << "Note that the value that Q points to is " << *Q.constant << "\n";
+   cout << "This matches with the value of q, which is " << q.constant << endl; 
+   cout << "Pointer Q has the time of " << time << "\n";
+
+
+
+   // Check operators are not issues
+
+   // QuadratureFunctionCoefficient q();
+   const real_t detJ = 3.141;
+   const real_t d_D = detJ * q.constant;  // auto d_D is a DeviceTensor made up of real_t elements; ERRORS
+
+   // cout << "d_D = " << d_D << endl;
+
+
+   // Check pointer behavior: 
+   double *f = new double[4];
+   double *F = f;
+   cout << "Size of F is " << sizeof (*F) << endl;
+
+   int i, j, k;
+   for (k=0;k<2;k++) {
+      for (j=0;j<6;j++) {
+         for (i=0;i<5;i++) {
+            cout << IDX2C(i,j,k,6) << " ";
+         }
+      }
+   }
+   cout << endl;
+
+   DenseTensor G(6,6,2);
+   Vector g(6);
+   Vector gy(6);
+
+   cout << "At first, G is " << endl;
+   for (k = 0; k < 2; k++) {
+      for (j = 0; j < 6; j++) {
+         for (i = 0; i < 6; i++) {
+            printf ("%7.0f", G.Data()[IDX2C(i,j,k,6)]);
+         }
+         printf("\n");
+      }
+      printf("\n");
+   }
+
+   for (k = 0; k < 2; k++) {
+      for (j = 0; j < 6; j++) {
+         for (i = 0; i < 6; i++) {
+            g.GetData()[i] = i;
+            G.Data()[IDX2C(i,j,k,6)] = (double)(IDX2C(i,j,k,6));
+         }
+         printf("\n");
+      }
+      printf("\n");
+   }
+   g.Print();
+   cout << g.GetData() << endl;
+
+   cout << "G is " << endl;
+   for (k = 0; k < 2; k++) {
+      for (j = 0; j < 6; j++) {
+         for (i = 0; i < 6; i++) {
+            printf ("%7.0f", G.Data()[IDX2C(i,j,k,6)]);
+         }
+         printf("\n");
+      }
+      printf("\n");
+   }
+   cout << "G at (6,6,2) is " << G.Data()[IDX2C(5,5,1,6)] << endl;
+   cout << &G.Data()[IDX2C(5,5,1,6)] << endl;
+
+   auto d_G = Reshape(G.Read(), 6, 6, 1);
+
+   gy.Print();
+   cout << endl;
+
+
+   double* y = 0;  // host device pointer
+   cout << "Size of y: " << sizeof(y) << endl;
+   y = (double *)malloc (6 * 1 * sizeof(*y));
+
+   for (j = 0; j < 6; j++) {
+      cout << y[j] << " ";
+   }
+   cout << endl;
+
+   Vector Y(y,6);
+   Y.SetData(y);
+   cout << "Vector Y is "; Y.Print();
+
+   cout << "Size of G: " << sizeof(*G.Data()) << endl;
+   cout << "Size of g: " << sizeof(g.GetData()) << endl;
+   cout << "Size of gy: " << sizeof(gy.GetData()) << endl;
+   cout << "Size of y: " << sizeof(y) << endl;
+
+
    return 0;
 }

fem/coefficient.cpp

@@ -1704,7 +1704,7 @@ void CoefficientVector::Project(Coefficient &coeff)
    {
       SetConstant(const_coeff->constant);
    }
-   else if (auto *qf_coeff = dynamic_cast<QuadratureFunctionCoefficient*>(&coeff))
+   else if (auto *qf_coeff = dynamic_cast<QuadratureFunctionCoefficient*>(&coeff)) 
    {
       MakeRef(qf_coeff->GetQuadFunction());
    }
@@ -1773,11 +1773,12 @@ void CoefficientVector::ProjectTranspose(MatrixCoefficient &coeff)
 void CoefficientVector::MakeRef(const QuadratureFunction &qf_)
 {
    vdim = qf_.GetVDim();
+   // cout << "Check this is 1: " << qf_.GetVDim() << endl;
    const QuadratureSpaceBase *qs2 = qf_.GetSpace();
    MFEM_CONTRACT_VAR(qs2); // qs2 used only for asserts
    MFEM_VERIFY(qs2 != NULL, "Invalid QuadratureSpace.")
    MFEM_VERIFY(qs2->GetMesh() == qs.GetMesh(), "Meshes differ.");
-   MFEM_VERIFY(qs2->GetOrder() == qs.GetOrder(), "Orders differ.");
+   MFEM_VERIFY(qs2->GetOrder() == qs.GetOrder(), "Orders differ. qs2 is "<<qs2->GetOrder()<<" but qs is "<<qs.GetOrder());
    Vector::MakeRef(const_cast<QuadratureFunction&>(qf_), 0, qf_.Size());
 }
 

fem/transfer.cpp

@@ -325,16 +325,16 @@ void L2ProjectionGridTransfer::L2Projection::ElemMixedMass(
 }
 
 L2ProjectionGridTransfer::L2ProjectionL2Space::L2ProjectionL2Space
-(const FiniteElementSpace &fes_ho_, const FiniteElementSpace &fes_lor_,
+(const FiniteElementSpace &fes_ho_, const FiniteElementSpace &fes_lor_, Coefficient *coeff_,
  const bool use_device_, const bool verify_solution_, MemoryType d_mt)
-   : L2Projection(fes_ho_, fes_lor_),
+   : L2Projection(fes_ho_, fes_lor_), 
      use_device(use_device_), verify_solution(verify_solution_),
-     d_mt_(d_mt)
+     d_mt_(d_mt), coeff(coeff_)
 {
 
    if (use_device)
    {
-      DeviceL2ProjectionL2Space(fes_ho_, fes_lor_);
+      DeviceL2ProjectionL2Space(fes_ho_, fes_lor_, coeff_);
       if (!verify_solution) {return;}
    }
 
@@ -478,10 +478,24 @@ L2ProjectionGridTransfer::L2ProjectionL2Space::L2ProjectionL2Space
    }
 }
 
-void L2ProjectionGridTransfer::L2ProjectionL2Space::DeviceL2ProjectionL2Space
-(const FiniteElementSpace &fes_ho_, const FiniteElementSpace &fes_lor_)
-{
 
+
+// START OF DEVICE IMPLEMENTATION
+
+void L2ProjectionGridTransfer::L2ProjectionL2Space::DeviceL2ProjectionL2Space
+(const FiniteElementSpace &fes_ho_, const FiniteElementSpace &fes_lor_, Coefficient *coeff_) // coeff : coefficient for weighted integration
+{
+   // dynamic_cast to check if QuadFuncCoeff; if yes, continue; if no, return error
+   auto qfunc_coeff = dynamic_cast<QuadratureFunctionCoefficient*>(coeff_);
+   if (qfunc_coeff == NULL)
+   {
+      mfem_error("Not a QuadratureFunctionCoefficient");
+   }
+
+   QuadratureFunction qfunc = qfunc_coeff->GetQuadFunction();
+   // Store the mixed mass matrix integration rule, which is assumed same on all elements
+   const IntegrationRule ir = qfunc.GetIntRule(0); 
+   
    Mesh *mesh_ho = fes_ho.GetMesh();
    Mesh *mesh_lor = fes_lor.GetMesh();
    int nel_ho = mesh_ho->GetNE();
@@ -567,6 +581,7 @@ void L2ProjectionGridTransfer::L2ProjectionL2Space::DeviceL2ProjectionL2Space
          const int dim = mesh_ho->Dimension();
 
          MFEM_ASSERT(nel_ho*nref == nel_lor, "we expect nel_ho*nref == nel_lor");
+         MFEM_VERIFY(D.TotalSize() == qfunc.Size(), "Dimensions don't match  "<<D.TotalSize()<<" "<<qfunc.Size()); 
 
          //*********************************
          // Setup data at quadrature points
@@ -574,9 +589,10 @@ void L2ProjectionGridTransfer::L2ProjectionL2Space::DeviceL2ProjectionL2Space
          if (dim == 1)
          {
 
-            const auto W = Reshape(ir->GetWeights().Read(), Q1D);
-            const auto J = Reshape(geo_facts->detJ.Read(), Q1D, nel_lor);
-            const auto d_D = Reshape(D.Write(), qPts, nref, nel_ho);
+            const auto W = Reshape(ir->GetWeights().Read(), Q1D);  // grabbing the weights of the integration rule 
+            const auto J = Reshape(geo_facts->detJ.Read(), Q1D, nel_lor);  // 
+            const auto d_D = Reshape(D.Write(), qPts, nref, nel_ho);  // diagonal at the quadrature points
+            const auto d_qfunc = Reshape(qfunc.Read(), qPts, nref, nel_ho);
 
             mfem::forall(nel_ho, [=] MFEM_HOST_DEVICE (int iho)
             {
@@ -585,22 +601,23 @@ void L2ProjectionGridTransfer::L2ProjectionL2Space::DeviceL2ProjectionL2Space
                   const int lo_el_id = iref + nref*iho;
                   for (int qx=0; qx<Q1D; ++qx)
                   {
-
                      const real_t detJ = J(qx, lo_el_id);
-                     d_D(qx, iref, iho) = W(qx) * detJ;
+                     d_D(qx, iref, iho) = W(qx) * detJ * d_qfunc(qx, iref, iho); 
                   }
                }
             });
 
          }
 
-         if (dim == 2)
+         if (dim == 2)  
          {
 
             const auto W = Reshape(ir->GetWeights().Read(), Q1D, Q1D);
             const auto J = Reshape(geo_facts->detJ.Read(), Q1D,Q1D, nel_lor);
             const auto d_D = Reshape(D.Write(), qPts, nref, nel_ho);
 
+            const auto d_qfunc = Reshape(qfunc.Read(), qPts, nref, nel_ho);
+
             mfem::forall(nel_ho, [=] MFEM_HOST_DEVICE (int iho)
             {
                for (int iref = 0; iref < nref; ++iref)
@@ -613,7 +630,7 @@ void L2ProjectionGridTransfer::L2ProjectionL2Space::DeviceL2ProjectionL2Space
                      {
                         const int q = qx + Q1D*qy;
                         const real_t detJ = J(qx, qy, lo_el_id);
-                        d_D(q, iref, iho) = W(qx, qy) * detJ;
+                        d_D(q, iref, iho) = W(qx, qy) * detJ * d_qfunc(q, iref, iho);
                      }
                   }
 
@@ -627,6 +644,7 @@ void L2ProjectionGridTransfer::L2ProjectionL2Space::DeviceL2ProjectionL2Space
             const auto W = Reshape(ir->GetWeights().Read(), Q1D, Q1D, Q1D);
             const auto J = Reshape(geo_facts->detJ.Read(), Q1D, Q1D, Q1D, nel_lor);
             const auto d_D = Reshape(D.Write(), qPts, nref, nel_ho);
+            const auto d_qfunc = Reshape(qfunc.Read(), qPts, nref, nel_ho);
 
             mfem::forall(nel_ho, [=] MFEM_HOST_DEVICE (int iho)
             {
@@ -641,10 +659,9 @@ void L2ProjectionGridTransfer::L2ProjectionL2Space::DeviceL2ProjectionL2Space
                      {
                         for (int qx=0; qx<Q1D; ++qx)
                         {
-
                            const int q = qx + Q1D*qy + Q1D*Q1D*qz;
                            const real_t detJ = J(qx, qy, qz, lo_el_id);
-                           d_D(q, iref, iho) = W(qx, qy, qz) * detJ;
+                           d_D(q, iref, iho) = W(qx, qy, qz) * detJ * d_qfunc(q, iref, iho); 
                         }
                      }
                   }
@@ -683,9 +700,8 @@ void L2ProjectionGridTransfer::L2ProjectionL2Space::DeviceL2ProjectionL2Space
    } //Competed setup of basis function and quadrature point
 
    //Assemble mixed mass matrix
-   Vector M_mixed_all;
-   Vector M_ea_lor;
-   int ndof_lor;
+   Vector M_ea_lor;  // declare lor mass matrix 
+   int ndof_lor;  
    int ndof_ho;
    int nref;
    {
@@ -746,9 +762,11 @@ void L2ProjectionGridTransfer::L2ProjectionL2Space::DeviceL2ProjectionL2Space
    } // end of mixed assembly mass matrix
 
 
+
    //R = inv(M_L) * M_mixed
-   //Need to compute M_L
-   MassIntegrator mi;
+   //Need to compute M_L. 
+   //Note: Using user-inputted M_LH IntegrationRule ir (higher order than needed) in order to re-use coeff 
+   MassIntegrator mi(*coeff, &ir); 
 
    const bool add = false;
    mi.AssembleEA(fes_lor, M_ea_lor, add);
@@ -1006,7 +1024,7 @@ void L2ProjectionGridTransfer::L2ProjectionL2Space::MultTranspose(
 
    if (use_device)
    {
-      DeviceMult(x,y);
+      DeviceMultTranspose(x,y);
       if (!verify_solution) {return;}
    }
 
@@ -1282,11 +1300,30 @@ void L2ProjectionGridTransfer::L2ProjectionL2Space::DeviceProlongateTranspose(
 }
 
 L2ProjectionGridTransfer::L2ProjectionH1Space::L2ProjectionH1Space(
-   const FiniteElementSpace& fes_ho_, const FiniteElementSpace& fes_lor_)
-   : L2Projection(fes_ho_, fes_lor_)
+   const FiniteElementSpace& fes_ho_, const FiniteElementSpace& fes_lor_, Coefficient* coeff_, 
+   const bool use_device_, const bool verify_solution_, MemoryType d_mt)
+   : L2Projection(fes_ho_, fes_lor_), 
+     use_device(use_device_), verify_solution(verify_solution_), 
+     d_mt_(d_mt), coeff(coeff_)
 {
+   if (use_device)
+   {
+      DeviceL2ProjectionH1Space(fes_ho_, fes_lor_, coeff_);
+      if (!verify_solution) {return;}
+   }
+
    std::unique_ptr<SparseMatrix> R_mat, M_LH_mat;
-   std::tie(R_mat, M_LH_mat) = ComputeSparseRAndM_LH();
+   bool GetM_LHError = false;
+   bool GetML_invError = false;
+   if (verify_solution)
+   {
+      // Check that M_LH is built correctly
+      GetM_LHError = true;
+
+      // Check that lumped(inv M_L) is built correctly 
+      GetML_invError = true;
+   }
+   std::tie(R_mat, M_LH_mat) = ComputeSparseRAndM_LH(GetM_LHError, GetML_invError);
 
    FiniteElementSpace fes_ho_scalar(fes_ho.GetMesh(), fes_ho.FEColl(), 1);
    FiniteElementSpace fes_lor_scalar(fes_lor.GetMesh(), fes_lor.FEColl(), 1);
@@ -1376,20 +1413,382 @@ L2ProjectionGridTransfer::L2ProjectionH1Space::L2ProjectionH1Space(
 
 void L2ProjectionGridTransfer::L2ProjectionH1Space::SetupPCG()
 {
-   // Basic PCG solver setup
-   pcg.SetPrintLevel(0);
-   // pcg.SetPrintLevel(IterativeSolver::PrintLevel().Summary());
-   pcg.SetMaxIter(1000);
-   // initial values for relative and absolute tolerance
-   pcg.SetRelTol(1e-13);
-   pcg.SetAbsTol(1e-13);
-   pcg.SetPreconditioner(*precon);
-   pcg.SetOperator(*RTxM_LH);
+    // Basic PCG solver setup
+    pcg.SetPrintLevel(0);
+    // pcg.SetPrintLevel(IterativeSolver::PrintLevel().Summary());
+    pcg.SetMaxIter(1000);
+    // initial values for relative and absolute tolerance
+    pcg.SetRelTol(1e-13);
+    pcg.SetAbsTol(1e-13);
+    pcg.SetPreconditioner(*precon);
+    pcg.SetOperator(*RTxM_LH);
+}
+
+void L2ProjectionGridTransfer::L2ProjectionH1Space::DeviceL2ProjectionH1Space(
+   const FiniteElementSpace &fes_ho_, const FiniteElementSpace &fes_lor_, Coefficient *coeff_)
+{
+   // TODO: implement matrix-free action of R 
+   // There are 4 parts: 
+   //  ( )  inv( lumped(M_L) ), which is a diagonal matrix (essentially a vector). Need to add QuadCoeffFunction ability.
+   //  (x)  ElementRestrictionOperator for LOR space, need to apply the transpose
+   //  (x)  M_{LH}, with QuadCoeffFunction ability
+   //  (x)  ElementRestrictionOperator for HO space
+
+   Mesh* mesh_ho = fes_ho.GetMesh();
+   Mesh* mesh_lor = fes_lor.GetMesh();
+   int nel_ho = mesh_ho->GetNE();
+   int nel_lor = mesh_lor->GetNE();
+   int ndof_ho = fes_ho.GetNDofs();
+   int ndof_lor = fes_lor.GetNDofs(); 
+
+   // If the local mesh is empty, skip all computations
+   if (nel_ho == 0)
+   {
+      return;
+   }
+
+   const CoarseFineTransformations& cf_tr = mesh_lor->GetRefinementTransforms();
+
+   int nref_max = 0;
+   Array<Geometry::Type> geoms;
+   mesh_ho->GetGeometries(mesh_ho->Dimension(), geoms); 
+   for (int ig = 0; ig < geoms.Size(); ++ig)
+   {
+      Geometry::Type geom = geoms[ig];
+      nref_max = std::max(nref_max, cf_tr.point_matrices[geom].SizeK());
+   }
+
+   BuildHo2Lor(nel_ho, nel_lor, cf_tr);
+
+   // **************************
+   // inv lumped M_L 
+   // **************************
+   // TODO: need to add quadratureFunctionCoefficient into this M_L assembly 
+
+   // ML_inv_ea contains the inverse lumped (row sum) mass matrix. Note that the
+   // method will also work with a full (consistent) mass matrix, though this is
+   // not implemented here. L refers to the low-order refined mesh
+   ML_inv_ea.SetSize(ndof_lor);
+   ML_inv_ea = 0.0;
+
+   // Compute ML_inv_ea
+   for (int iho = 0; iho < nel_ho; ++iho)
+   {
+      Array<int> lor_els;
+      ho2lor.GetRow(iho, lor_els);
+      int nref = ho2lor.RowSize(iho);
+
+      Geometry::Type geom = mesh_ho->GetElementBaseGeometry(iho);
+      const FiniteElement& fe_lor = *fes_lor.GetFE(lor_els[0]);
+      int nedof_lor = fe_lor.GetDof();
+
+      // Instead of using a MassIntegrator, manually loop over integration
+      // points so we can row sum and store the diagonal as a Vector.
+      Vector ML_el(nedof_lor);
+      Vector shape_lor(nedof_lor);
+      Array<int> dofs_lor(nedof_lor);
+
+      for (int iref = 0; iref < nref; ++iref)
+      {
+         int ilor = lor_els[iref];
+         ElementTransformation* el_tr = fes_lor.GetElementTransformation(ilor);
+
+         int order = 2 * fe_lor.GetOrder() + el_tr->OrderW();
+         const IntegrationRule* ir = &IntRules.Get(geom, order);
+         ML_el = 0.0;
+         for (int i = 0; i < ir->GetNPoints(); ++i)
+         {
+            const IntegrationPoint& ip_lor = ir->IntPoint(i);
+            fe_lor.CalcShape(ip_lor, shape_lor);
+            el_tr->SetIntPoint(&ip_lor);
+            ML_el += (shape_lor *= (el_tr->Weight() * ip_lor.weight));
+         }
+         fes_lor.GetElementDofs(ilor, dofs_lor);
+         ML_inv_ea.AddElementVector(dofs_lor, ML_el);
+      }
+   }
+   // DOF by DOF inverse of non-zero entries
+   LumpedMassInverse(ML_inv_ea); 
+
+   // **************************
+   // ElementRestrictionOperator for HO
+   // **************************
+   elem_restrict_h = fes_ho.GetElementRestriction(ElementDofOrdering::NATIVE);
+
+   // **************************
+   // mixed mass M_LH from L2Space
+   // **************************
+
+   // dynamic_cast to check if QuadFuncCoeff; if yes, continue; if no, return error
+   auto qfunc_coeff = dynamic_cast<QuadratureFunctionCoefficient*>(coeff_);
+   if (qfunc_coeff == NULL)
+   {
+      mfem_error("Not a QuadratureFunctionCoefficient");
+   }
+   QuadratureFunction qfunc = qfunc_coeff->GetQuadFunction();
+   // Store the mixed mass matrix integration rule, which is assumed same on all elements
+   const IntegrationRule ir = qfunc.GetIntRule(0); 
+
+   offsets.SetSize(nel_ho+1);
+   offsets[0] = 0;
+   for (int iho = 0; iho < nel_ho; ++iho)
+   {
+      int nref = ho2lor.RowSize(iho);
+      const FiniteElement &fe_ho = *fes_ho.GetFE(iho);
+      const FiniteElement &fe_lor = *fes_lor.GetFE(ho2lor.GetRow(iho)[0]);
+      offsets[iho+1] = offsets[iho] + fe_ho.GetDof()*fe_lor.GetDof()*nref;
+   }
+
+   IntegrationPointTransformation ip_tr;
+   IsoparametricTransformation &emb_tr = ip_tr.Transf;
+
+   //Gather basis functions (B_L, B_HO)
+   //and data at quadrature points
+   DenseTensor B_L, B_H, D;
+   {
+      //Assume all HO elements are LOR in the same way
+      const int iho = 0;
+      {
+         Array<int> lor_els;
+         ho2lor.GetRow(iho, lor_els);
+         int nref = ho2lor.RowSize(iho);
+
+         Geometry::Type geom = mesh_ho->GetElementBaseGeometry(iho);
+         const FiniteElement &fe_ho = *fes_ho.GetFE(iho);
+         const FiniteElement &fe_lor = *fes_lor.GetFE(lor_els[0]);
+
+         //Allocate space for DenseTensors
+         ElementTransformation *el_tr = fes_lor.GetElementTransformation(0);
+         int order = fe_lor.GetOrder() + fe_ho.GetOrder() + el_tr->OrderW();
+         const IntegrationRule* ir = &IntRules.Get(geom, order);
+         int qPts = ir->GetNPoints();
+
+         //Containers for the basis functions sampled
+         //at quadrature points
+         B_L.SetSize(qPts, fe_lor.GetDof(), nref, d_mt_);
+         B_H.SetSize(qPts, fe_ho.GetDof(), nref, d_mt_);
+         D.SetSize(qPts, nref, nel_ho, d_mt_);
+
+         const DofToQuad *maps_lor = &fe_lor.GetDofToQuad(*ir, DofToQuad::TENSOR);
+
+         const GeometricFactors *geo_facts =
+            mesh_lor->GetGeometricFactors(*ir, GeometricFactors::DETERMINANTS);
+
+         const int Q1D         = maps_lor->nqpt;
+
+         const int dim = mesh_ho->Dimension();
+
+         MFEM_ASSERT(nel_ho*nref == nel_lor, "we expect nel_ho*nref == nel_lor");
+         MFEM_VERIFY(D.TotalSize() == qfunc.Size(), "Dimensions don't match  "<<D.TotalSize()<<" "<<qfunc.Size()); 
+
+         //*********************************
+         // Setup data at quadrature points
+         //*********************************
+         if (dim == 1)
+         {
+
+            const auto W = Reshape(ir->GetWeights().Read(), Q1D);  // grabbing the weights of the integration rule 
+            const auto J = Reshape(geo_facts->detJ.Read(), Q1D, nel_lor);  // 
+            const auto d_D = Reshape(D.Write(), qPts, nref, nel_ho);  // diagonal at the quadrature points
+            const auto d_qfunc = Reshape(qfunc.Read(), qPts, nref, nel_ho);
+
+            mfem::forall(nel_ho, [=] MFEM_HOST_DEVICE (int iho)
+            {
+               for (int iref = 0; iref < nref; ++iref)
+               {
+                  const int lo_el_id = iref + nref*iho;
+                  for (int qx=0; qx<Q1D; ++qx)
+                  {
+                     const real_t detJ = J(qx, lo_el_id);
+                     d_D(qx, iref, iho) = W(qx) * detJ * d_qfunc(qx, iref, iho); 
+                  }
+               }
+            });
+
+         }
+
+         if (dim == 2)  
+         {
+
+            const auto W = Reshape(ir->GetWeights().Read(), Q1D, Q1D);
+            const auto J = Reshape(geo_facts->detJ.Read(), Q1D,Q1D, nel_lor);
+            const auto d_D = Reshape(D.Write(), qPts, nref, nel_ho);
+            const auto d_qfunc = Reshape(qfunc.Read(), qPts, nref, nel_ho);
+
+            mfem::forall(nel_ho, [=] MFEM_HOST_DEVICE (int iho)
+            {
+               for (int iref = 0; iref < nref; ++iref)
+               {
+
+                  const int lo_el_id = iref + nref*iho;
+                  for (int qy=0; qy<Q1D; ++qy)
+                  {
+                     for (int qx=0; qx<Q1D; ++qx)
+                     {
+                        const int q = qx + Q1D*qy;
+                        const real_t detJ = J(qx, qy, lo_el_id);
+                        d_D(q, iref, iho) = W(qx, qy) * detJ * d_qfunc(q, iref, iho);
+                     }
+                  }
+
+               }
+            });
+         }
+
+         if (dim == 3)
+         {
+
+            const auto W = Reshape(ir->GetWeights().Read(), Q1D, Q1D, Q1D);
+            const auto J = Reshape(geo_facts->detJ.Read(), Q1D, Q1D, Q1D, nel_lor);
+            const auto d_D = Reshape(D.Write(), qPts, nref, nel_ho);
+            const auto d_qfunc = Reshape(qfunc.Read(), qPts, nref, nel_ho);
+
+            mfem::forall(nel_ho, [=] MFEM_HOST_DEVICE (int iho)
+            {
+               for (int iref = 0; iref < nref; ++iref)
+               {
+
+                  const int lo_el_id = iref + nref*iho;
+
+                  for (int qz=0; qz<Q1D; qz++)
+                  {
+                     for (int qy=0; qy<Q1D; ++qy)
+                     {
+                        for (int qx=0; qx<Q1D; ++qx)
+                        {
+                           const int q = qx + Q1D*qy + Q1D*Q1D*qz;
+                           const real_t detJ = J(qx, qy, qz, lo_el_id);
+                           d_D(q, iref, iho) = W(qx, qy, qz) * detJ * d_qfunc(q, iref, iho); 
+                        }
+                     }
+                  }
+
+               }
+            });
+         }
+
+         emb_tr.SetIdentityTransformation(geom);
+         const DenseTensor &pmats = cf_tr.point_matrices[geom];
+
+         //Collect the basis functions
+         for (int iref = 0; iref < nref; ++iref)
+         {
+
+            int ilor = lor_els[iref];
+
+            // Now assemble the block-row of the mixed mass matrix associated
+            // with integrating HO functions against LOR functions on the LOR
+            // sub-element.
+
+            // Create the transformation that embeds the fine low-order element
+            // within the coarse high-order element in reference space
+            emb_tr.SetPointMat(pmats(cf_tr.embeddings[ilor].matrix));
+
+            DenseMatrix &b_lo = B_L(ilor);
+            DenseMatrix &b_ho = B_H(ilor);
+
+            ElemMixedMass(geom, fe_ho, fe_lor, el_tr, ip_tr, b_lo, b_ho);
+
+         }//loop over subcells of ho element
+
+         //-------[End of quadrature point setup]-----
+      }//
+
+   } //Competed setup of basis function and quadrature point
+
+   //Assemble mixed mass matrix
+   {
+      int iho = 0;
+      Array<int> lor_els;
+      ho2lor.GetRow(iho, lor_els);
+      int nref = ho2lor.RowSize(iho);
+
+      const FiniteElement &fe_ho = *fes_ho.GetFE(iho);
+      const FiniteElement &fe_lor = *fes_lor.GetFE(lor_els[0]);
+      ndof_ho = fe_ho.GetDof();
+      ndof_lor = fe_lor.GetDof();
+
+      const int qPts = D.SizeI();
+
+      M_mixed_all_ea.SetSize(ndof_lor*ndof_ho*nref*nel_ho, d_mt_);
+      M_mixed_all_ea = 0.0;
+
+      //Rows x columns
+      //Recall MFEM is column major
+
+      //rows x columns is inverted  - matrix is ndof_lor x ndof_ho
+      auto v_M_mixed_all_ea = mfem::Reshape(M_mixed_all_ea.Write(), ndof_lor, ndof_ho, nref,
+                                         nel_ho);
+
+      const int fe_ho_ndof = fe_ho.GetDof();
+      const int fe_lor_ndof = fe_lor.GetDof();
+
+      auto d_B_L = mfem::Reshape(B_L.Read(), qPts, fe_lor_ndof, nref);
+      auto d_B_H = mfem::Reshape(B_H.Read(), qPts, fe_ho_ndof, nref);
+      auto d_D   = mfem::Reshape(D.Read(), qPts, nref, nel_ho);
+
+      mfem::forall(nel_ho, [=] MFEM_HOST_DEVICE (int iho)
+      {
+
+         for (int iref = 0; iref < nref; ++iref)
+         {
+            // (B_lo_dofs x Q) x (Q x B_ho_dofs)
+            for (int bh=0; bh<fe_ho_ndof; ++bh)
+            {
+               for (int bl=0; bl<fe_lor_ndof; ++bl)
+               {
+
+                  real_t dot = 0.0;
+                  for (int qi=0; qi<qPts; ++qi)
+                  {
+                     dot += d_B_L(qi, bl, iref) *  d_D(qi, iref, iho) * d_B_H(qi, bh, iref);
+                  }
+
+                  //column major storange
+                  v_M_mixed_all_ea(bl , bh, iref, iho) = dot;
+               }
+            }
+
+         }
+      });
+   } // end of mixed assembly mass matrix
+
+   // M_mixed_all_ea = PullL2SpaceDeviceM_LH(fes_ho, fes_lor);
+   
+   // printf("Successfully build M_LH \n"); 
+   // std::cout << "M_LH has size " << M_mixed_all_ea.Size() << " and is " << std::endl;
+   // int iho = 19;  // 20
+   // int nref = 16;  // 16
+   // int iref = 15;
+   // for (int i = 0; i < ndof_lor; ++i) {
+   //    for (int j = 0; j < ndof_ho; ++j) {
+   //       printf("%+14.6e", M_mixed_all_ea.HostRead()[i+j*ndof_lor+iho*nref*ndof_lor*ndof_ho+iref*ndof_lor*ndof_ho]);
+   //       if (j%4 == 3) {printf("\n");}
+   //    }
+   //    printf("\n");
+   // }
+   // for (int i = 0; i < 100; ++i){
+   //    printf("%9.5f", M_mixed_all_ea.HostRead()[i]);
+   //    if (i%5 == 4) {printf("\n");}
+   // }
+
+
+   // **************************
+   // ElementRestrictionOperator for LOR
+   // **************************
+   elem_restrict_l = fes_lor.GetElementRestriction(ElementDofOrdering::NATIVE);
+
 }
 
 void L2ProjectionGridTransfer::L2ProjectionH1Space::Mult(
    const Vector& x, Vector& y) const
 {
+   if (use_device)
+   {
+      DeviceMult(x,y);
+      if (!verify_solution) {return;}
+   }
+
    Vector X(fes_ho.GetTrueVSize());
    Vector X_dim(R->Width());
 
@@ -1410,11 +1809,125 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::Mult(
    }
 
    SetFromTDofs(fes_lor, Y, y);
+
+   if (verify_solution)
+   {
+      Vector y_temp(y.Size());
+      DeviceMult(x, y_temp);
+      y_temp -= y;
+      real_t error = y_temp.Norml2();
+      if (error > ho_lor_tol)
+      {
+         MFEM_VERIFY(false, "Mult difference too high = "<<error);
+      }
+   }
+}
+
+void L2ProjectionGridTransfer::L2ProjectionH1Space::DeviceMult(
+   const Vector& x, Vector& y) const
+{
+   const int vdim = fes_ho.GetVDim();
+   const int iho = 0;
+   const int nref = ho2lor.RowSize(iho);
+   const int ndof_ho = fes_ho.GetFE(iho)->GetDof();
+   const int ndof_lor = fes_lor.GetFE(ho2lor.GetRow(iho)[0])->GetDof();
+   const Mesh *mesh_ho = fes_ho.GetMesh();
+   const int nel_ho = mesh_ho->GetNE();
+
+   Vector tempx(elem_restrict_h->Height());
+   tempx = 0.0;
+   elem_restrict_h->Mult(x, tempx);
+
+   Vector tempy(ndof_lor*nref*vdim*nel_ho);
+
+   auto v_M_mixed_ea = mfem::Reshape(M_mixed_all_ea.Read(), ndof_lor, ndof_ho, nref, nel_ho); 
+   auto v_tempx    = mfem::Reshape(tempx.Read(), ndof_ho, vdim, nel_ho);
+   auto v_tempy    = mfem::Reshape(tempy.Write(), ndof_lor, nref, vdim, nel_ho);
+
+   mfem::forall(nel_ho, [=] MFEM_HOST_DEVICE (int iho)
+   {
+
+      for (int v=0; v<vdim; ++v)
+      {
+
+         for (int i=0; i<nref; ++i)
+         {
+            for (int j=0; j<ndof_lor; ++j)
+            {
+
+               real_t dot = 0.0;
+               for (int k=0; k<ndof_ho; ++k)
+               {
+                  dot += v_M_mixed_ea(j, k, i, iho) * v_tempx(k, v, iho);
+               }
+
+               v_tempy(j, i, v, iho) = dot;
+            }
+         }
+      }
+   });
+
+   Vector tempy2(y.Size());
+   tempy2 = 0.0;
+   elem_restrict_l->MultTranspose(tempy, tempy2);
+
+   auto v_ML_inv_ea  = mfem::Reshape(ML_inv_ea.Read(), ML_inv_ea.Size());
+   auto v_tempy2  = mfem::Reshape(tempy2.Read(), tempy2.Size());
+   auto v_y   = mfem::Reshape(y.Write(), y.Size());
+   MFEM_ASSERT(ML_inv_ea.Size() == tempy2.Size(), "sizes not the same");
+
+   mfem::forall(ML_inv_ea.Size(), [=] MFEM_HOST_DEVICE (int i)
+   {
+      v_y(i) = v_ML_inv_ea(i) * v_tempy2(i);
+   });
+}
+
+Vector L2ProjectionGridTransfer::L2ProjectionH1Space::PullL2SpaceDeviceM_LH(
+                        const FiniteElementSpace& coarse_fes_,
+                        const FiniteElementSpace& fine_fes_)
+{
+   Mesh* mesh_ho = coarse_fes_.GetMesh();
+   int dim = mesh_ho->Dimension();
+   Mesh* mesh_lor = fine_fes_.GetMesh();
+
+   int order = coarse_fes_.FEColl()->GetOrder();
+   int lorder = fine_fes_.FEColl()->GetOrder();
+
+   FiniteElementCollection *fec_l2, *fec_lor_l2;
+   fec_l2 = new L2_FECollection(order, dim);
+   fec_lor_l2 = new L2_FECollection(lorder, dim);
+
+   FiniteElementSpace fespace_l2(mesh_ho, fec_l2);
+   FiniteElementSpace fespace_lor_l2(mesh_lor, fec_lor_l2);
+
+   Geometry::Type geom = mesh_ho->GetElementBaseGeometry(0); 
+   const FiniteElement &fe = *fespace_l2.GetFE(0);
+   const FiniteElement &fe_lor = *fespace_lor_l2.GetFE(0); 
+   ElementTransformation *el_tr = fespace_lor_l2.GetElementTransformation(0);
+   int qorder = fe_lor.GetOrder() + fe.GetOrder() + el_tr->OrderW(); 
+   const IntegrationRule* ir = &IntRules.Get(geom, qorder);
+
+   QuadratureSpace qspace(*mesh_lor, *ir); 
+   QuadratureFunction qfunc(&qspace);
+   qfunc = 1.0; 
+   QuadratureFunctionCoefficient coeff_l2(qfunc); 
+
+   L2ProjectionL2Space *gt_l2;
+   gt_l2 = new L2ProjectionL2Space(fespace_l2, fespace_lor_l2, &coeff_l2, true, true);
+
+   Vector M_LH_l2 = gt_l2->M_mixed_all;
+
+   return M_LH_l2;
 }
 
 void L2ProjectionGridTransfer::L2ProjectionH1Space::MultTranspose(
    const Vector& x, Vector& y) const
 {
+   if (use_device)
+   {
+      DeviceMultTranspose(x,y); 
+      if (!verify_solution) {return;}
+   }
    Vector X(fes_lor.GetTrueVSize());
    Vector X_dim(R->Height());
 
@@ -1435,6 +1948,24 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::MultTranspose(
    }
 
    SetFromTDofsTranspose(fes_ho, Y, y);
+
+   if (verify_solution)
+   {
+      Vector y_temp(y.Size());
+      DeviceMultTranspose(x, y_temp);
+      y_temp -= y;
+      real_t error = y_temp.Norml2();
+      if (error > ho_lor_tol)
+      {
+         MFEM_VERIFY(false, "MultTranspose difference too high = "<<error);
+      }
+   }
+}
+
+void L2ProjectionGridTransfer::L2ProjectionH1Space::DeviceMultTranspose(
+   const Vector& x, Vector& y) const
+{
+   // TODO: Transpose of action of R
 }
 
 void L2ProjectionGridTransfer::L2ProjectionH1Space::Prolongate(
@@ -1508,7 +2039,7 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::SetAbsTol(real_t p_atol_)
 std::pair<
 std::unique_ptr<SparseMatrix>,
 std::unique_ptr<SparseMatrix>>
-                            L2ProjectionGridTransfer::L2ProjectionH1Space::ComputeSparseRAndM_LH()
+                            L2ProjectionGridTransfer::L2ProjectionH1Space::ComputeSparseRAndM_LH(bool GetM_LHError, bool GetML_invError)
 {
    std::pair<std::unique_ptr<SparseMatrix>,
        std::unique_ptr<SparseMatrix>> r_and_mlh;
@@ -1609,6 +2140,7 @@ std::unique_ptr<SparseMatrix>>
    IsoparametricTransformation& emb_tr = ip_tr.Transf;
 
    // Compute M_LH and R
+   real_t error = 0.0;
    for (int iho = 0; iho < nel_ho; ++iho)
    {
       Array<int> lor_els;
@@ -1627,6 +2159,8 @@ std::unique_ptr<SparseMatrix>>
       DenseMatrix M_LH_el(nedof_lor, nedof_ho);
       DenseMatrix R_el(nedof_lor, nedof_ho);
 
+      // auto M_mixed_all_ea_mat = Reshape(M_mixed_all_ea.Read(), nedof_lor, nedof_ho, nref, nel_ho);
+
       for (int iref = 0; iref < nref; ++iref)
       {
          int ilor = lor_els[iref];
@@ -1637,6 +2171,7 @@ std::unique_ptr<SparseMatrix>>
          emb_tr.SetPointMat(pmats(cf_tr.embeddings[ilor].matrix));
 
          ElemMixedMass(geom, fe_ho, fe_lor, el_tr, ip_tr, M_LH_el);
+         // if (iref > 0) {M_LH_el = 0.0;}
 
          Array<int> dofs_lor(nedof_lor);
          fes_lor.GetElementDofs(ilor, dofs_lor);
@@ -1650,6 +2185,20 @@ std::unique_ptr<SparseMatrix>>
          fes_ho.GetElementDofs(iho, dofs_ho);
          r_and_mlh.second->AddSubMatrix(dofs_lor, dofs_ho, M_LH_el);
          r_and_mlh.first->AddSubMatrix(dofs_lor, dofs_ho, R_el);
+
+         if (GetM_LHError == true) {
+            real_t mat_error = 0.0;
+            for (int i = 0; i < nedof_lor; ++i) {
+               for (int j = 0; j < nedof_ho; ++j) {
+                  mat_error += (M_mixed_all_ea.Read()[i+j*nedof_lor+offsets[iho]+nedof_lor*nedof_ho*iref] - M_LH_el(i,j))*(M_mixed_all_ea.Read()[i+j*nedof_lor+offsets[iho]+nedof_lor*nedof_ho*iref] - M_LH_el(i,j));
+               }
+            }
+            error += std::sqrt(mat_error);
+            if (error > ho_lor_tol)
+            {
+               MFEM_VERIFY(false, "error in M_LH difference too high = " << error);
+            }
+         }
       }
    }
 
@@ -1878,7 +2427,8 @@ void L2ProjectionGridTransfer::BuildF()
    {
       if (!Parallel())
       {
-         F = new L2ProjectionH1Space(dom_fes, ran_fes);
+         F = new L2ProjectionH1Space(dom_fes, ran_fes, coeff, 
+                                     use_device, verify_solution);
       }
       else
       {
@@ -1892,9 +2442,9 @@ void L2ProjectionGridTransfer::BuildF()
       }
    }
    else
-   {
-      F = new L2ProjectionL2Space(dom_fes, ran_fes, use_device,
-                                  verify_solution);
+   {  
+      F = new L2ProjectionL2Space(dom_fes, ran_fes, coeff, 
+                                  use_device, verify_solution);
    }
 }
 

fem/transfer.hpp

@@ -245,17 +245,21 @@ public:
 
       const bool use_device, verify_solution;
       MemoryType d_mt_;
+      Coefficient *coeff;
 
    public:
+      Vector M_mixed_all;
       L2ProjectionL2Space(const FiniteElementSpace& fes_ho_,
                           const FiniteElementSpace& fes_lor_,
+                          Coefficient* coeff_,
                           const bool use_device_,
                           const bool verify_solution_,
                           MemoryType d_mt_ = MemoryType::DEFAULT);
 
       /*Same as above but assembles and stores R_ea, P_ea */
       void DeviceL2ProjectionL2Space(const FiniteElementSpace& fes_ho_,
-                                     const FiniteElementSpace& fes_lor_);
+                                     const FiniteElementSpace& fes_lor_,
+                                     Coefficient* coeff_);
 
       /// Maps <tt>x</tt>, primal field coefficients defined on a coarse mesh
       /// with a higher order L2 finite element space, to <tt>y</tt>, primal
@@ -312,18 +316,40 @@ protected:
 
    friend class L2ProjectionL2Space;
 
+   //Class below must be public as we now have device code
+public:
+
    /** Projection operator between a H1 high-order finite element space on a
        coarse mesh, and a H1 low-order finite element space on a refined mesh
        (LOR). */
    class L2ProjectionH1Space : public L2Projection
    {
+      const bool use_device, verify_solution;
+      MemoryType d_mt_;
+      Coefficient* coeff;
+      Array<int> offsets;
+
+      const ElementRestrictionOperator* elem_restrict_h;
+      Vector M_mixed_all_ea;
+      const ElementRestrictionOperator* elem_restrict_l;
+      Vector ML_inv_ea;
+
+
    public:
       L2ProjectionH1Space(const FiniteElementSpace &fes_ho_,
-                          const FiniteElementSpace &fes_lor_);
-#ifdef MFEM_USE_MPI
-      L2ProjectionH1Space(const ParFiniteElementSpace &pfes_ho_,
-                          const ParFiniteElementSpace &pfes_lor_);
-#endif
+                          const FiniteElementSpace &fes_lor_,
+                          Coefficient *coeff_, 
+                          const bool use_device_,
+                          const bool verify_solution_, 
+                          MemoryType d_mt_ = MemoryType::DEFAULT);
+      #ifdef MFEM_USE_MPI
+            L2ProjectionH1Space(const ParFiniteElementSpace &pfes_ho_,
+                              const ParFiniteElementSpace &pfes_lor_);
+      #endif
+      /* Same as above but assembles action of R through ElementRestrictionOperator */
+      void DeviceL2ProjectionH1Space(const FiniteElementSpace &fes_ho_,
+                                     const FiniteElementSpace &fes_lor_, 
+                                     Coefficient* coeff_);
       /// Maps <tt>x</tt>, primal field coefficients defined on a coarse mesh
       /// with a higher order H1 finite element space, to <tt>y</tt>, primal
       /// field coefficients defined on a refined mesh with a low order H1
@@ -331,6 +357,10 @@ protected:
       /// the coarse mesh. Coefficients are computed through minimization of L2
       /// error between the fields.
       virtual void Mult(const Vector& x, Vector& y) const;
+
+      // Perform mult on the device (same as above)
+      void DeviceMult(const Vector&x, Vector& y) const;
+
       /// Maps <tt>x</tt>, dual field coefficients defined on a refined mesh
       /// with a low order H1 finite element space, to <tt>y</tt>, dual field
       /// coefficients defined on a coarse mesh with a higher order H1 finite
@@ -339,6 +369,9 @@ protected:
       /// error between the primal fields. Note, if the <tt>x</tt>-coefficients
       /// come from ProlongateTranspose, then mass is conserved.
       virtual void MultTranspose(const Vector& x, Vector& y) const;
+
+      void DeviceMultTranspose(const Vector& x, Vector& y) const;
+      
       /// Maps <tt>x</tt>, primal field coefficients defined on a refined mesh
       /// with a low order H1 finite element space, to <tt>y</tt>, primal field
       /// coefficients defined on a coarse mesh with a higher order H1 finite
@@ -357,12 +390,15 @@ protected:
       virtual void ProlongateTranspose(const Vector& x, Vector& y) const;
       virtual void SetRelTol(real_t p_rtol_);
       virtual void SetAbsTol(real_t p_atol_);
+
+      virtual Vector PullL2SpaceDeviceM_LH(const FiniteElementSpace& coarse_fes_,
+                        const FiniteElementSpace& fine_fes_);
    protected:
       /// Sets up the PCG solver (sets parameters, operator, and preconditioner)
       void SetupPCG();
       /// Computes on-rank R and M_LH matrices.
       std::pair<std::unique_ptr<SparseMatrix>,
-          std::unique_ptr<SparseMatrix>> ComputeSparseRAndM_LH();
+          std::unique_ptr<SparseMatrix>> ComputeSparseRAndM_LH(bool GetM_LHError, bool getML_invError);
       /// @brief Recovers vector of tdofs given a vector of dofs and a finite
       /// element space
       void GetTDofs(const FiniteElementSpace& fes, const Vector& x, Vector& X) const;
@@ -402,8 +438,11 @@ protected:
       // Used to compute P = (RT*M_LH)^(-1) M_LH^T
       std::unique_ptr<Operator> M_LH;
       std::unique_ptr<Operator> RTxM_LH;
+
+      friend class L2ProjectionL2Space;
    };
 
+protected:
    /** Mass-conservative prolongation operator going in the opposite direction
        as L2Projection. This operator is a left inverse to the L2Projection. */
    class L2Prolongation : public Operator
@@ -429,12 +468,16 @@ protected:
    bool force_l2_space;
 
 public:
+   // Coefficient for weighted integration in mass matrices; allows for spatial variation
+   Coefficient *coeff; 
+      
    L2ProjectionGridTransfer(FiniteElementSpace &coarse_fes_,
                             FiniteElementSpace &fine_fes_,
+                            Coefficient *coeff_,
                             bool force_l2_space_ = false,
                             MemoryType d_mt = MemoryType::DEFAULT)
       : GridTransfer(coarse_fes_, fine_fes_, d_mt),
-        F(NULL), B(NULL), force_l2_space(force_l2_space_)
+        F(NULL), B(NULL), force_l2_space(force_l2_space_), coeff(coeff_)
    { }
    virtual ~L2ProjectionGridTransfer();
 

linalg/batchlinalg.cpp

@@ -146,7 +146,6 @@ void BatchSolver::ComputeLU()
 
       MFEM_VERIFY(status == MFEM_CU_or_HIP(BLAS_STATUS_SUCCESS),
                   "Failed at blasDgetrfBatched");
-
    }
    else
 #endif
@@ -270,27 +269,61 @@ void ApplyBlkMult(const DenseTensor &Mat, const Vector &x,
                   Vector &y)
 {
    const int ndof = Mat.SizeI();
-   const int NE   = Mat.SizeK();
-   auto X         = Reshape(x.Read(), ndof, NE);
-   auto Y         = Reshape(y.Write(), ndof, NE);
-   auto Me        = Reshape(Mat.Read(), ndof, ndof, NE);
-
-   //Takes row major format
-   mfem::forall(ndof* NE, [=] MFEM_HOST_DEVICE (int tid)
+   MFEM_VERIFY(Mat.SizeI() == Mat.SizeJ(), "Batcched matrices are not square: not invertible");
+   const int NE = Mat.SizeK();
+   
+#if defined(MFEM_USE_CUDA_OR_HIP)
+   if (Device::Allows(Backend::DEVICE_MASK))
    {
-
-      const int c = tid % ndof;
-      const int e = tid / ndof;
-
-      {
-         double dot = 0;
-         for (int r = 0; r < ndof; ++r)
-         {
-            dot += Me(r, c, e) * X(r, e);
-         }
-         Y(c, e) = dot;
+      Array<double *> Mat_ptr(NE); 
+      Array<double *> x_ptr(NE);
+      Array<double *> y_ptr(NE);
+      for (int k = 0; k < NE; k++) {
+         Mat_ptr[k] = &const_cast<DenseTensor&>(Mat).ReadWrite()[ndof*ndof*k];
+         x_ptr[k] = &const_cast<Vector&>(x).ReadWrite()[ndof*k];
+         y_ptr[k] = &y.ReadWrite()[ndof*k];
       }
-   });
+
+      double alpha = 1.0;
+      double beta = 0.0; 
+      MFEM_cu_or_hip(blasStatus_t)
+      status = MFEM_cu_or_hip(blasDgemvBatched)(DeviceBlasHandle(), 
+                                                MFEM_CU_or_HIP(BLAS_OP_N), 
+                                                ndof, ndof,
+                                                &alpha, 
+                                                Mat_ptr.Read(), ndof, 
+                                                x_ptr.Read(), 1,
+                                                &beta, 
+                                                y_ptr.ReadWrite(), 1, 
+                                                NE);
+      MFEM_VERIFY(status == MFEM_CU_or_HIP(BLAS_STATUS_SUCCESS), 
+                  "Failed at blasDgemvBatched");
+   }
+   else
+#endif
+   {
+      // Hand written version
+      auto X         = Reshape(x.Read(), ndof, NE);
+      auto Y         = Reshape(y.Write(), ndof, NE);
+      auto Me        = Reshape(Mat.Read(), ndof, ndof, NE);
+
+      //Takes row major format
+      mfem::forall(ndof* NE, [=] MFEM_HOST_DEVICE (int tid)
+      {
+
+         const int c = tid % ndof;
+         const int e = tid / ndof;
+
+         {
+            double dot = 0;
+            for (int r = 0; r < ndof; ++r)
+            {
+               dot += Me(r, c, e) * X(r, e);
+            }
+            Y(c, e) = dot;
+         }
+      });
+   }
 }
 
 

miniapps/tools/lor-transfer.cpp

@@ -133,6 +133,23 @@ int main(int argc, char *argv[])
    FiniteElementSpace fespace(&mesh, fec);
    FiniteElementSpace fespace_lor(&mesh_lor, fec_lor);
 
+   // Build the integration rule that matches with quadrature on mixed mass matrix,  
+   // assuming HO elements are the same, and that all HO are LOR in the same way 
+   Geometry::Type geom = mesh.GetElementBaseGeometry(0); 
+   const FiniteElement &fe = *fespace.GetFE(0);
+   const FiniteElement &fe_lor = *fespace_lor.GetFE(0); 
+   ElementTransformation *el_tr = fespace_lor.GetElementTransformation(0);
+   int qorder = fe_lor.GetOrder() + fe.GetOrder() + el_tr->OrderW(); // 0 + 3 + 1
+   const IntegrationRule* ir = &IntRules.Get(geom, qorder);
+
+   QuadratureSpace qspace(mesh_lor, *ir); 
+   QuadratureFunction qfunc(&qspace);
+   qfunc = 1.0; 
+   // qfunc(2) = 7; // does not pass verify_solution
+   // qfunc(7) = 333.000001; // does not pass verify_solution
+   // qfunc(7) = 333.0000001; // passes verify_solution
+   QuadratureFunctionCoefficient coeff(qfunc); 
+
    GridFunction rho(&fespace);
    GridFunction rho_lor(&fespace_lor);
 
@@ -170,17 +187,18 @@ int main(int argc, char *argv[])
    }
    else
    {
-      gt = new L2ProjectionGridTransfer(fespace, fespace_lor);
+      gt = new L2ProjectionGridTransfer(fespace, fespace_lor, &coeff); 
    }
 
    gt->UseDevice(true);
    gt->VerifySolution(true);
 
    const Operator &R = gt->ForwardOperator();
+   // printf("Get past forward operator call \n");
 
    // HO->LOR restriction
    direction = "HO -> LOR @ LOR";
-   R.Mult(rho, rho_lor);
+   R.Mult(rho, rho_lor);  
    compute_mass(&fespace_lor, ho_mass, LOR_dc, "R(HO)    ");
    if (vis) { visualize(LOR_dc, "R(HO)", Wx, Wy); Wx += offx; }