Finished DeviceMult for H1Space

Co-authored-by: Tom Stitt <stitt4@llnl.gov>

config/defaults.mk

@@ -497,11 +497,11 @@ GSLIB_LIB = -L$(GSLIB_DIR)/lib -lgs
 
 # CUDA library configuration
 CUDA_OPT =
-CUDA_LIB = -lcusparse
+CUDA_LIB = -lcusparse -lcusolver -lcublas
 
 # HIP library configuration
 HIP_OPT =
-HIP_LIB = -L$(HIP_DIR)/lib $(XLINKER)-rpath,$(HIP_DIR)/lib -lhipsparse
+HIP_LIB = -L$(HIP_DIR)/lib $(XLINKER)-rpath,$(HIP_DIR)/lib -lhipsparse -lhipblas
 
 # OCCA library configuration
 OCCA_DIR = @MFEM_DIR@/../occa

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/ceed/interface/coefficient.hpp

@@ -295,7 +295,7 @@ void InitCoefficientWithIndices(mfem::Coefficient *Q, mfem::Mesh &mesh,
          const int e = d_indices[sub_e];
          out(q, sub_e) = in(q, e);
       });
-      m_indices.DeleteDevice();
+      m_indices.ReleaseDeviceMemory();
       InitVector(ceedCoeff->coeff, ceedCoeff->coeffVector);
       coeff_ptr = ceedCoeff;
    }
@@ -388,7 +388,7 @@ void InitCoefficientWithIndices(mfem::VectorCoefficient *VQ, mfem::Mesh &mesh,
             out(d, q, sub_e) = in(d, q, e);
          }
       });
-      m_indices.DeleteDevice();
+      m_indices.ReleaseDeviceMemory();
       InitVector(ceedCoeff->coeff, ceedCoeff->coeffVector);
       coeff_ptr = ceedCoeff;
    }

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/tmop.cpp

@@ -2855,13 +2855,13 @@ void TMOP_Integrator::ReleasePADeviceMemory(bool copy_to_host)
 {
    if (PA.enabled)
    {
-      PA.H.GetMemory().DeleteDevice(copy_to_host);
-      PA.H0.GetMemory().DeleteDevice(copy_to_host);
+      PA.H.GetMemory().ReleaseDeviceMemory(copy_to_host);
+      PA.H0.GetMemory().ReleaseDeviceMemory(copy_to_host);
       if (!copy_to_host && !PA.Jtr.GetMemory().HostIsValid())
       {
          PA.Jtr_needs_update = true;
       }
-      PA.Jtr.GetMemory().DeleteDevice(copy_to_host);
+      PA.Jtr.GetMemory().ReleaseDeviceMemory(copy_to_host);
    }
 }
 

fem/transfer.hpp

@@ -40,6 +40,11 @@ protected:
    OperatorHandle fw_t_oper; ///< Forward true-dof operator
    OperatorHandle bw_t_oper; ///< Backward true-dof operator
 
+   bool use_device;
+   bool verify_solution;
+
+   MemoryType d_mt_;
+
 #ifdef MFEM_USE_MPI
    bool parallel;
 #endif
@@ -60,11 +65,17 @@ protected:
 public:
    /** Construct a transfer algorithm between the domain, @a dom_fes_, and
        range, @a ran_fes_, FE spaces. */
-   GridTransfer(FiniteElementSpace &dom_fes_, FiniteElementSpace &ran_fes_);
+   GridTransfer(FiniteElementSpace &dom_fes_,
+                FiniteElementSpace &ran_fes_,
+                MemoryType d_mt_ = MemoryType::DEFAULT);
 
    /// Virtual destructor
    virtual ~GridTransfer() { }
 
+   void UseDevice(bool use_device_ = true) { use_device = use_device_;}
+
+   void VerifySolution(bool verify) { verify_solution = verify;}
+
    /** @brief Set the desired Operator::Type for the construction of all
        operators defined by the underlying transfer algorithm. */
    /** The default value is Operator::ANY_TYPE which typically corresponds to a
@@ -170,6 +181,7 @@ public:
 class L2ProjectionGridTransfer : public GridTransfer
 {
 protected:
+
    /** Abstract class representing projection operator between a high-order
        finite element space on a coarse mesh, and a low-order finite element
        space on a refined mesh (LOR). We assume that the low-order space,
@@ -206,8 +218,16 @@ protected:
                          const FiniteElement& fe_lor, ElementTransformation* el_tr,
                          IntegrationPointTransformation& ip_tr,
                          DenseMatrix& M_mixed_el) const;
+
+      void ElemMixedMass(Geometry::Type geom, const FiniteElement& fe_ho,
+                         const FiniteElement& fe_lor, ElementTransformation* el_tr,
+                         IntegrationPointTransformation& ip_tr,
+                         DenseMatrix& B_L, DenseMatrix& B_H) const;
    };
 
+   //Class below must be public as we now have device code
+public:
+
    /** Class for projection operator between a L2 high-order finite element
        space on a coarse mesh, and a L2 low-order finite element space on a
        refined mesh (LOR). */
@@ -219,11 +239,28 @@ protected:
       // arrays. The entries of the i'th high-order element are stored at the
       // index given by offsets[i].
       mutable Array<real_t> R, P;
+      mutable Array<real_t> R_ea, P_ea;
+
       Array<int> offsets;
 
+      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_);
+                          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_,
+                                     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
       /// field coefficients defined on a refined mesh with a low order L2
@@ -231,6 +268,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 L2 finite element space, to <tt>y</tt>, dual field
       /// coefficients defined on a coarse mesh with a higher order L2 finite
@@ -239,6 +280,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 L2 finite element space, to <tt>y</tt>, primal field
       /// coefficients defined on a coarse mesh with a higher order L2 finite
@@ -247,6 +291,9 @@ protected:
       /// left-inverse prolongation operation. This functionality is also
       /// provided as an Operator by L2Prolongation.
       virtual void Prolongate(const Vector& x, Vector& y) const;
+
+      void DeviceProlongate(const Vector& x, Vector& y) const;
+
       /// Maps <tt>x</tt>, dual field coefficients defined on a coarse mesh with
       /// a higher order L2 finite element space, to <tt>y</tt>, dual field
       /// coefficients defined on a refined mesh with a low order L2 finite
@@ -255,22 +302,54 @@ protected:
       /// conservative left-inverse prolongation operation. This functionality
       /// is also provided as an Operator by L2Prolongation.
       virtual void ProlongateTranspose(const Vector& x, Vector& y) const;
+
+      void DeviceProlongateTranspose(const Vector& x, Vector& y) const;
+
+
       virtual void SetRelTol(real_t p_rtol_) { } ///< No-op.
       virtual void SetAbsTol(real_t p_atol_) { } ///< No-op.
+
+      //friend class L2ProjectionGridTransfer;
    };
 
+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
@@ -278,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
@@ -286,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
@@ -304,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;
@@ -349,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
@@ -376,11 +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_,
-                            bool force_l2_space_ = false)
-      : GridTransfer(coarse_fes_, fine_fes_),
-        F(NULL), B(NULL), force_l2_space(force_l2_space_)
+                            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_), coeff(coeff_)
    { }
    virtual ~L2ProjectionGridTransfer();
 
@@ -391,6 +488,7 @@ public:
    virtual bool SupportsBackwardsOperator() const;
 private:
    void BuildF();
+
 };
 
 /// Matrix-free transfer operator between finite element spaces

general/backends.hpp

@@ -63,6 +63,28 @@
 #define MFEM_FOREACH_THREAD(i,k,N) for(int i=0; i<N; i++)
 #endif
 
+#if defined(MFEM_USE_CUDA)
+#include <cublas_v2.h>
+#include <cusolverDn.h>
+#define MFEM_cu_or_hip(stub) cu##stub
+#define MFEM_Cu_or_Hip(stub) Cu##stub
+#define MFEM_CU_or_HIP(stub) CU##stub
+#define MFEM_CUDA_or_HIP(stub) CUDA##stub
+#elif defined(MFEM_USE_HIP)
+#include <hipblas/hipblas.h>
+#include <hip/hip_runtime.h>
+#include <hip/hip_runtime_api.h>
+#define MFEM_cu_or_hip(stub) hip##stub
+#define MFEM_Cu_or_Hip(stub) Hip##stub
+#define MFEM_CU_or_HIP(stub) HIP##stub
+#define MFEM_CUDA_or_HIP(stub) HIP##stub
+#else
+#define MFEM_cu_or_hip(stub)
+#define MFEM_Cu_or_Hip(stub)
+#define MFEM_CU_or_HIP(stub)
+#define MFEM_CUDA_or_HIP(stub)
+#endif
+
 // 'double' atomicAdd implementation for previous versions of CUDA
 #if defined(MFEM_USE_CUDA) && defined(__CUDA_ARCH__) && __CUDA_ARCH__ < 600
 MFEM_DEVICE inline real_t atomicAdd(real_t *add, real_t val)

general/mem_manager.cpp

@@ -954,7 +954,7 @@ void MemoryManager::Delete_(void *h_ptr, MemoryType h_mt, unsigned flags)
    }
 }
 
-void MemoryManager::DeleteDevice_(void *h_ptr, unsigned & flags)
+void MemoryManager::ReleaseDeviceMemory_(void *h_ptr, unsigned & flags)
 {
    const bool owns_device = flags & Mem::OWNS_DEVICE;
    if (owns_device)

general/mem_manager.hpp

@@ -423,7 +423,7 @@ public:
    /** @brief Delete the device pointer, if owned. If @a copy_to_host is true
        and the data is valid only on device, move it to host before deleting.
        Invalidates the device memory. */
-   inline void DeleteDevice(bool copy_to_host = true);
+   inline void ReleaseDeviceMemory(bool copy_to_host = true);
 
    /// Array subscript operator for host memory.
    inline T &operator[](int idx);
@@ -674,7 +674,7 @@ private: // Static methods used by the Memory<T> class
                                                 unsigned flags);
 
    /// Free device memory identified by its host pointer
-   static void DeleteDevice_(void *h_ptr, unsigned & flags);
+   static void ReleaseDeviceMemory_(void *h_ptr, unsigned & flags);
 
    /// Check if the memory types given the memory class are valid
    static bool MemoryClassCheck_(MemoryClass mc, void *h_ptr,
@@ -1108,12 +1108,12 @@ inline void Memory<T>::Delete()
 }
 
 template <typename T>
-inline void Memory<T>::DeleteDevice(bool copy_to_host)
+inline void Memory<T>::ReleaseDeviceMemory(bool copy_to_host)
 {
    if (flags & Registered)
    {
       if (copy_to_host) { Read(MemoryClass::HOST, capacity); }
-      MemoryManager::DeleteDevice_((void*)h_ptr, flags);
+      MemoryManager::ReleaseDeviceMemory_((void*)h_ptr, flags);
    }
 }
 

linalg/CMakeLists.txt

@@ -11,6 +11,7 @@
 
 list(APPEND SRCS
   auxiliary.cpp
+  batchlinalg.cpp
   blockmatrix.cpp
   blockoperator.cpp
   blockvector.cpp
@@ -31,6 +32,7 @@ list(APPEND SRCS
 
 list(APPEND HDRS
   auxiliary.hpp
+  batchlinalg.hpp
   blockmatrix.hpp
   blockoperator.hpp
   blockvector.hpp

linalg/batchlinalg.cpp

@@ -0,0 +1,394 @@
+// Copyright (c) 2010-2024, 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.
+
+
+// Implementation of batchlinalg class
+
+#include "batchlinalg.hpp"
+#include "../general/forall.hpp"
+#include "../general/backends.hpp"
+
+namespace mfem
+{
+
+#if defined(MFEM_USE_CUDA_OR_HIP)
+static MFEM_cu_or_hip(blasHandle_t) device_blas_handle = nullptr;
+
+const MFEM_cu_or_hip(blasHandle_t) & DeviceBlasHandle()
+{
+   if (!device_blas_handle)
+   {
+      auto status = MFEM_cu_or_hip(blasCreate)(&device_blas_handle);
+      MFEM_VERIFY(status == MFEM_CU_or_HIP(BLAS_STATUS_SUCCESS),
+                  "Cannot initialize GPU BLAS");
+      atexit([]()
+      {
+         MFEM_cu_or_hip(blasDestroy)(device_blas_handle);
+         device_blas_handle = nullptr;
+      });
+   }
+   return device_blas_handle;
+}
+#endif
+
+
+BatchSolver::BatchSolver(const DenseTensor &MatrixBatch,
+                         const SolveMode mode, MemoryType d_mt)
+   : mode_(mode)
+     // TODO: should this really be a copy?
+   , LUMatrixBatch_(MatrixBatch)
+   , d_mt_(d_mt)
+{
+   if (!setup_)
+   {
+      Setup();
+   }
+}
+
+BatchSolver::BatchSolver(const SolveMode mode, MemoryType d_mt) : mode_(mode),
+   d_mt_(d_mt) {}
+
+void BatchSolver::AssignMatrices(const DenseTensor &MatrixBatch)
+{
+   // TODO: should this really be a copy?
+   LUMatrixBatch_ = MatrixBatch;
+   setup_         = false;  //Setup is now false as the matrices have changed
+   lu_valid_      = false;
+
+   //Need to always call setup since the matrices have changed
+   if (!setup_)
+   {
+      Setup();
+   }
+}
+
+void BatchSolver::AssignMatrices(const Vector &vMatrixBatch,
+                                 const int size,
+                                 const int num_matrices)
+{
+   const int totalSize = size * size * num_matrices;
+   LUMatrixBatch_.SetSize(size, size,
+                          num_matrices, d_mt_);
+   double *d_LUMatrixBatch      = LUMatrixBatch_.Write();
+   const double *d_vMatrixBatch = vMatrixBatch.Read();
+
+   mfem::forall(totalSize, [=] MFEM_HOST_DEVICE (int i) { d_LUMatrixBatch[i] = d_vMatrixBatch[i]; });
+
+   AssignMatrices(LUMatrixBatch_);
+}
+
+void BatchSolver::GetInverse(DenseTensor &InvMatBatch) const
+{
+
+   /*
+    if (mode_ == SolveMode::INVERSE)
+    {
+       MFEM_WARNING("GetInverse with SolveMode::Inverse involves and extra memory copy, consider "
+                    "GetInverse(M, M_inv) instead");
+    }
+   */
+
+   if (!setup_)
+   {
+      mfem_error("BatchSolver has not been setup");
+   }
+
+   // use existing inverse
+   if (mode_ == SolveMode::INVERSE)
+   {
+      MFEM_VERIFY(InvMatrixBatch_.TotalSize() == InvMatBatch.TotalSize(),
+                  "Internal error, InvMatrixBatch_.TotalSize() != InvMatBatch.TotalSize()");
+
+      const double *d_M_inv = InvMatrixBatch_.Read();
+      double *d_out         = InvMatBatch.Write();
+
+      mfem::forall(InvMatrixBatch_.TotalSize(), [=] MFEM_HOST_DEVICE (int i) { d_out[i] = d_M_inv[i]; });
+   }
+   else if (mode_ == SolveMode::LU)
+   {
+      return ComputeInverse(InvMatBatch);
+   }
+   else
+   {
+      mfem_error("unsupported mode");
+   }
+}
+
+void BatchSolver::ComputeLU()
+{
+   if (lu_valid_)
+   {
+      return;
+   }
+
+#if defined(MFEM_USE_CUDA_OR_HIP)
+   if (Device::Allows(Backend::DEVICE_MASK))
+   {
+
+      Array<int> info_array(num_matrices_); // need to move to temp mem
+
+      MFEM_cu_or_hip(blasStatus_t)
+      status = MFEM_cu_or_hip(blasDgetrfBatched)(DeviceBlasHandle(),
+                                                 matrix_size_,
+                                                 lu_ptr_array_.ReadWrite(),
+                                                 matrix_size_,
+                                                 P_.Write(),
+                                                 info_array.Write(),
+                                                 num_matrices_);
+
+      MFEM_VERIFY(status == MFEM_CU_or_HIP(BLAS_STATUS_SUCCESS),
+                  "Failed at blasDgetrfBatched");
+   }
+   else
+#endif
+   {
+      //Hand written version
+      BatchLUFactor(LUMatrixBatch_, P_);
+   }
+   lu_valid_ = true;
+}
+
+
+void BatchSolver::ComputeInverse(DenseTensor &InvMatBatch) const
+{
+
+   MFEM_VERIFY(lu_valid_, "LU must be valid");
+
+#if defined(MFEM_USE_CUDA_OR_HIP)
+   if (Device::Allows(Backend::DEVICE_MASK))
+   {
+
+      Array<double *> inv_ptr_array(num_matrices_, d_mt_);
+
+      double *inv_ptr_base = InvMatBatch.Write();
+      double **d_inv_ptr_array = inv_ptr_array.Write();
+      const int matrix_size    = matrix_size_;
+      mfem::forall(num_matrices_, [=] MFEM_HOST_DEVICE (int i)
+      {
+         d_inv_ptr_array[i] = inv_ptr_base + i * matrix_size * matrix_size;
+      });
+
+      Array<int> info_array(num_matrices_, d_mt_);
+
+      //Invert matrices
+      MFEM_cu_or_hip(blasStatus_t) status =
+         MFEM_cu_or_hip(blasDgetriBatched)(DeviceBlasHandle(),
+                                           matrix_size_,
+                                           lu_ptr_array_.Read(),
+                                           matrix_size_,
+                                           // from hipblas.h: @param[in] ipiv
+                                           // we can const_cast safely because it's an "in" variable
+                                           const_cast<int *>(P_.Read()),
+                                           inv_ptr_array.ReadWrite(),
+                                           matrix_size_,
+                                           info_array.Write(),
+                                           num_matrices_);
+
+      MFEM_VERIFY(status == MFEM_CU_or_HIP(BLAS_STATUS_SUCCESS),
+                  "Failed at blasDgetriBatched");
+   }
+   else
+#endif
+   {
+      BatchInverseMatrix(LUMatrixBatch_, P_, InvMatBatch);
+   }
+}
+
+void BatchSolver::SolveLU(const Vector &b, Vector &x) const
+{
+
+   x = b;
+#if defined(MFEM_USE_CUDA_OR_HIP)
+   if (Device::Allows(Backend::DEVICE_MASK))
+   {
+      Array<double *> vector_array(num_matrices_, d_mt_);
+
+      double *x_ptr_base = x.ReadWrite();
+
+      double alpha = 1.0;
+
+      double **d_vector_array = vector_array.Write();
+      const int matrix_size   = matrix_size_;
+      mfem::forall(num_matrices_, [=] MFEM_HOST_DEVICE (int i) { d_vector_array[i] = x_ptr_base + i * matrix_size; });
+
+
+      MFEM_cu_or_hip(blasStatus_t)
+      status_lo = MFEM_cu_or_hip(blasDtrsmBatched)(DeviceBlasHandle(),
+                                                   MFEM_CU_or_HIP(BLAS_SIDE_LEFT),
+                                                   MFEM_CU_or_HIP(BLAS_FILL_MODE_LOWER),
+                                                   MFEM_CU_or_HIP(BLAS_OP_N),
+                                                   MFEM_CU_or_HIP(BLAS_DIAG_UNIT),
+                                                   matrix_size_,
+                                                   1,
+                                                   &alpha,
+                                                   lu_ptr_array_.Read(),
+                                                   matrix_size_,
+                                                   vector_array.ReadWrite(),
+                                                   matrix_size_,
+                                                   num_matrices_);
+
+      MFEM_VERIFY(status_lo == MFEM_CU_or_HIP(BLAS_STATUS_SUCCESS),
+                  "Failed at blasDtrsmBatched lo");
+
+      MFEM_cu_or_hip(blasStatus_t)
+      status_upp = MFEM_cu_or_hip(blasDtrsmBatched)(DeviceBlasHandle(),
+                                                    MFEM_CU_or_HIP(BLAS_SIDE_LEFT),
+                                                    MFEM_CU_or_HIP(BLAS_FILL_MODE_UPPER),
+                                                    MFEM_CU_or_HIP(BLAS_OP_N),
+                                                    MFEM_CU_or_HIP(BLAS_DIAG_NON_UNIT),
+                                                    matrix_size_,
+                                                    1,
+                                                    &alpha,
+                                                    lu_ptr_array_.Read(),
+                                                    matrix_size_,
+                                                    vector_array.ReadWrite(),
+                                                    matrix_size_,
+                                                    num_matrices_);
+
+      MFEM_VERIFY(status_upp == MFEM_CU_or_HIP(BLAS_STATUS_SUCCESS),
+                  "Failed at blasDtrsmBatched upper");
+   }
+   else
+#endif
+   {
+      BatchLUSolve(LUMatrixBatch_, P_, x);
+   }
+
+}
+
+//Hand rolled -- TODO replace with vendor call
+void ApplyBlkMult(const DenseTensor &Mat, const Vector &x,
+                  Vector &y)
+{
+   const int ndof = Mat.SizeI();
+   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))
+   {
+      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;
+         }
+      });
+   }
+}
+
+
+void BatchSolver::ApplyInverse(const Vector &b, Vector &x) const
+{
+   //Extend with vendor library capabilities
+   ApplyBlkMult(InvMatrixBatch_, b, x);
+}
+
+
+void BatchSolver::Setup()
+{
+   matrix_size_  = LUMatrixBatch_.SizeI();
+   num_matrices_ = LUMatrixBatch_.SizeK();
+
+   P_.SetSize(matrix_size_ * num_matrices_, d_mt_);
+   lu_ptr_array_.SetSize(num_matrices_, d_mt_);
+
+   // TODO: can this just be a Write?
+   double *lu_ptr_base = LUMatrixBatch_.ReadWrite();
+
+   const int matrix_size   = matrix_size_;
+   double **d_lu_ptr_array = lu_ptr_array_.Write();
+
+   mfem::forall(num_matrices_, [=] MFEM_HOST_DEVICE (int i)
+   {
+      d_lu_ptr_array[i] = lu_ptr_base + i * matrix_size * matrix_size;
+   });
+
+   switch (mode_)
+   {
+      case SolveMode::LU: ComputeLU(); break;
+
+      case SolveMode::INVERSE:
+         ComputeLU();
+         InvMatrixBatch_.SetSize(matrix_size_,
+                                 matrix_size_,
+                                 num_matrices_, d_mt_);
+         ComputeInverse(InvMatrixBatch_);
+         break;
+
+      default: mfem_error("Case not supported");
+   }
+   setup_ = true;
+}
+
+void BatchSolver::Mult(const Vector &b, Vector &x) const
+{
+   switch (mode_)
+   {
+      case SolveMode::LU: return SolveLU(b, x);
+
+      case SolveMode::INVERSE: return ApplyInverse(b, x);
+
+      default: mfem_error("Case not supported");
+   }
+}
+
+void BatchSolver::ReleaseMemory()
+{
+   LUMatrixBatch_.GetMemory().ReleaseDeviceMemory(false);
+   InvMatrixBatch_.GetMemory().ReleaseDeviceMemory(false);
+   P_.GetMemory().ReleaseDeviceMemory(false);
+   lu_ptr_array_.GetMemory().ReleaseDeviceMemory(false);
+}
+
+
+} // namespace mfem

linalg/batchlinalg.hpp

@@ -0,0 +1,85 @@
+// Copyright (c) 2010-2024, 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_BATCHLINALG
+#define MFEM_BATCHLINALG
+
+#include "../config/config.hpp"
+#include "../general/globals.hpp"
+#include "matrix.hpp"
+#include "densemat.hpp"
+
+namespace mfem
+{
+
+class BatchSolver
+{
+public:
+   enum struct SolveMode : int
+   {
+      LU,
+      INVERSE
+   };
+
+   // Compute inverse of a batch of matrices in M; layout must be (height, width, num_mats)
+   static void GetInverse(const DenseTensor &M, DenseTensor &M_inv);
+
+private:
+   BatchSolver::SolveMode mode_;
+
+   bool setup_    = false;
+   bool lu_valid_ = false;
+
+   DenseTensor LUMatrixBatch_, InvMatrixBatch_;
+   Array<int> P_;
+   Array<double *> lu_ptr_array_;
+   MemoryType d_mt_;
+
+   int num_matrices_, matrix_size_;
+
+   void ApplyInverse(const Vector &b, Vector &x) const;
+
+   // for nvcc
+public:
+   //Compute LU or Inverse
+   void Setup();
+
+   void ComputeLU();
+
+   void SolveLU(const Vector &b, Vector &x) const;
+
+   void ComputeInverse(DenseTensor &InvMatBatch) const;
+
+public:
+   BatchSolver() = delete;
+
+   BatchSolver(const SolveMode mode, MemoryType d_mt = MemoryType::DEFAULT);
+
+   BatchSolver(const DenseTensor &MatrixBatch, const SolveMode mode,
+               MemoryType d_mt = MemoryType::DEFAULT);
+
+   void AssignMatrices(const DenseTensor &MatrixBatch);
+
+   void AssignMatrices(const Vector &vMatrixBatch, const int size,
+                       const int num_matrices);
+
+   void GetInverse(DenseTensor &InvMatBatch) const;
+
+   //Solve linear system Ax = b
+   void Mult(const Vector &b, Vector &x) const;
+
+   void ReleaseMemory();
+
+};
+
+} // namespace mfem
+
+#endif

linalg/densemat.cpp

@@ -4633,16 +4633,40 @@ DenseTensor &DenseTensor::operator=(const DenseTensor &other)
 void BatchLUFactor(DenseTensor &Mlu, Array<int> &P, const real_t TOL)
 {
    const int m = Mlu.SizeI();
-   const int NE = Mlu.SizeK();
-   P.SetSize(m*NE);
+   const int len = Mlu.SizeK();
 
-   auto data_all = mfem::Reshape(Mlu.ReadWrite(), m, m, NE);
-   auto ipiv_all = mfem::Reshape(P.Write(), m, NE);
+   Vector vec(Mlu.ReadWrite(), m*m*len);
+   BatchLUFactor(m, len, vec, P, TOL);
+}
+
+void BatchLUSolve(const DenseTensor &Mlu, const Array<int> &P, Vector &X)
+{
+
+   const int m = Mlu.SizeI();
+   const int len = Mlu.SizeK();
+
+   auto data_all = mfem::Reshape(Mlu.Read(), m, m, len);
+   auto piv_all = mfem::Reshape(P.Read(), m, len);
+   auto x_all = mfem::Reshape(X.ReadWrite(), m, len);
+
+   mfem::forall(len, [=] MFEM_HOST_DEVICE (int e)
+   {
+      kernels::LUSolve(&data_all(0, 0,e), m, &piv_all(0, e), &x_all(0,e));
+   });
+
+}
+
+void BatchLUFactor(const int m, const int len,
+                   mfem::Vector &A, mfem::Array<int> &P, const double TOL)
+{
+   P.SetSize(m * len);
+   auto data_all = mfem::Reshape(A.ReadWrite(), m, m, len);
+   auto ipiv_all = mfem::Reshape(P.Write(), m, len);
    Array<bool> pivot_flag(1);
    pivot_flag[0] = true;
    bool *d_pivot_flag = pivot_flag.ReadWrite();
 
-   mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
+   mfem::forall(len, [=] MFEM_HOST_DEVICE (int e)
    {
       for (int i = 0; i < m; i++)
       {
@@ -4695,23 +4719,214 @@ void BatchLUFactor(DenseTensor &Mlu, Array<int> &P, const real_t TOL)
    });
 
    MFEM_ASSERT(pivot_flag.HostRead()[0], "Batch LU factorization failed \n");
+
 }
 
-void BatchLUSolve(const DenseTensor &Mlu, const Array<int> &P, Vector &X)
+
+void BatchLUSolve(mfem::Vector &Minv, int m, int len, mfem::Array<int> &P,
+                  mfem::Vector &X)
 {
+   auto data_all = mfem::Reshape(Minv.Read(), m, m, len);
+   auto piv_all  = mfem::Reshape(P.Read(), m, len);
+   auto x_all    = mfem::Reshape(X.ReadWrite(), m, len);
 
-   const int m = Mlu.SizeI();
-   const int NE = Mlu.SizeK();
 
-   auto data_all = mfem::Reshape(Mlu.Read(), m, m, NE);
-   auto piv_all = mfem::Reshape(P.Read(), m, NE);
-   auto x_all = mfem::Reshape(X.ReadWrite(), m, NE);
-
-   mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
+   mfem::forall(len, [=] MFEM_HOST_DEVICE (int e)
    {
-      kernels::LUSolve(&data_all(0, 0,e), m, &piv_all(0, e), &x_all(0,e));
-   });
 
+      const int *ipiv = &piv_all(0, e);
+      double *x       = &x_all(0, e);
+
+      // X <- P X
+      for (int i = 0; i < m; i++)
+      {
+         mfem::kernels::internal::Swap<double>(x[i], x[ipiv[i]]);
+      }
+
+      // X <- L^{-1} X
+      for (int j = 0; j < m; j++)
+      {
+         const double x_j = x[j];
+         for (int i = j + 1; i < m; i++)
+         {
+            x[i] -= data_all(i, j, e) * x_j;
+         }
+      }
+
+      // X <- U^{-1} X
+      for (int j = m - 1; j >= 0; j--)
+      {
+         const double x_j = (x[j] /= data_all(j, j, e));
+         for (int i = 0; i < j; i++)
+         {
+            x[i] -= data_all(i, j, e) * x_j;
+         }
+      }
+   });
+}
+
+void BatchInverseMatrix(const mfem::Vector &LU,
+                        const int m,
+                        const int len,
+                        const mfem::Array<int> &P,
+                        mfem::Vector &INV)
+{
+   auto data_all = mfem::Reshape(LU.Read(), m, m, len);
+   auto piv_all  = mfem::Reshape(P.Read(), m, len);
+   auto inv_all  = mfem::Reshape(INV.ReadWrite(), m, m, len);
+
+
+   mfem::forall(len, [=] MFEM_HOST_DEVICE (int e)
+   {
+      // A^{-1} = U^{-1} L^{-1} P
+      // X <- U^{-1} (set only the upper triangular part of X)
+      double *X          = &inv_all(0, 0, e);
+      double *x          = X;
+      const double *data = &data_all(0, 0, e);
+      const int *ipiv    = &piv_all(0, e);
+
+      for (int k = 0; k < m; k++)
+      {
+         const double minus_x_k = -(x[k] = 1.0 / data[k + k * m]);
+         for (int i = 0; i < k; i++)
+         {
+            x[i] = data[i + k * m] * minus_x_k;
+         }
+         for (int j = k - 1; j >= 0; j--)
+         {
+            const double x_j = (x[j] /= data[j + j * m]);
+            for (int i = 0; i < j; i++)
+            {
+               x[i] -= data[i + j * m] * x_j;
+            }
+         }
+         x += m;
+      }
+
+      // X <- X L^{-1} (use input only from the upper triangular part of X)
+      {
+         int k = m - 1;
+         for (int j = 0; j < k; j++)
+         {
+            const double minus_L_kj = -data[k + j * m];
+            for (int i = 0; i <= j; i++)
+            {
+               X[i + j * m] += X[i + k * m] * minus_L_kj;
+            }
+            for (int i = j + 1; i < m; i++)
+            {
+               X[i + j * m] = X[i + k * m] * minus_L_kj;
+            }
+         }
+      }
+      for (int k = m - 2; k >= 0; k--)
+      {
+         for (int j = 0; j < k; j++)
+         {
+            const double L_kj = data[k + j * m];
+            for (int i = 0; i < m; i++)
+            {
+               X[i + j * m] -= X[i + k * m] * L_kj;
+            }
+         }
+      }
+
+      // X <- X P
+      for (int k = m - 1; k >= 0; k--)
+      {
+         const int piv_k = ipiv[k];
+         if (k != piv_k)
+         {
+            for (int i = 0; i < m; i++)
+            {
+               //Swap<double>(X[i+k*m], X[i+piv_k*m]);
+               mfem::kernels::internal::Swap<double>(X[i + k * m], X[i + piv_k * m]);
+            }
+         }
+      }
+   });
+}
+
+void BatchInverseMatrix(const DenseTensor &LU,
+                        const Array<int> &P,
+                        DenseTensor &INV)
+{
+   const int m = LU.SizeI();
+   const int len = LU.SizeK();
+   auto data_all = mfem::Reshape(LU.Read(), m, m, len);
+   auto piv_all  = mfem::Reshape(P.Read(), m, len);
+   auto inv_all  = mfem::Reshape(INV.ReadWrite(), m, m, len);
+
+
+   mfem::forall(len, [=] MFEM_HOST_DEVICE (int e)
+   {
+      // A^{-1} = U^{-1} L^{-1} P
+      // X <- U^{-1} (set only the upper triangular part of X)
+      double *X          = &inv_all(0, 0, e);
+      double *x          = X;
+      const double *data = &data_all(0, 0, e);
+      const int *ipiv    = &piv_all(0, e);
+
+      for (int k = 0; k < m; k++)
+      {
+         const double minus_x_k = -(x[k] = 1.0 / data[k + k * m]);
+         for (int i = 0; i < k; i++)
+         {
+            x[i] = data[i + k * m] * minus_x_k;
+         }
+         for (int j = k - 1; j >= 0; j--)
+         {
+            const double x_j = (x[j] /= data[j + j * m]);
+            for (int i = 0; i < j; i++)
+            {
+               x[i] -= data[i + j * m] * x_j;
+            }
+         }
+         x += m;
+      }
+
+      // X <- X L^{-1} (use input only from the upper triangular part of X)
+      {
+         int k = m - 1;
+         for (int j = 0; j < k; j++)
+         {
+            const double minus_L_kj = -data[k + j * m];
+            for (int i = 0; i <= j; i++)
+            {
+               X[i + j * m] += X[i + k * m] * minus_L_kj;
+            }
+            for (int i = j + 1; i < m; i++)
+            {
+               X[i + j * m] = X[i + k * m] * minus_L_kj;
+            }
+         }
+      }
+      for (int k = m - 2; k >= 0; k--)
+      {
+         for (int j = 0; j < k; j++)
+         {
+            const double L_kj = data[k + j * m];
+            for (int i = 0; i < m; i++)
+            {
+               X[i + j * m] -= X[i + k * m] * L_kj;
+            }
+         }
+      }
+
+      // X <- X P
+      for (int k = m - 1; k >= 0; k--)
+      {
+         const int piv_k = ipiv[k];
+         if (k != piv_k)
+         {
+            for (int i = 0; i < m; i++)
+            {
+               //Swap<double>(X[i+k*m], X[i+piv_k*m]);
+               mfem::kernels::internal::Swap<double>(X[i + k * m], X[i + piv_k * m]);
+            }
+         }
+      }
+   });
 }
 
 } // namespace mfem

linalg/densemat.hpp

@@ -1290,6 +1290,22 @@ void BatchLUFactor(DenseTensor &Mlu, Array<int> &P, const real_t TOL = 0.0);
 void BatchLUSolve(const DenseTensor &Mlu, const Array<int> &P, Vector &X);
 
 
+void BatchLUFactor(const int m, const int len,
+                   mfem::Vector &A, mfem::Array<int> &P, const double TOL = 0.0);
+
+void BatchLUSolve(mfem::Vector &Minv, int m, int len, mfem::Array<int> &P,
+                  mfem::Vector &X);
+
+void BatchInverseMatrix(const Vector &LU,
+                        const int m,
+                        const int len,
+                        const Array<int> &P,
+                        Vector &INV);
+
+void BatchInverseMatrix(const DenseTensor &LU,
+                        const Array<int> &P,
+                        DenseTensor &INV);
+
 // Inline methods
 
 inline real_t &DenseMatrix::operator()(int i, int j)

linalg/linalg.hpp

@@ -32,6 +32,7 @@
 #include "invariants.hpp"
 #include "constraints.hpp"
 #include "auxiliary.hpp"
+#include "batchlinalg.hpp"
 
 #ifdef MFEM_USE_AMGX
 #include "amgxsolver.hpp"

linalg/sparsemat.cpp

@@ -25,10 +25,6 @@
 #include <cstring>
 
 #if defined(MFEM_USE_CUDA)
-#define MFEM_cu_or_hip(stub) cu##stub
-#define MFEM_Cu_or_Hip(stub) Cu##stub
-#define MFEM_CU_or_HIP(stub) CU##stub
-#define MFEM_CUDA_or_HIP(stub) CUDA##stub
 
 #if CUSPARSE_VERSION >=  11400
 #define MFEM_GPUSPARSE_ALG CUSPARSE_SPMV_CSR_ALG1
@@ -37,10 +33,6 @@
 #endif // CUSPARSE_VERSION >= 11400
 
 #elif defined(MFEM_USE_HIP)
-#define MFEM_cu_or_hip(stub) hip##stub
-#define MFEM_Cu_or_Hip(stub) Hip##stub
-#define MFEM_CU_or_HIP(stub) HIP##stub
-#define MFEM_CUDA_or_HIP(stub) HIP##stub
 
 // https://hipsparse.readthedocs.io/en/latest/usermanual.html#hipsparsespmvalg-t
 #define MFEM_GPUSPARSE_ALG HIPSPARSE_CSRMV_ALG1

linalg/vector.hpp

@@ -211,8 +211,8 @@ public:
    /** @brief Delete the device pointer, if owned. If @a copy_to_host is true
        and the data is valid only on device, move it to host before deleting.
        Invalidates the device memory. */
-   void DeleteDevice(bool copy_to_host = true)
-   { data.DeleteDevice(copy_to_host); }
+   void ReleaseDeviceMemory(bool copy_to_host = true)
+   { data.ReleaseDeviceMemory(copy_to_host); }
 
    /// Returns the size of the vector.
    inline int Size() const { return size; }

miniapps/tools/lor-transfer.cpp

@@ -33,6 +33,8 @@
 //
 // Sample runs:  lor-transfer
 //               lor-transfer -h1
+//               lor-transfer -d 'cuda'
+//               lor-transfer -d 'hip'
 //               lor-transfer -t
 //               lor-transfer -m ../../data/star-q2.mesh -lref 5 -p 4
 //               lor-transfer -m ../../data/star-mixed.mesh -lref 3 -p 2
@@ -75,6 +77,7 @@ int main(int argc, char *argv[])
    bool vis = true;
    bool useH1 = false;
    bool use_pointwise_transfer = false;
+   const char *device_config = "cpu";
 
    OptionsParser args(argc, argv);
    args.AddOption(&mesh_file, "-m", "--mesh",
@@ -95,6 +98,8 @@ int main(int argc, char *argv[])
    args.AddOption(&use_pointwise_transfer, "-t", "--use-pointwise-transfer",
                   "-no-t", "--dont-use-pointwise-transfer",
                   "Use pointwise transfer operators instead of L2 projection.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
    args.ParseCheck();
 
    // Read the mesh from the given mesh file.
@@ -128,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);
 
@@ -165,13 +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; }
 

tests/unit/general/test_umpire_mem.cpp

@@ -225,17 +225,17 @@ static void test_umpire_device_memory()
    SPLIT();
 
    //
-   // Check DeleteDevice with temporary device buffers
+   // Check ReleaseDeviceMemory with temporary device buffers
    //
 
    // remove from temporary memory
    // don't copy to host, verify that the value is still the "host" value
-   host_temp.DeleteDevice(false);
+   host_temp.ReleaseDeviceMemory(false);
    REQUIRE(host_temp[0] == host_val);
    // copy to host, verify that the value is the "device" value
-   dev_temp.DeleteDevice();
+   dev_temp.ReleaseDeviceMemory();
    REQUIRE(dev_temp[0] == dev_val);
-   pinned_host_temp.DeleteDevice();
+   pinned_host_temp.ReleaseDeviceMemory();
 
    printf("Delete all temporary memory: ");
    CHECK_PERM(num_bytes*3);