Merge remote-tracking branch 'origin/hypre-runtime-compute-policy' into apps-dev

examples/ex2p.cpp

@@ -100,6 +100,21 @@ int main(int argc, char *argv[])
    Device device(device_config);
    if (myid == 0) { device.Print(); }
 
+   if (mfem::Device::Allows(mfem::Backend::DEVICE_MASK))
+   {
+      HYPRE_SetMemoryLocation(HYPRE_MEMORY_DEVICE);
+      HYPRE_SetExecutionPolicy(HYPRE_EXEC_DEVICE);
+      HYPRE_DeviceInitialize();
+   }
+   else
+   {
+      HYPRE_SetMemoryLocation(HYPRE_MEMORY_HOST);
+      HYPRE_SetExecutionPolicy(HYPRE_EXEC_HOST);
+   }
+
+   auto loc = mfem::GetHypreMemoryLocation();
+   auto exec = mfem::GetHypreExecutionPolicy();
+
    // 4. Read the (serial) mesh from the given mesh file on all processors.  We
    //    can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
    //    and volume meshes with the same code.

fem/coefficient.cpp

@@ -1592,6 +1592,10 @@ void VectorQuadratureFunctionCoefficient::Eval(Vector &V,
    QuadF.HostRead();
 
    const int el_idx = QuadF.GetSpace()->GetEntityIndex(T);
+   // Handle the case of "interior boundary elements" and FaceQuadratureSpace
+   // with FaceType::Boundary.
+   if (el_idx < 0) { V = 0.0; return; }
+
    const int ip_idx = QuadF.GetSpace()->GetPermutedIndex(el_idx, ip.index);
 
    if (index == 0 && vdim == QuadF.GetVDim())
@@ -1629,6 +1633,9 @@ double QuadratureFunctionCoefficient::Eval(ElementTransformation &T,
    QuadF.HostRead();
    Vector temp(1);
    const int el_idx = QuadF.GetSpace()->GetEntityIndex(T);
+   // Handle the case of "interior boundary elements" and FaceQuadratureSpace
+   // with FaceType::Boundary.
+   if (el_idx < 0) { return 0.0; }
    const int ip_idx = QuadF.GetSpace()->GetPermutedIndex(el_idx, ip.index);
    QuadF.GetValues(el_idx, ip_idx, temp);
    return temp[0];

fem/complex_fem.cpp

@@ -1243,25 +1243,28 @@ ParSesquilinearForm::FormLinearSystem(const Array<int> &ess_tdof_list,
          HypreParMatrix * Ah;
          A_i.Get(Ah);
          hypre_ParCSRMatrix *Aih = *Ah;
-#if !defined(HYPRE_USING_GPU)
-         ess_tdof_list.HostRead();
-         for (int k = 0; k < n; k++)
+         if (!HypreUsingGPU())
          {
-            const int j = ess_tdof_list[k];
-            Aih->diag->data[Aih->diag->i[j]] = 0.0;
+            ess_tdof_list.HostRead();
+            for (int k = 0; k < n; k++)
+            {
+               const int j = ess_tdof_list[k];
+               Aih->diag->data[Aih->diag->i[j]] = 0.0;
+            }
          }
-#else
-         Ah->HypreReadWrite();
-         const int *d_ess_tdof_list =
-            ess_tdof_list.GetMemory().Read(MemoryClass::DEVICE, n);
-         const int *d_diag_i = Aih->diag->i;
-         double *d_diag_data = Aih->diag->data;
-         MFEM_GPU_FORALL(k, n,
+         else
          {
-            const int j = d_ess_tdof_list[k];
-            d_diag_data[d_diag_i[j]] = 0.0;
-         });
-#endif
+            Ah->HypreReadWrite();
+            const int *d_ess_tdof_list =
+               ess_tdof_list.GetMemory().Read(MemoryClass::DEVICE, n);
+            const int *d_diag_i = Aih->diag->i;
+            double *d_diag_data = Aih->diag->data;
+            MFEM_GPU_FORALL(k, n,
+            {
+               const int j = d_ess_tdof_list[k];
+               d_diag_data[d_diag_i[j]] = 0.0;
+            });
+         }
       }
       else
       {

fem/integ/bilininteg_mass_pa.cpp

@@ -128,7 +128,7 @@ void MassIntegrator::AssemblePABoundary(const FiniteElementSpace &fes)
 
    int map_type = el.GetMapType();
    dim = el.GetDim(); // Dimension of the boundary element, *not* the mesh
-   ne = fes.GetMesh()->GetNBE();
+   ne = fes.GetMesh()->GetNFbyType(FaceType::Boundary);
    nq = ir->GetNPoints();
    face_geom = mesh->GetFaceGeometricFactors(*ir, GeometricFactors::DETERMINANTS,
                                              FaceType::Boundary, mt);

fem/lor/lor_ams.cpp

@@ -291,7 +291,8 @@ void BatchedLOR_AMS::FormCoordinateVectors(const Vector &X_vert)
    const auto ltdof_ldof = HypreRead(R->GetMemoryJ());
 
    // Go from E-vector format directly to T-vector format
-   MFEM_HYPRE_FORALL(i, ntdofs,
+   //MFEM_HYPRE_FORALL(i, ntdofs,
+   mfem::forall_switch(HypreUsingGPU(), ntdofs, [=] MFEM_HOST_DEVICE (int i)
    {
       const int j = d_offsets[ltdof_ldof[i]];
       for (int c = 0; c < sdim; ++c)

fem/qspace.cpp

@@ -164,13 +164,20 @@ int FaceQuadratureSpace::GetPermutedIndex(int idx, int iq) const
 
 int FaceQuadratureSpace::GetEntityIndex(const ElementTransformation &T) const
 {
+   auto get_face_index = [this](const int idx)
+   {
+      const auto it = face_indices_inv.find(idx);
+      if (it == face_indices_inv.end()) { return -1; }
+      else { return it->second; }
+   };
+
    switch (T.ElementType)
    {
       case ElementTransformation::FACE:
-         return face_indices_inv.at(T.ElementNo);
+         return get_face_index(T.ElementNo);
       case ElementTransformation::BDR_ELEMENT:
       case ElementTransformation::BDR_FACE:
-         return face_indices_inv.at(mesh.GetBdrElementEdgeIndex(T.ElementNo));
+         return get_face_index(mesh.GetBdrElementEdgeIndex(T.ElementNo));
       default:
          MFEM_ABORT("Invalid element type.");
          return -1;

fem/qspace.hpp

@@ -87,7 +87,8 @@ public:
    ///
    /// For a QuadratureSpace defined on elements, this just returns the element
    /// index. For FaceQuadratureSpace, the returned index depends on the chosen
-   /// FaceType.
+   /// FaceType. If the entity is not found (for example, if @a T represents an
+   /// interior face, and the space has FaceType::Boundary) then -1 is returned.
    virtual int GetEntityIndex(const ElementTransformation &T) const = 0;
 
    /// Write the QuadratureSpace to the stream @a out.

general/mem_manager.hpp

@@ -17,9 +17,6 @@
 #include <cstring> // std::memcpy
 #include <type_traits> // std::is_const
 #include <cstddef> // std::max_align_t
-#ifdef MFEM_USE_MPI
-#include <HYPRE_config.h> // HYPRE_USING_GPU
-#endif
 
 namespace mfem
 {
@@ -994,17 +991,7 @@ inline void Memory<T>::MakeAlias(const Memory &base, int offset, int size)
    h_ptr = base.h_ptr + offset;
    if (!(base.flags & Registered))
    {
-      if (
-#if !defined(HYPRE_USING_GPU)
-         // If the following condition is true then MemoryManager::Exists()
-         // should also be true:
-         IsDeviceMemory(MemoryManager::GetDeviceMemoryType())
-#else
-         // When HYPRE_USING_GPU is defined we always register the 'base' if
-         // the MemoryManager::Exists():
-         MemoryManager::Exists()
-#endif
-      )
+      if (IsDeviceMemory(MemoryManager::GetDeviceMemoryType()))
       {
          // Register 'base':
          MemoryManager::Register_(base.h_ptr, nullptr, base.capacity*sizeof(T),

linalg/hypre.cpp

@@ -121,22 +121,27 @@ bool CanShallowCopy(const Memory<T> &src, MemoryClass mc)
 inline void HypreParVector::_SetDataAndSize_()
 {
    hypre_Vector *x_loc = hypre_ParVectorLocalVector(x);
-#if !defined(HYPRE_USING_GPU)
-   SetDataAndSize(hypre_VectorData(x_loc),
-                  internal::to_int(hypre_VectorSize(x_loc)));
-#else
-   size = internal::to_int(hypre_VectorSize(x_loc));
-   MemoryType mt = (hypre_VectorMemoryLocation(x_loc) == HYPRE_MEMORY_HOST
-                    ? MemoryType::HOST : GetHypreMemoryType());
-   if (hypre_VectorData(x_loc) != NULL)
+#if defined(HYPRE_USING_GPU)
+   if (HypreUsingGPU())
    {
-      data.Wrap(hypre_VectorData(x_loc), size, mt, false);
+      size = internal::to_int(hypre_VectorSize(x_loc));
+      MemoryType mt = (hypre_VectorMemoryLocation(x_loc) == HYPRE_MEMORY_HOST
+                       ? MemoryType::HOST : GetHypreMemoryType());
+      if (hypre_VectorData(x_loc) != NULL)
+      {
+         data.Wrap(hypre_VectorData(x_loc), size, mt, false);
+      }
+      else
+      {
+         data.Reset();
+      }
    }
    else
-   {
-      data.Reset();
-   }
 #endif
+   {
+      SetDataAndSize(hypre_VectorData(x_loc),
+                     internal::to_int(hypre_VectorSize(x_loc)));
+   }
 }
 
 HypreParVector::HypreParVector(MPI_Comm comm, HYPRE_BigInt glob_size,
@@ -316,7 +321,7 @@ void HypreParVector::HypreRead() const
    hypre_VectorData(x_loc) =
       const_cast<double*>(data.Read(GetHypreMemoryClass(), size));
 #ifdef HYPRE_USING_GPU
-   hypre_VectorMemoryLocation(x_loc) = HYPRE_MEMORY_DEVICE;
+   hypre_VectorMemoryLocation(x_loc) = mfem::GetHypreMemoryLocation();
 #endif
 }
 
@@ -325,7 +330,7 @@ void HypreParVector::HypreReadWrite()
    hypre_Vector *x_loc = hypre_ParVectorLocalVector(x);
    hypre_VectorData(x_loc) = data.ReadWrite(GetHypreMemoryClass(), size);
 #ifdef HYPRE_USING_GPU
-   hypre_VectorMemoryLocation(x_loc) = HYPRE_MEMORY_DEVICE;
+   hypre_VectorMemoryLocation(x_loc) = mfem::GetHypreMemoryLocation();
 #endif
 }
 
@@ -334,7 +339,7 @@ void HypreParVector::HypreWrite()
    hypre_Vector *x_loc = hypre_ParVectorLocalVector(x);
    hypre_VectorData(x_loc) = data.Write(GetHypreMemoryClass(), size);
 #ifdef HYPRE_USING_GPU
-   hypre_VectorMemoryLocation(x_loc) = HYPRE_MEMORY_DEVICE;
+   hypre_VectorMemoryLocation(x_loc) = mfem::GetHypreMemoryLocation();
 #endif
 }
 
@@ -348,7 +353,7 @@ void HypreParVector::WrapMemoryRead(const Memory<double> &mem)
    hypre_VectorData(x_loc) =
       const_cast<double*>(mem.Read(GetHypreMemoryClass(), size));
 #ifdef HYPRE_USING_GPU
-   hypre_VectorMemoryLocation(x_loc) = HYPRE_MEMORY_DEVICE;
+   hypre_VectorMemoryLocation(x_loc) = mfem::GetHypreMemoryLocation();
 #endif
    data.MakeAlias(mem, 0, size);
 }
@@ -362,7 +367,7 @@ void HypreParVector::WrapMemoryReadWrite(Memory<double> &mem)
    hypre_Vector *x_loc = hypre_ParVectorLocalVector(x);
    hypre_VectorData(x_loc) = mem.ReadWrite(GetHypreMemoryClass(), size);
 #ifdef HYPRE_USING_GPU
-   hypre_VectorMemoryLocation(x_loc) = HYPRE_MEMORY_DEVICE;
+   hypre_VectorMemoryLocation(x_loc) = mfem::GetHypreMemoryLocation();
 #endif
    data.MakeAlias(mem, 0, size);
 }
@@ -376,7 +381,7 @@ void HypreParVector::WrapMemoryWrite(Memory<double> &mem)
    hypre_Vector *x_loc = hypre_ParVectorLocalVector(x);
    hypre_VectorData(x_loc) = mem.Write(GetHypreMemoryClass(), size);
 #ifdef HYPRE_USING_GPU
-   hypre_VectorMemoryLocation(x_loc) = HYPRE_MEMORY_DEVICE;
+   hypre_VectorMemoryLocation(x_loc) = mfem::GetHypreMemoryLocation();
 #endif
    data.MakeAlias(mem, 0, size);
 }
@@ -536,6 +541,10 @@ void HypreParMatrix::Init()
 
 void HypreParMatrix::Read(MemoryClass mc) const
 {
+   if (GetHypreMemoryLocation() == HYPRE_MEMORY_HOST)
+   {
+      mc = Device::GetHostMemoryClass();
+   }
    hypre_CSRMatrix *diag = hypre_ParCSRMatrixDiag(A);
    hypre_CSRMatrix *offd = hypre_ParCSRMatrixOffd(A);
    const int num_rows = NumRows();
@@ -548,8 +557,8 @@ void HypreParMatrix::Read(MemoryClass mc) const
    offd->j = const_cast<HYPRE_Int*>(mem_offd.J.Read(mc, offd_nnz));
    offd->data = const_cast<double*>(mem_offd.data.Read(mc, offd_nnz));
 #if MFEM_HYPRE_VERSION >= 21800
-   decltype(diag->memory_location) ml =
-      (mc != GetHypreMemoryClass() ? HYPRE_MEMORY_HOST : HYPRE_MEMORY_DEVICE);
+   decltype(diag->memory_location) ml = (mc == MemoryClass::HOST) ?
+                                        HYPRE_MEMORY_HOST : GetHypreMemoryLocation();
    diag->memory_location = ml;
    offd->memory_location = ml;
 #endif
@@ -557,6 +566,10 @@ void HypreParMatrix::Read(MemoryClass mc) const
 
 void HypreParMatrix::ReadWrite(MemoryClass mc)
 {
+   if (GetHypreMemoryLocation() == HYPRE_MEMORY_HOST)
+   {
+      mc = Device::GetHostMemoryClass();
+   }
    hypre_CSRMatrix *diag = hypre_ParCSRMatrixDiag(A);
    hypre_CSRMatrix *offd = hypre_ParCSRMatrixOffd(A);
    const int num_rows = NumRows();
@@ -569,8 +582,8 @@ void HypreParMatrix::ReadWrite(MemoryClass mc)
    offd->j = mem_offd.J.ReadWrite(mc, offd_nnz);
    offd->data = mem_offd.data.ReadWrite(mc, offd_nnz);
 #if MFEM_HYPRE_VERSION >= 21800
-   decltype(diag->memory_location) ml =
-      (mc != GetHypreMemoryClass() ? HYPRE_MEMORY_HOST : HYPRE_MEMORY_DEVICE);
+   decltype(diag->memory_location) ml = (mc == MemoryClass::HOST) ?
+                                        HYPRE_MEMORY_HOST : GetHypreMemoryLocation();
    diag->memory_location = ml;
    offd->memory_location = ml;
 #endif
@@ -578,6 +591,10 @@ void HypreParMatrix::ReadWrite(MemoryClass mc)
 
 void HypreParMatrix::Write(MemoryClass mc, bool set_diag, bool set_offd)
 {
+   if (GetHypreMemoryLocation() == HYPRE_MEMORY_HOST)
+   {
+      mc = Device::GetHostMemoryClass();
+   }
    hypre_CSRMatrix *diag = hypre_ParCSRMatrixDiag(A);
    hypre_CSRMatrix *offd = hypre_ParCSRMatrixOffd(A);
    if (set_diag)
@@ -593,8 +610,8 @@ void HypreParMatrix::Write(MemoryClass mc, bool set_diag, bool set_offd)
       offd->data = mem_offd.data.Write(mc, mem_offd.data.Capacity());
    }
 #if MFEM_HYPRE_VERSION >= 21800
-   decltype(diag->memory_location) ml =
-      (mc != GetHypreMemoryClass() ? HYPRE_MEMORY_HOST : HYPRE_MEMORY_DEVICE);
+   decltype(diag->memory_location) ml = (mc == MemoryClass::HOST) ?
+                                        HYPRE_MEMORY_HOST : GetHypreMemoryLocation();
    if (set_diag) { diag->memory_location = ml; }
    if (set_offd) { offd->memory_location = ml; }
 #endif
@@ -760,7 +777,7 @@ signed char HypreParMatrix::HypreCsrToMem(hypre_CSRMatrix *h_mat,
       h_mat->data = mem.data.ReadWrite(hypre_mc, nnz);
       h_mat->owns_data = 0;
 #if MFEM_HYPRE_VERSION >= 21800
-      h_mat->memory_location = HYPRE_MEMORY_DEVICE;
+      h_mat->memory_location = mfem::GetHypreMemoryLocation();
 #endif
       return 3;
    }
@@ -1361,8 +1378,11 @@ hypre_ParCSRMatrix* HypreParMatrix::StealData()
    MFEM_ASSERT(ParCSROwner, "");
    hypre_ParCSRMatrix *R = A;
 #ifdef HYPRE_USING_GPU
-   if (diagOwner == -1) { HostReadWrite(); }
-   else { HypreReadWrite(); }
+   if (HypreUsingGPU())
+   {
+      if (diagOwner == -1) { HostReadWrite(); }
+      else { HypreReadWrite(); }
+   }
 #endif
    ParCSROwner = false;
    Destroy();
@@ -1710,7 +1730,10 @@ void HypreParMatrix::EnsureMultTranspose() const
 #if (MFEM_HYPRE_VERSION == 22500 && HYPRE_DEVELOP_NUMBER >= 1) || \
     (MFEM_HYPRE_VERSION > 22500)
 #ifdef HYPRE_USING_GPU
-   hypre_ParCSRMatrixLocalTranspose(A);
+   if (HypreUsingGPU())
+   {
+      hypre_ParCSRMatrixLocalTranspose(A);
+   }
 #endif
 #endif
 }
@@ -1720,15 +1743,18 @@ void HypreParMatrix::ResetTranspose() const
 #if (MFEM_HYPRE_VERSION == 22500 && HYPRE_DEVELOP_NUMBER >= 1) || \
     (MFEM_HYPRE_VERSION > 22500)
 #ifdef HYPRE_USING_GPU
-   if (A->diagT)
+   if (HypreUsingGPU())
    {
-      hypre_CSRMatrixDestroy(A->diagT);
-      A->diagT = NULL;
-   }
-   if (A->offdT)
-   {
-      hypre_CSRMatrixDestroy(A->offdT);
-      A->offdT = NULL;
+      if (A->diagT)
+      {
+         hypre_CSRMatrixDestroy(A->diagT);
+         A->diagT = NULL;
+      }
+      if (A->offdT)
+      {
+         hypre_CSRMatrixDestroy(A->offdT);
+         A->offdT = NULL;
+      }
    }
 #endif
 #endif
@@ -2417,8 +2443,10 @@ void HypreParMatrix::EliminateBC(const Array<int> &ess_dofs,
       }
 
       // Which of the local rows are to be eliminated?
-      MFEM_HYPRE_FORALL(i, diag_nrows, eliminate_row[i] = 0; );
-      MFEM_HYPRE_FORALL(i, n_ess_dofs, eliminate_row[ess_dofs_d[i]] = 1; );
+      //MFEM_HYPRE_FORALL(i, diag_nrows, eliminate_row[i] = 0; ); // TODO TMS
+      mfem::forall_switch(HypreUsingGPU(), diag_nrows, [=] MFEM_HOST_DEVICE (int i) { eliminate_row[i] = 0; });
+      //MFEM_HYPRE_FORALL(i, n_ess_dofs, eliminate_row[ess_dofs_d[i]] = 1; ); // TODO TMS
+      mfem::forall_switch(HypreUsingGPU(), n_ess_dofs, [=] MFEM_HOST_DEVICE (int i) { eliminate_row[ess_dofs_d[i]] = 1; });
 
       // Use a matvec communication pattern to find (in eliminate_col) which of
       // the local offd columns are to be eliminated
@@ -2429,26 +2457,37 @@ void HypreParMatrix::EliminateBC(const Array<int> &ess_dofs,
 
       HYPRE_Int *send_map_elmts;
 #if defined(HYPRE_USING_GPU)
-      hypre_ParCSRCommPkgCopySendMapElmtsToDevice(comm_pkg);
-      send_map_elmts = hypre_ParCSRCommPkgDeviceSendMapElmts(comm_pkg);
-#else
-      send_map_elmts = hypre_ParCSRCommPkgSendMapElmts(comm_pkg);
+      if (HypreUsingGPU())
+      {
+         hypre_ParCSRCommPkgCopySendMapElmtsToDevice(comm_pkg);
+         send_map_elmts = hypre_ParCSRCommPkgDeviceSendMapElmts(comm_pkg);
+      }
+      else
 #endif
-      MFEM_HYPRE_FORALL(i, int_buf_sz,
+      {
+         send_map_elmts = hypre_ParCSRCommPkgSendMapElmts(comm_pkg);
+      }
+      //MFEM_HYPRE_FORALL(i, int_buf_sz,
+      mfem::forall_switch(HypreUsingGPU(), int_buf_sz, [=] MFEM_HOST_DEVICE (int i)
       {
          int k = send_map_elmts[i];
          int_buf_data[i] = eliminate_row[k];
       });
 
 #if defined(HYPRE_USING_GPU)
-      // Try to use device-aware MPI for the communication if available
-      comm_handle = hypre_ParCSRCommHandleCreate_v2(
-                       11, comm_pkg, HYPRE_MEMORY_DEVICE, int_buf_data,
-                       HYPRE_MEMORY_DEVICE, eliminate_col);
-#else
-      comm_handle = hypre_ParCSRCommHandleCreate(
-                       11, comm_pkg, int_buf_data, eliminate_col );
+      if (HypreUsingGPU())
+      {
+         // Try to use device-aware MPI for the communication if available
+         comm_handle = hypre_ParCSRCommHandleCreate_v2(
+                          11, comm_pkg, HYPRE_MEMORY_DEVICE, int_buf_data,
+                          HYPRE_MEMORY_DEVICE, eliminate_col);
+      }
+      else
 #endif
+      {
+         comm_handle = hypre_ParCSRCommHandleCreate(
+                          11, comm_pkg, int_buf_data, eliminate_col );
+      }
    }
 
    // Eliminate rows and columns in the diagonal block
@@ -2457,7 +2496,8 @@ void HypreParMatrix::EliminateBC(const Array<int> &ess_dofs,
       const auto J = diag->j;
       auto data = diag->data;
 
-      MFEM_HYPRE_FORALL(i, n_ess_dofs,
+      //MFEM_HYPRE_FORALL(i, n_ess_dofs,
+      mfem::forall_switch(HypreUsingGPU(), n_ess_dofs, [=] MFEM_HOST_DEVICE (int i)
       {
          const int idof = ess_dofs_d[i];
          for (int j=I[idof]; j<I[idof+1]; ++j)
@@ -2495,7 +2535,8 @@ void HypreParMatrix::EliminateBC(const Array<int> &ess_dofs,
    {
       const auto I = offd->i;
       auto data = offd->data;
-      MFEM_HYPRE_FORALL(i, n_ess_dofs,
+      //MFEM_HYPRE_FORALL(i, n_ess_dofs,
+      mfem::forall_switch(HypreUsingGPU(), n_ess_dofs, [=] MFEM_HOST_DEVICE (int i)
       {
          const int idof = ess_dofs_d[i];
          for (int j=I[idof]; j<I[idof+1]; ++j)
@@ -2516,7 +2557,8 @@ void HypreParMatrix::EliminateBC(const Array<int> &ess_dofs,
       const auto I = offd->i;
       const auto J = offd->j;
       auto data = offd->data;
-      MFEM_HYPRE_FORALL(i, nrows_offd,
+      //MFEM_HYPRE_FORALL(i, nrows_offd,
+      mfem::forall_switch(HypreUsingGPU(), nrows_offd, [=] MFEM_HOST_DEVICE (int i)
       {
          for (int j=I[i]; j<I[i+1]; ++j)
          {
@@ -2692,7 +2734,7 @@ void HypreParMatrix::Destroy()
    if (A == NULL) { return; }
 
 #ifdef HYPRE_USING_GPU
-   if (ParCSROwner && (diagOwner < 0 || offdOwner < 0))
+   if (HypreUsingGPU() && ParCSROwner && (diagOwner < 0 || offdOwner < 0))
    {
       // Put the "host" or "hypre" pointers in {i,j,data} of A->{diag,offd}, so
       // that they can be destroyed by hypre when hypre_ParCSRMatrixDestroy(A)
@@ -2856,10 +2898,15 @@ HypreParMatrix * ParMult(const HypreParMatrix *A, const HypreParMatrix *B,
 {
    hypre_ParCSRMatrix * ab;
 #ifdef HYPRE_USING_GPU
-   ab = hypre_ParCSRMatMat(*A, *B);
-#else
-   ab = hypre_ParMatmul(*A,*B);
+   if (HypreUsingGPU())
+   {
+      ab = hypre_ParCSRMatMat(*A, *B);
+   }
+   else
 #endif
+   {
+      ab = hypre_ParMatmul(*A,*B);
+   }
    hypre_ParCSRMatrixSetNumNonzeros(ab);
 
    if (!hypre_ParCSRMatrixCommPkg(ab)) { hypre_MatvecCommPkgCreate(ab); }
@@ -2883,6 +2930,7 @@ HypreParMatrix * RAP(const HypreParMatrix *A, const HypreParMatrix *P)
    //        in ex28p.
    // Quick fix: add a diagonal matrix with 0 diagonal.
    // Maybe use hypre_CSRMatrixCheckDiagFirst to see if we need the fix.
+   if (HypreUsingGPU())
    {
       hypre_ParCSRMatrix *Q = hypre_ParCSRMatMat(*A,*P);
       const bool keepTranspose = false;
@@ -2892,25 +2940,27 @@ HypreParMatrix * RAP(const HypreParMatrix *A, const HypreParMatrix *P)
       // alternative:
       // hypre_ParCSRMatrixRAPKT
    }
-#else
-#if MFEM_HYPRE_VERSION <= 22200
-   HYPRE_Int P_owns_its_col_starts =
-      hypre_ParCSRMatrixOwnsColStarts((hypre_ParCSRMatrix*)(*P));
+   else
 #endif
-
-   hypre_BoomerAMGBuildCoarseOperator(*P,*A,*P,&rap);
-
-#if MFEM_HYPRE_VERSION <= 22200
-   /* Warning: hypre_BoomerAMGBuildCoarseOperator steals the col_starts
-      from P (even if it does not own them)! */
-   hypre_ParCSRMatrixSetRowStartsOwner(rap,0);
-   hypre_ParCSRMatrixSetColStartsOwner(rap,0);
-   if (P_owns_its_col_starts)
    {
-      hypre_ParCSRMatrixSetColStartsOwner(*P, 1);
+#if MFEM_HYPRE_VERSION <= 22200
+      HYPRE_Int P_owns_its_col_starts =
+         hypre_ParCSRMatrixOwnsColStarts((hypre_ParCSRMatrix*)(*P));
+#endif
+
+      hypre_BoomerAMGBuildCoarseOperator(*P,*A,*P,&rap);
+
+#if MFEM_HYPRE_VERSION <= 22200
+      /* Warning: hypre_BoomerAMGBuildCoarseOperator steals the col_starts
+         from P (even if it does not own them)! */
+      hypre_ParCSRMatrixSetRowStartsOwner(rap,0);
+      hypre_ParCSRMatrixSetColStartsOwner(rap,0);
+      if (P_owns_its_col_starts)
+      {
+         hypre_ParCSRMatrixSetColStartsOwner(*P, 1);
+      }
+#endif
    }
-#endif
-#endif
 
    hypre_ParCSRMatrixSetNumNonzeros(rap);
    // hypre_MatvecCommPkgCreate(rap);
@@ -2924,36 +2974,39 @@ HypreParMatrix * RAP(const HypreParMatrix * Rt, const HypreParMatrix *A,
    hypre_ParCSRMatrix * rap;
 
 #ifdef HYPRE_USING_GPU
+   if (HypreUsingGPU())
    {
       hypre_ParCSRMatrix *Q = hypre_ParCSRMatMat(*A,*P);
       rap = hypre_ParCSRTMatMat(*Rt,Q);
       hypre_ParCSRMatrixDestroy(Q);
    }
-#else
+   else
+#endif
+   {
 #if MFEM_HYPRE_VERSION <= 22200
-   HYPRE_Int P_owns_its_col_starts =
-      hypre_ParCSRMatrixOwnsColStarts((hypre_ParCSRMatrix*)(*P));
-   HYPRE_Int Rt_owns_its_col_starts =
-      hypre_ParCSRMatrixOwnsColStarts((hypre_ParCSRMatrix*)(*Rt));
+      HYPRE_Int P_owns_its_col_starts =
+         hypre_ParCSRMatrixOwnsColStarts((hypre_ParCSRMatrix*)(*P));
+      HYPRE_Int Rt_owns_its_col_starts =
+         hypre_ParCSRMatrixOwnsColStarts((hypre_ParCSRMatrix*)(*Rt));
 #endif
 
-   hypre_BoomerAMGBuildCoarseOperator(*Rt,*A,*P,&rap);
+      hypre_BoomerAMGBuildCoarseOperator(*Rt,*A,*P,&rap);
 
 #if MFEM_HYPRE_VERSION <= 22200
-   /* Warning: hypre_BoomerAMGBuildCoarseOperator steals the col_starts
-      from Rt and P (even if they do not own them)! */
-   hypre_ParCSRMatrixSetRowStartsOwner(rap,0);
-   hypre_ParCSRMatrixSetColStartsOwner(rap,0);
-   if (P_owns_its_col_starts)
-   {
-      hypre_ParCSRMatrixSetColStartsOwner(*P, 1);
-   }
-   if (Rt_owns_its_col_starts)
-   {
-      hypre_ParCSRMatrixSetColStartsOwner(*Rt, 1);
-   }
-#endif
+      /* Warning: hypre_BoomerAMGBuildCoarseOperator steals the col_starts
+         from Rt and P (even if they do not own them)! */
+      hypre_ParCSRMatrixSetRowStartsOwner(rap,0);
+      hypre_ParCSRMatrixSetColStartsOwner(rap,0);
+      if (P_owns_its_col_starts)
+      {
+         hypre_ParCSRMatrixSetColStartsOwner(*P, 1);
+      }
+      if (Rt_owns_its_col_starts)
+      {
+         hypre_ParCSRMatrixSetColStartsOwner(*Rt, 1);
+      }
 #endif
+   }
 
    hypre_ParCSRMatrixSetNumNonzeros(rap);
    // hypre_MatvecCommPkgCreate(rap);
@@ -3383,7 +3436,7 @@ int ParCSRRelax_FIR(hypre_ParCSRMatrix *A, // matrix to relax with
 
 HypreSmoother::HypreSmoother() : Solver()
 {
-   type = default_type;
+   type = DefaultType();
    relax_times = 1;
    relax_weight = 1.0;
    omega = 1.0;
@@ -3408,7 +3461,7 @@ HypreSmoother::HypreSmoother(const HypreParMatrix &A_, int type_,
                              double omega_, int poly_order_,
                              double poly_fraction_, int eig_est_cg_iter_)
 {
-   type = type_;
+   type = (type_ == -1) ? DefaultType() : type_;
    relax_times = relax_times_;
    relax_weight = relax_weight_;
    omega = omega_;
@@ -3523,17 +3576,22 @@ void HypreSmoother::SetOperator(const Operator &op)
    if (l1_norms && pos_l1_norms)
    {
 #if defined(HYPRE_USING_GPU)
-      double *d_l1_norms = l1_norms;  // avoid *this capture
-      MFEM_GPU_FORALL(i, height,
+      if (HypreUsingGPU())
       {
-         d_l1_norms[i] = std::abs(d_l1_norms[i]);
-      });
-#else
-      for (int i = 0; i < height; i++)
-      {
-         l1_norms[i] = std::abs(l1_norms[i]);
+         double *d_l1_norms = l1_norms;  // avoid *this capture
+         MFEM_GPU_FORALL(i, height,
+         {
+            d_l1_norms[i] = std::abs(d_l1_norms[i]);
+         });
       }
+      else
 #endif
+      {
+         for (int i = 0; i < height; i++)
+         {
+            l1_norms[i] = std::abs(l1_norms[i]);
+         }
+      }
    }
 
    if (type == 16)
@@ -4809,41 +4867,49 @@ HypreBoomerAMG::HypreBoomerAMG(const HypreParMatrix &A) : HypreSolver(&A)
 
 void HypreBoomerAMG::SetDefaultOptions()
 {
-#if !defined(HYPRE_USING_GPU)
-   // AMG coarsening options:
-   int coarsen_type = 10;   // 10 = HMIS, 8 = PMIS, 6 = Falgout, 0 = CLJP
-   int agg_levels   = 1;    // number of aggressive coarsening levels
-   double theta     = 0.25; // strength threshold: 0.25, 0.5, 0.8
-
    // AMG interpolation options:
-   int interp_type  = 6;    // 6 = extended+i, 0 = classical
-   int Pmax         = 4;    // max number of elements per row in P
+   int coarsen_type, agg_levels, interp_type, Pmax, relax_type, relax_sweeps,
+       print_level, max_levels;
+   double theta;
 
-   // AMG relaxation options:
-   int relax_type   = 8;    // 8 = l1-GS, 6 = symm. GS, 3 = GS, 18 = l1-Jacobi
-   int relax_sweeps = 1;    // relaxation sweeps on each level
+   if (!HypreUsingGPU())
+   {
+      // AMG coarsening options:
+      coarsen_type = 10;   // 10 = HMIS, 8 = PMIS, 6 = Falgout, 0 = CLJP
+      agg_levels   = 1;    // number of aggressive coarsening levels
+      theta     = 0.25; // strength threshold: 0.25, 0.5, 0.8
 
-   // Additional options:
-   int print_level  = 1;    // print AMG iterations? 1 = no, 2 = yes
-   int max_levels   = 25;   // max number of levels in AMG hierarchy
-#else
-   // AMG coarsening options:
-   int coarsen_type = 8;    // 10 = HMIS, 8 = PMIS, 6 = Falgout, 0 = CLJP
-   int agg_levels   = 0;    // number of aggressive coarsening levels
-   double theta     = 0.25; // strength threshold: 0.25, 0.5, 0.8
+      // AMG interpolation options:
+      interp_type  = 6;    // 6 = extended+i, 0 = classical
+      Pmax         = 4;    // max number of elements per row in P
 
-   // AMG interpolation options:
-   int interp_type  = 6;    // 6 = extended+i, or 18 = extended+e
-   int Pmax         = 4;    // max number of elements per row in P
+      // AMG relaxation options:
+      relax_type   = 8;    // 8 = l1-GS, 6 = symm. GS, 3 = GS, 18 = l1-Jacobi
+      relax_sweeps = 1;    // relaxation sweeps on each level
 
-   // AMG relaxation options:
-   int relax_type   = 18;   // 18 = l1-Jacobi, or 16 = Chebyshev
-   int relax_sweeps = 1;    // relaxation sweeps on each level
+      // Additional options:
+      print_level  = 1;    // print AMG iterations? 1 = no, 2 = yes
+      max_levels   = 25;   // max number of levels in AMG hierarchy
+   }
+   else
+   {
+      // AMG coarsening options:
+      coarsen_type = 8;    // 10 = HMIS, 8 = PMIS, 6 = Falgout, 0 = CLJP
+      agg_levels   = 0;    // number of aggressive coarsening levels
+      theta     = 0.25; // strength threshold: 0.25, 0.5, 0.8
 
-   // Additional options:
-   int print_level  = 1;    // print AMG iterations? 1 = no, 2 = yes
-   int max_levels   = 25;   // max number of levels in AMG hierarchy
-#endif
+      // AMG interpolation options:
+      interp_type  = 6;    // 6 = extended+i, or 18 = extended+e
+      Pmax         = 4;    // max number of elements per row in P
+
+      // AMG relaxation options:
+      relax_type   = 18;   // 18 = l1-Jacobi, or 16 = Chebyshev
+      relax_sweeps = 1;    // relaxation sweeps on each level
+
+      // Additional options:
+      print_level  = 1;    // print AMG iterations? 1 = no, 2 = yes
+      max_levels   = 25;   // max number of levels in AMG hierarchy
+   }
 
    HYPRE_BoomerAMGSetCoarsenType(amg_precond, coarsen_type);
    HYPRE_BoomerAMGSetAggNumLevels(amg_precond, agg_levels);
@@ -4980,14 +5046,20 @@ void HypreBoomerAMG::SetSystemsOptions(int dim, bool order_bynodes)
       // After the addition of hypre_IntArray, mapping is assumed
       // to be a device pointer. Previously, it was assumed to be
       // a host pointer.
+      HYPRE_Int *mapping = nullptr;
 #if defined(hypre_IntArrayData) && defined(HYPRE_USING_GPU)
-      HYPRE_Int *mapping = mfem_hypre_CTAlloc(HYPRE_Int, height);
-      hypre_TMemcpy(mapping, h_mapping, HYPRE_Int, height,
-                    HYPRE_MEMORY_DEVICE, HYPRE_MEMORY_HOST);
-      mfem_hypre_TFree_host(h_mapping);
-#else
-      HYPRE_Int *mapping = h_mapping;
+      if (HypreUsingGPU())
+      {
+         mapping = mfem_hypre_CTAlloc(HYPRE_Int, height);
+         hypre_TMemcpy(mapping, h_mapping, HYPRE_Int, height,
+                       HYPRE_MEMORY_DEVICE, HYPRE_MEMORY_HOST);
+         mfem_hypre_TFree_host(h_mapping);
+      }
+      else
 #endif
+      {
+         mapping = h_mapping;
+      }
 
       // hypre actually deletes the mapping pointer in HYPRE_BoomerAMGDestroy,
       // so we don't need to track it
@@ -5079,7 +5151,10 @@ void HypreBoomerAMG::RecomputeRBMs()
 void HypreBoomerAMG::SetElasticityOptions(ParFiniteElementSpace *fespace_)
 {
 #ifdef HYPRE_USING_GPU
-   MFEM_ABORT("this method is not supported in hypre built with GPU support");
+   if (HypreUsingGPU())
+   {
+      MFEM_ABORT("this method is not supported in hypre built with GPU support");
+   }
 #endif
 
    // Save the finite element space to support multiple calls to SetOperator()
@@ -5315,23 +5390,14 @@ void HypreAMS::MakeSolver(int sdim, int cycle_type)
    int rlx_sweeps       = 1;
    double rlx_weight    = 1.0;
    double rlx_omega     = 1.0;
-#if !defined(HYPRE_USING_GPU)
-   int amg_coarsen_type = 10;
-   int amg_agg_levels   = 1;
-   int amg_rlx_type     = 8;
-   int rlx_type         = 2;
+   const bool hypre_gpu = HypreUsingGPU();
+   int amg_coarsen_type = hypre_gpu ? 8 : 10;
+   int amg_agg_levels   = hypre_gpu ? 0 : 1;
+   int amg_rlx_type     = hypre_gpu ? 18 : 8;
+   int rlx_type         = hypre_gpu ? 1: 2;
    double theta         = 0.25;
    int amg_interp_type  = 6;
    int amg_Pmax         = 4;
-#else
-   int amg_coarsen_type = 8;
-   int amg_agg_levels   = 0;
-   int amg_rlx_type     = 18;
-   int rlx_type         = 1;
-   double theta         = 0.25;
-   int amg_interp_type  = 6;
-   int amg_Pmax         = 4;
-#endif
 
    space_dim = sdim;
    ams_cycle_type = cycle_type;
@@ -5692,23 +5758,14 @@ void HypreADS::MakeSolver()
    int rlx_sweeps       = 1;
    double rlx_weight    = 1.0;
    double rlx_omega     = 1.0;
-#if !defined(HYPRE_USING_GPU)
-   int rlx_type         = 2;
-   int amg_coarsen_type = 10;
-   int amg_agg_levels   = 1;
-   int amg_rlx_type     = 8;
+   const bool hypre_gpu = HypreUsingGPU();
+   int rlx_type         = hypre_gpu ? 1 : 2;
+   int amg_coarsen_type = hypre_gpu ? 8 : 10;
+   int amg_agg_levels   = hypre_gpu ? 0 : 1;
+   int amg_rlx_type     = hypre_gpu ? 18 : 8;
    double theta         = 0.25;
    int amg_interp_type  = 6;
    int amg_Pmax         = 4;
-#else
-   int rlx_type         = 1;
-   int amg_coarsen_type = 8;
-   int amg_agg_levels   = 0;
-   int amg_rlx_type     = 18;
-   double theta         = 0.25;
-   int amg_interp_type  = 6;
-   int amg_Pmax         = 4;
-#endif
 
    HYPRE_ADSCreate(&ads);
 

linalg/hypre.hpp

@@ -46,15 +46,15 @@
 // MFEM_HYPRE_FORALL is a macro similar to mfem::forall, but it executes on the
 // device that hypre was configured with (no matter what device was selected
 // in MFEM's runtime configuration).
-#if defined(HYPRE_USING_CUDA)
-#define MFEM_HYPRE_FORALL(i, N,...) CuWrap1D(N, [=] MFEM_DEVICE      \
-                                       (int i) {__VA_ARGS__})
-#elif defined(HYPRE_USING_HIP)
-#define MFEM_HYPRE_FORALL(i, N,...) HipWrap1D(N, [=] MFEM_DEVICE     \
-                                        (int i) {__VA_ARGS__})
-#else
-#define MFEM_HYPRE_FORALL(i, N,...) for (int i = 0; i < N; i++) { __VA_ARGS__ }
-#endif
+//#if defined(HYPRE_USING_CUDA)
+//#define MFEM_HYPRE_FORALL(i, N,...) CuWrap1D(N, [=] MFEM_DEVICE      \
+//                                       (int i) {__VA_ARGS__})
+//#elif defined(HYPRE_USING_HIP)
+//#define MFEM_HYPRE_FORALL(i, N,...) HipWrap1D(N, [=] MFEM_DEVICE     \
+//                                        (int i) {__VA_ARGS__})
+//#else
+//#define MFEM_HYPRE_FORALL(i, N,...) for (int i = 0; i < N; i++) { __VA_ARGS__ }
+//#endif
 
 #include "sparsemat.hpp"
 #include "hypre_parcsr.hpp"
@@ -134,14 +134,15 @@ inline int to_int(HYPRE_Int i)
 
 
 /// The MemoryClass used by Hypre objects.
-inline constexpr MemoryClass GetHypreMemoryClass()
+inline MemoryClass GetHypreMemoryClass()
 {
 #if !defined(HYPRE_USING_GPU)
    return MemoryClass::HOST;
 #elif defined(HYPRE_USING_UNIFIED_MEMORY)
    return MemoryClass::MANAGED;
 #else
-   return MemoryClass::DEVICE;
+   return (GetHypreMemoryLocation() == HYPRE_MEMORY_DEVICE) ? MemoryClass::DEVICE :
+          MemoryClass::HOST;
 #endif
 }
 
@@ -153,7 +154,19 @@ inline MemoryType GetHypreMemoryType()
 #elif defined(HYPRE_USING_UNIFIED_MEMORY)
    return MemoryType::MANAGED;
 #else
-   return MemoryType::DEVICE;
+   return (GetHypreMemoryLocation() == HYPRE_MEMORY_DEVICE) ? MemoryType::DEVICE :
+          Device::GetHostMemoryType();
+#endif
+}
+
+inline bool HypreUsingGPU()
+{
+#ifdef HYPRE_USING_GPU
+   HYPRE_MemoryLocation loc;
+   HYPRE_GetMemoryLocation(&loc);
+   return loc == HYPRE_MEMORY_DEVICE;
+#else
+   return false;
 #endif
 }
 
@@ -1026,6 +1039,7 @@ protected:
 
 public:
    /** Hypre smoother types:
+      -1    = Undefined, replaced with DefaultType() in constructors
        0    = Jacobi
        1    = l1-scaled Jacobi
        2    = l1-scaled block Gauss-Seidel/SSOR
@@ -1036,18 +1050,18 @@ public:
        16   = Chebyshev
        1001 = Taubin polynomial smoother
        1002 = FIR polynomial smoother. */
-   enum Type { Jacobi = 0, l1Jacobi = 1, l1GS = 2, l1GStr = 4, lumpedJacobi = 5,
+   enum Type { Undefined = -1, Jacobi = 0, l1Jacobi = 1, l1GS = 2, l1GStr = 4, lumpedJacobi = 5,
                GS = 6, OPFS = 10, Chebyshev = 16, Taubin = 1001, FIR = 1002
              };
-#if !defined(HYPRE_USING_GPU)
-   static constexpr Type default_type = l1GS;
-#else
-   static constexpr Type default_type = l1Jacobi;
-#endif
+
+   Type DefaultType()
+   {
+      return HypreUsingGPU() ? l1Jacobi : l1GS;
+   }
 
    HypreSmoother();
 
-   HypreSmoother(const HypreParMatrix &A_, int type = default_type,
+   HypreSmoother(const HypreParMatrix &A_, int type = Undefined,
                  int relax_times = 1, double relax_weight = 1.0,
                  double omega = 1.0, int poly_order = 2,
                  double poly_fraction = .3, int eig_est_cg_iter = 10);

linalg/hypre_parcsr.hpp

@@ -29,6 +29,34 @@ typedef HYPRE_Int HYPRE_BigInt;
 #define HYPRE_MPI_BIG_INT HYPRE_MPI_INT
 #endif
 
+
+namespace mfem
+{
+#if MFEM_HYPRE_VERSION >= 22600
+inline HYPRE_MemoryLocation GetHypreMemoryLocation()
+{
+   HYPRE_MemoryLocation loc;
+   HYPRE_GetMemoryLocation(&loc);
+   return loc;
+}
+inline HYPRE_ExecutionPolicy GetHypreExecutionPolicy()
+{
+   HYPRE_ExecutionPolicy pol;
+   HYPRE_GetExecutionPolicy(&pol);
+   return pol;
+}
+#else
+inline HYPRE_MemoryLocation GetHypreMemoryLocation()
+{
+#ifdef HYPRE_USING_GPU
+   return HYPRE_MEMORY_DEVICE;
+#else
+   return HYPRE_MEMORY_HOST;
+#endif // HYPRE_USING_GPU
+}
+#endif
+};
+
 // Define macro wrappers for hypre_TAlloc, hypre_CTAlloc and hypre_TFree:
 // mfem_hypre_TAlloc, mfem_hypre_CTAlloc, and mfem_hypre_TFree, respectively.
 // Note: these macros are used in hypre.cpp, hypre_parcsr.cpp, and perhaps
@@ -46,10 +74,10 @@ typedef HYPRE_Int HYPRE_BigInt;
 #else // MFEM_HYPRE_VERSION >= 21400
 
 #define mfem_hypre_TAlloc(type, size) \
-   hypre_TAlloc(type, size, HYPRE_MEMORY_DEVICE)
+   hypre_TAlloc(type, size, mfem::GetHypreMemoryLocation())
 #define mfem_hypre_CTAlloc(type, size) \
-   hypre_CTAlloc(type, size, HYPRE_MEMORY_DEVICE)
-#define mfem_hypre_TFree(ptr) hypre_TFree(ptr, HYPRE_MEMORY_DEVICE)
+   hypre_CTAlloc(type, size, mfem::GetHypreMemoryLocation())
+#define mfem_hypre_TFree(ptr) hypre_TFree(ptr, mfem::GetHypreMemoryLocation())
 
 #define mfem_hypre_TAlloc_host(type, size) \
    hypre_TAlloc(type, size, HYPRE_MEMORY_HOST)

miniapps/hdiv-linear-solver/discrete_divergence.cpp

@@ -85,7 +85,8 @@ void EliminateColumns(HypreParMatrix &D, const Array<int> &ess_dofs)
       const auto I = diag->i;
       const auto J = diag->j;
       auto data = diag->data;
-      MFEM_HYPRE_FORALL(i, nrows_diag,
+      //MFEM_HYPRE_FORALL(i, nrows_diag,
+      mfem::forall_switch(HypreUsingGPU(), nrows_diag, [=] MFEM_HOST_DEVICE (int i)
       {
          for (int jj=I[i]; jj<I[i+1]; ++jj)
          {
@@ -109,7 +110,8 @@ void EliminateColumns(HypreParMatrix &D, const Array<int> &ess_dofs)
       const auto I = offd->i;
       const auto J = offd->j;
       auto data = offd->data;
-      MFEM_HYPRE_FORALL(i, nrows_offd,
+      //MFEM_HYPRE_FORALL(i, nrows_offd,
+      mfem::forall_switch(HypreUsingGPU(), nrows_offd, [=] MFEM_HOST_DEVICE (int i)
       {
          for (int jj=I[i]; jj<I[i+1]; ++jj)
          {