AlgebraicCeedMultigrid: allow custom AmgX config in constructor

AlgebraicSpaceHierarchy.

This should allow more dynamic decisions about amount to coarsen...

CHANGELOG

@@ -167,6 +167,12 @@ Linear and nonlinear solvers
   matrix with the function HypreParMatrixFromBlocks. This could be useful for
   solving block systems with parallel direct solvers such as STRUMPACK.
 
+- Added AlgebraicCeedSolver that does matrix-free algebraic p-multigrid for
+  diffusion problems with the Ceed backend.
+
+- Added CUDA support for SUNDIALS ODE integrators. See the updated SUNDIALS
+  modification of Example 9/9p.
+
 - Added wrappers for hypre's flexible GMRES solver and the new parallel ILU
   preconditioner. The latter requires hypre version 2.19.0 or later.
 

examples/ex1.cpp

@@ -35,6 +35,7 @@
 //               ex1 -pa -d raja-omp
 //               ex1 -pa -d occa-omp
 //               ex1 -pa -d ceed-cpu
+//               ex1 -pa -d ceed-cpu -o 4 -a
 //             * ex1 -pa -d ceed-cuda
 //             * ex1 -pa -d ceed-hip
 //               ex1 -pa -d ceed-cuda:/gpu/cuda/shared
@@ -73,6 +74,7 @@ int main(int argc, char *argv[])
    bool pa = false;
    const char *device_config = "cpu";
    bool visualization = true;
+   bool algebraic_ceed = false;
 
    OptionsParser args(argc, argv);
    args.AddOption(&mesh_file, "-m", "--mesh",
@@ -86,6 +88,8 @@ int main(int argc, char *argv[])
                   "--no-partial-assembly", "Enable Partial Assembly.");
    args.AddOption(&device_config, "-d", "--device",
                   "Device configuration string, see Device::Configure().");
+   args.AddOption(&algebraic_ceed, "-a", "--algebraic", "-no-a", "--no-algebraic",
+                  "Use algebraic Ceed solver");
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                   "--no-visualization",
                   "Enable or disable GLVis visualization.");
@@ -207,17 +211,25 @@ int main(int argc, char *argv[])
       umf_solver.Mult(B, X);
 #endif
    }
-   else // Jacobi preconditioning in partial assembly mode
+   else
    {
-      if (UsesTensorBasis(fespace))
+#ifdef MFEM_USE_CEED
+      if (DeviceCanUseCeed() && algebraic_ceed)
       {
-         OperatorJacobiSmoother M(a, ess_tdof_list);
+         AlgebraicCeedSolver M(a, ess_tdof_list);
          PCG(*A, M, B, X, 1, 400, 1e-12, 0.0);
       }
       else
-      {
-         CG(*A, B, X, 1, 400, 1e-12, 0.0);
-      }
+#endif
+         if (UsesTensorBasis(fespace))
+         {
+            OperatorJacobiSmoother M(a, ess_tdof_list);
+            PCG(*A, M, B, X, 1, 400, 1e-12, 0.0);
+         }
+         else
+         {
+            CG(*A, B, X, 1, 400, 1e-12, 0.0);
+         }
    }
 
    // 12. Recover the solution as a finite element grid function.

examples/ex1p.cpp

@@ -32,6 +32,7 @@
 //               mpirun -np 4 ex1p -pa -d occa-cuda
 //               mpirun -np 4 ex1p -pa -d raja-omp
 //               mpirun -np 4 ex1p -pa -d ceed-cpu
+//               mpirun -np 4 ex1p -pa -d ceed-cpu -o 4 -a
 //             * mpirun -np 4 ex1p -pa -d ceed-cuda
 //             * mpirun -np 4 ex1p -pa -d ceed-hip
 //               mpirun -np 4 ex1p -pa -d ceed-cuda:/gpu/cuda/shared
@@ -53,6 +54,7 @@
 //               optional connection to the GLVis tool for visualization.
 
 #include "mfem.hpp"
+
 #include <fstream>
 #include <iostream>
 
@@ -62,10 +64,9 @@ using namespace mfem;
 int main(int argc, char *argv[])
 {
    // 1. Initialize MPI.
-   int num_procs, myid;
-   MPI_Init(&argc, &argv);
-   MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
-   MPI_Comm_rank(MPI_COMM_WORLD, &myid);
+   MPI_Session mpi;
+   int num_procs = mpi.WorldSize();
+   int myid = mpi.WorldRank();
 
    // 2. Parse command-line options.
    const char *mesh_file = "../data/star.mesh";
@@ -74,6 +75,7 @@ int main(int argc, char *argv[])
    bool pa = false;
    const char *device_config = "cpu";
    bool visualization = true;
+   bool algebraic_ceed = false;
 
    OptionsParser args(argc, argv);
    args.AddOption(&mesh_file, "-m", "--mesh",
@@ -87,6 +89,8 @@ int main(int argc, char *argv[])
                   "--no-partial-assembly", "Enable Partial Assembly.");
    args.AddOption(&device_config, "-d", "--device",
                   "Device configuration string, see Device::Configure().");
+   args.AddOption(&algebraic_ceed, "-a", "--algebraic", "-no-a", "--no-algebraic",
+                  "Use algebraic Ceed solver");
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                   "--no-visualization",
                   "Enable or disable GLVis visualization.");
@@ -204,7 +208,9 @@ int main(int argc, char *argv[])
    //     domain integrator.
    ParBilinearForm a(&fespace);
    if (pa) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
+
    a.AddDomainIntegrator(new DiffusionIntegrator(one));
+   // a.AddDomainIntegrator(new MassIntegrator(one));
 
    // 12. Assemble the parallel bilinear form and the corresponding linear
    //     system, applying any necessary transformations such as: parallel
@@ -225,13 +231,23 @@ int main(int argc, char *argv[])
    {
       if (UsesTensorBasis(fespace))
       {
-         prec = new OperatorJacobiSmoother(a, ess_tdof_list);
+#ifdef MFEM_USE_CEED
+         if (DeviceCanUseCeed() && algebraic_ceed)
+         {
+            prec = new AlgebraicCeedSolver(a, ess_tdof_list);
+         }
+         else
+#endif
+         {
+            prec = new OperatorJacobiSmoother(a, ess_tdof_list);
+         }
       }
    }
    else
    {
       prec = new HypreBoomerAMG;
    }
+
    CGSolver cg(MPI_COMM_WORLD);
    cg.SetRelTol(1e-12);
    cg.SetMaxIter(2000);
@@ -277,7 +293,6 @@ int main(int argc, char *argv[])
    {
       delete fec;
    }
-   MPI_Finalize();
 
    return 0;
 }

fem/CMakeLists.txt

@@ -34,6 +34,7 @@ set(SRCS
   bilininteg_vecdiffusion_mf.cpp
   bilininteg_vecmass.cpp
   bilininteg_vecmass_mf.cpp
+  ceed_algebraic.cpp
   coefficient.cpp
   complex_fem.cpp
   convergence.cpp
@@ -50,6 +51,12 @@ set(SRCS
   libceed/ceed.cpp
   libceed/diffusion.cpp
   libceed/mass.cpp
+  libceed/ceed-assemble.cpp
+  libceed/ceedsolvers-atpmg.cpp
+  libceed/ceedsolvers-interpolation.cpp
+  libceed/ceedsolvers-utility.cpp
+  libceed/ceed-wrappers.cpp
+  libceed/ceedsolvers-qcoarsen.cpp
   linearform.cpp
   lininteg.cpp
   multigrid.cpp
@@ -72,6 +79,7 @@ set(HDRS
   bilinearform.hpp
   bilinearform_ext.hpp
   bilininteg.hpp
+  ceed_algebraic.hpp
   coefficient.hpp
   complex_fem.hpp
   convergence.hpp
@@ -89,6 +97,12 @@ set(HDRS
   libceed/ceed.hpp
   libceed/diffusion.hpp
   libceed/mass.hpp
+  libceed/ceed-assemble.hpp
+  libceed/ceedsolvers-atpmg.h
+  libceed/ceedsolvers-interpolation.h
+  libceed/ceedsolvers-utility.h
+  libceed/ceed-wrappers.hpp
+  libceed/ceedsolvers-qcoarsen.h
   linearform.hpp
   lininteg.hpp
   multigrid.hpp

fem/bilininteg.hpp

@@ -2049,6 +2049,8 @@ public:
 
    static const IntegrationRule &GetRule(const FiniteElement &trial_fe,
                                          const FiniteElement &test_fe);
+
+   CeedData * GetCeedData() { return ceedDataPtr; }
 };
 
 /** Class for local mass matrix assembling a(u,v) := (Q u, v) */
@@ -2115,6 +2117,8 @@ public:
    static const IntegrationRule &GetRule(const FiniteElement &trial_fe,
                                          const FiniteElement &test_fe,
                                          ElementTransformation &Trans);
+
+   CeedData * GetCeedData() { return ceedDataPtr; }
 };
 
 /** Mass integrator (u, v) restricted to the boundary of a domain */
@@ -2271,6 +2275,10 @@ public:
    virtual void AssembleDiagonalMF(Vector &diag);
    virtual void AddMultPA(const Vector &x, Vector &y) const;
    virtual void AddMultMF(const Vector &x, Vector &y) const;
+
+#ifdef MFEM_USE_CEED
+   CeedData * GetCeedData() { return ceedDataPtr; }
+#endif
 };
 
 
@@ -2728,6 +2736,10 @@ public:
    virtual void AssembleDiagonalMF(Vector &diag);
    virtual void AddMultPA(const Vector &x, Vector &y) const;
    virtual void AddMultMF(const Vector &x, Vector &y) const;
+
+#ifdef MFEM_USE_CEED
+   CeedData * GetCeedData() { return ceedDataPtr; }
+#endif
 };
 
 /** Integrator for the linear elasticity form:

fem/ceed_algebraic.cpp

@@ -0,0 +1,815 @@
+// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "ceed_algebraic.hpp"
+
+#ifdef MFEM_USE_CEED
+#include "../fem/bilinearform.hpp"
+#include "../fem/fespace.hpp"
+#include "../fem/libceed/ceedsolvers-atpmg.h"
+#include "../fem/libceed/ceedsolvers-interpolation.h"
+#include "../fem/libceed/ceed-assemble.hpp"
+#include "../fem/libceed/ceedsolvers-qcoarsen.h"
+#include "../fem/libceed/ceedsolvers-sparsify.h"
+#include "../fem/pfespace.hpp"
+
+namespace mfem
+{
+
+Solver *BuildSmootherFromCeed(MFEMCeedOperator &op, bool chebyshev)
+{
+   CeedOperator ceed_op = op.GetCeedOperator();
+   const Array<int> &ess_tdofs = op.GetEssentialTrueDofs();
+   const Operator *P = op.GetProlongation();
+   // Assemble the a local diagonal, in the sense of L-vector
+   CeedVector diagceed;
+   CeedInt length;
+   CeedOperatorGetSize(ceed_op, &length);
+   CeedVectorCreate(internal::ceed, length, &diagceed);
+   CeedMemType mem;
+   CeedGetPreferredMemType(internal::ceed, &mem);
+   if (!Device::Allows(Backend::CUDA) || mem != CEED_MEM_DEVICE)
+   {
+      mem = CEED_MEM_HOST;
+   }
+   Vector local_diag(length);
+   CeedScalar *ptr = (mem == CEED_MEM_HOST) ? local_diag.HostWrite() :
+                     local_diag.Write(true);
+   CeedVectorSetArray(diagceed, mem, CEED_USE_POINTER, ptr);
+   CeedOperatorLinearAssembleDiagonal(ceed_op, diagceed, CEED_REQUEST_IMMEDIATE);
+   CeedVectorTakeArray(diagceed, mem, NULL);
+
+   Vector t_diag;
+   if (P)
+   {
+      t_diag.SetSize(P->Width());
+      P->MultTranspose(local_diag, t_diag);
+   }
+   else
+   {
+      t_diag.NewMemoryAndSize(local_diag.GetMemory(), length, false);
+   }
+   Solver *out = NULL;
+   if (chebyshev)
+   {
+      const int cheb_order = 3;
+      out = new OperatorChebyshevSmoother(&op, t_diag, ess_tdofs, cheb_order);
+   }
+   else
+   {
+      const double jacobi_scale = 0.65;
+      out = new OperatorJacobiSmoother(t_diag, ess_tdofs, jacobi_scale);
+   }
+   CeedVectorDestroy(&diagceed);
+   return out;
+}
+
+#ifdef MFEM_USE_MPI
+
+class CeedAMG : public Solver
+{
+public:
+   CeedAMG(MFEMCeedOperator &oper, HypreParMatrix *P, bool amgx=false,
+           const std::string amgx_config_file="")
+   {
+      MFEM_ASSERT(P != NULL, "");
+      const Array<int> ess_tdofs = oper.GetEssentialTrueDofs();
+      height = width = oper.Height();
+
+      CeedOperatorFullAssemble(oper.GetCeedOperator(), &mat_local);
+
+      {
+         HypreParMatrix hypre_local(
+            P->GetComm(), P->GetGlobalNumRows(), P->RowPart(), mat_local);
+         op_assembled = RAP(&hypre_local, P);
+      }
+      HypreParMatrix *mat_e = op_assembled->EliminateRowsCols(ess_tdofs);
+      delete mat_e;
+
+#ifdef MFEM_USE_AMGX
+      if (amgx)
+      {
+         if (amgx_config_file == "")
+         {
+            bool amgx_verbose = false;
+            amg = new AmgXSolver(op_assembled->GetComm(),
+                                 AmgXSolver::PRECONDITIONER, amgx_verbose);
+         }
+         else
+         {
+            AmgXSolver * amgx_prec = new AmgXSolver;
+            amgx_prec->ReadParameters(amgx_config_file, AmgXSolver::EXTERNAL);
+            amgx_prec->InitExclusiveGPU(MPI_COMM_WORLD);
+            amgx_prec->SetOperator(*op_assembled);
+            amg = amgx_prec;
+         }
+      }
+      else
+#endif
+      {
+         HypreBoomerAMG * hypre_amg = new HypreBoomerAMG(*op_assembled);
+         hypre_amg->SetPrintLevel(0);
+         amg = hypre_amg;
+      }
+   }
+
+   void SetOperator(const Operator &op) { amg->SetOperator(op); }
+   void Mult(const Vector &x, Vector &y) const { amg->Mult(x, y); }
+
+   ~CeedAMG()
+   {
+      delete op_assembled;
+      delete amg;
+      delete mat_local;
+   }
+
+private:
+   SparseMatrix *mat_local;
+   HypreParMatrix *op_assembled;
+   Solver *amg;
+};
+
+/**
+   Too much copied code; this and CeedAMG should probably inherit
+   from a common base class.
+*/
+class CeedSparsifyAMG : public Solver
+{
+public:
+   CeedSparsifyAMG(MFEMCeedOperator &oper, HypreParMatrix *P, bool amgx=false)
+   {
+      MFEM_ASSERT(P != NULL, "");
+      const Array<int> ess_tdofs = oper.GetEssentialTrueDofs();
+      height = width = oper.Height();
+
+      CeedSparsifySimple(oper.GetCeedOperator(), &sparse_basis, &sparse_oper);
+      CeedOperatorFullAssemble(sparse_oper, &mat_local);
+
+      {
+         HypreParMatrix hypre_local(
+            P->GetComm(), P->GetGlobalNumRows(), P->RowPart(), mat_local);
+         op_assembled = RAP(&hypre_local, P);
+      }
+      HypreParMatrix *mat_e = op_assembled->EliminateRowsCols(ess_tdofs);
+      delete mat_e;
+
+#ifdef MFEM_USE_AMGX
+      if (amgx)
+      {
+         bool amgx_verbose = false;
+         amg = new AmgXSolver(op_assembled->GetComm(),
+                              AmgXSolver::PRECONDITIONER, amgx_verbose);
+         amg->SetOperator(*op_assembled);
+      }
+      else
+#endif
+      {
+         HypreBoomerAMG * hypre_amg = new HypreBoomerAMG(*op_assembled);
+         hypre_amg->SetPrintLevel(0);
+         amg = hypre_amg;
+      }
+   }
+
+   void SetOperator(const Operator &op) { amg->SetOperator(op); }
+   void Mult(const Vector &x, Vector &y) const { amg->Mult(x, y); }
+
+   ~CeedSparsifyAMG()
+   {
+      CeedBasisDestroy(&sparse_basis);
+      CeedOperatorDestroy(&sparse_oper);
+
+      delete op_assembled;
+      delete amg;
+      delete mat_local;
+   }
+
+private:
+   CeedOperator sparse_oper;
+   CeedBasis sparse_basis;
+   SparseMatrix *mat_local;
+   HypreParMatrix *op_assembled;
+   Solver *amg;
+};
+
+#endif
+
+void CoarsenEssentialDofs(const Operator &interp,
+                          const Array<int> &ho_ess_tdofs,
+                          Array<int> &alg_lo_ess_tdofs)
+{
+   Vector ho_boundary_ones(interp.Height());
+   ho_boundary_ones = 0.0;
+   const int *ho_ess_tdofs_h = ho_ess_tdofs.HostRead();
+   for (int i=0; i<ho_ess_tdofs.Size(); ++i)
+   {
+      ho_boundary_ones[ho_ess_tdofs_h[i]] = 1.0;
+   }
+   Vector lo_boundary_ones(interp.Width());
+   interp.MultTranspose(ho_boundary_ones, lo_boundary_ones);
+   auto lobo = lo_boundary_ones.HostRead();
+   for (int i = 0; i < lo_boundary_ones.Size(); ++i)
+   {
+      if (lobo[i] > 0.9)
+      {
+         alg_lo_ess_tdofs.Append(i);
+      }
+   }
+}
+
+template <typename INTEG>
+int TryToAddCeedSubOperator(BilinearFormIntegrator *integ_in, CeedOperator op)
+{
+   INTEG *integ = dynamic_cast<INTEG*>(integ_in);
+   if (integ != NULL)
+   {
+      CeedCompositeOperatorAddSub(op, integ->GetCeedData()->oper);
+      return 1;
+   }
+   return 0;
+}
+
+CeedOperator CreateCeedCompositeOperatorFromBilinearForm(BilinearForm &form)
+{
+   CeedOperator op;
+   CeedCompositeOperatorCreate(internal::ceed, &op);
+
+   // Get the domain bilinear form integrators (DBFIs)
+   Array<BilinearFormIntegrator*> *bffis = form.GetDBFI();
+   int num_integrators = bffis->Size();
+
+   int count = 0;
+   for (int i = 0; i < num_integrators; ++i)
+   {
+      BilinearFormIntegrator *integ = (*bffis)[i];
+      count += TryToAddCeedSubOperator<DiffusionIntegrator>(integ, op);
+      count += TryToAddCeedSubOperator<MassIntegrator>(integ, op);
+      count += TryToAddCeedSubOperator<VectorDiffusionIntegrator>(integ, op);
+      count += TryToAddCeedSubOperator<VectorMassIntegrator>(integ, op);
+   }
+   if (count != num_integrators)
+   {
+      mfem_error("Some integrator does not support Ceed!");
+   }
+   return op;
+}
+
+CeedOperator CoarsenCeedCompositeOperator(
+   CeedOperator op, CeedElemRestriction er,
+   CeedBasis c2f, int order_reduction,
+   int qorder_reduction,
+   CeedQuadMode fine_qmode, CeedQuadMode coarse_qmode
+)
+{
+   bool isComposite;
+   CeedOperatorIsComposite(op, &isComposite);
+   MFEM_ASSERT(isComposite, "");
+
+   CeedOperator op_coarse;
+   CeedCompositeOperatorCreate(internal::ceed, &op_coarse);
+
+   int nsub;
+   CeedOperatorGetNumSub(op, &nsub);
+   CeedOperator *subops;
+   CeedOperatorGetSubList(op, &subops);
+   for (int isub=0; isub<nsub; ++isub)
+   {
+      CeedOperator subop = subops[isub];
+      CeedBasis basis_coarse, basis_c2f;
+      CeedOperator subop_coarse, t_subop_coarse;
+      CeedATPMGOperator(subop, order_reduction, er, &basis_coarse, &basis_c2f,
+                        &t_subop_coarse);
+      if (qorder_reduction == 0)
+      {
+         subop_coarse = t_subop_coarse;
+      }
+      else
+      {
+         CeedVector qcoarsen_assembledqf;
+         CeedQFunctionContext qcoarsen_context;
+         CeedOperatorQCoarsen(t_subop_coarse, qorder_reduction, &subop_coarse,
+                              &qcoarsen_assembledqf, &qcoarsen_context,
+                              fine_qmode, coarse_qmode);
+         CeedVectorDestroy(&qcoarsen_assembledqf); // todo: delete inside previous function?
+         CeedQFunctionContextDestroy(&qcoarsen_context);
+         CeedOperatorDestroy(&t_subop_coarse);
+      }
+      CeedBasisDestroy(&basis_coarse); // refcounted by subop_coarse
+      CeedBasisDestroy(&basis_c2f);
+      CeedCompositeOperatorAddSub(op_coarse, subop_coarse);
+      CeedOperatorDestroy(&subop_coarse); // refcounted by composite operator
+   }
+   return op_coarse;
+}
+
+AlgebraicCeedMultigrid::AlgebraicCeedMultigrid(
+   AlgebraicSpaceHierarchy &hierarchy,
+   BilinearForm &form,
+   const Array<int> &ess_tdofs,
+   int print_level,
+   double contrast_threshold,
+   int switch_amg_order,
+   bool collocate_coarse,
+   bool sparsification,
+   const std::string amgx_config_file
+) : GeometricMultigrid(hierarchy)
+{
+   // Construct finest level
+   ceed_operators.Prepend(CreateCeedCompositeOperatorFromBilinearForm(form));
+   essentialTrueDofs.Prepend(new Array<int>);
+   *essentialTrueDofs[0] = ess_tdofs;
+
+   int current_order = hierarchy.GetFESpaceAtLevel(0).GetOrder(0);
+
+   // Construct interpolation, operators, at all levels of hierarchy by coarsening
+   int level_counter = 0;
+   while (current_order > 1)
+   {
+      double minq, maxq, absmin;
+      CeedOperatorGetHeuristics(ceed_operators[0], &minq, &maxq, &absmin);
+      // TODO: in principle we need to communicate heuristics across
+      // processors!
+      double heuristic = std::max(std::abs(minq), std::abs(maxq)) / absmin;
+
+      int order_reduction;
+      if (heuristic > contrast_threshold && current_order <= switch_amg_order)
+      {
+         // assemble at this level
+         break;
+      }
+      else if (heuristic > contrast_threshold || current_order == 3)
+      {
+         // coarsening directly from 3 to 1 appears to be bad
+         order_reduction = 1;
+      }
+      else
+      {
+         order_reduction = current_order - (current_order/2);
+      }
+      if (print_level > 0)
+      {
+         std::cout << "  lc: " << level_counter << " heuristic = " << heuristic
+                   << ", coarsening from order " << current_order 
+                   << " to " << current_order - order_reduction << std::endl;
+      }
+      hierarchy.PrependPCoarsenedLevel(current_order, order_reduction);
+      current_order = current_order - order_reduction;
+
+      AlgebraicCoarseSpace &space = hierarchy.GetAlgebraicCoarseSpace(0);
+      // int qor = (level_counter == 0) ? order_reduction + 1 : order_reduction;
+      int qor = order_reduction;
+      if (collocate_coarse)
+      {
+         if (level_counter == 0)
+         {
+            ceed_operators.Prepend(
+               CoarsenCeedCompositeOperator(
+                  ceed_operators[0], space.GetCeedElemRestriction(),
+                  space.GetCeedCoarseToFine(), space.GetOrderReduction(),
+                  qor, CEED_GAUSS, CEED_GAUSS_LOBATTO)
+               );
+         }
+         else
+         {
+            ceed_operators.Prepend(
+               CoarsenCeedCompositeOperator(
+                  ceed_operators[0], space.GetCeedElemRestriction(),
+                  space.GetCeedCoarseToFine(), space.GetOrderReduction(),
+                  qor, CEED_GAUSS_LOBATTO, CEED_GAUSS_LOBATTO)
+               );
+         }
+      }
+      else
+      {
+         ceed_operators.Prepend(
+            CoarsenCeedCompositeOperator(
+               ceed_operators[0], space.GetCeedElemRestriction(),
+               space.GetCeedCoarseToFine(), space.GetOrderReduction(),
+               qor, CEED_GAUSS, CEED_GAUSS)
+            );
+      }
+
+      Operator *P = hierarchy.GetProlongationAtLevel(0);
+      essentialTrueDofs.Prepend(new Array<int>);
+      CoarsenEssentialDofs(*P, *essentialTrueDofs[1],
+                           *essentialTrueDofs[0]);
+
+      level_counter++;
+   }
+
+   int nlevels = fespaces.GetNumLevels();
+   // Add the operators and smoothers to the hierarchy, from coarse to fine
+   for (int ilevel=0; ilevel<nlevels; ++ilevel)
+   {
+      FiniteElementSpace &space = hierarchy.GetFESpaceAtLevel(ilevel);
+      const Operator *P = space.GetProlongationMatrix();
+      MFEMCeedOperator *op = new MFEMCeedOperator(
+         ceed_operators[ilevel], *essentialTrueDofs[ilevel], P);
+      Solver *smoother;
+      if (ilevel != 0)
+      {
+         smoother = BuildSmootherFromCeed(*op, true);
+      }
+      else
+      {
+         bool assemble_matrix = false;
+#ifdef MFEM_USE_MPI
+#ifdef MFEM_USE_AMGX
+	 assemble_matrix = true;
+#else
+         if (!Device::Allows(Backend::CUDA)) { assemble_matrix = true; }
+#endif
+         HypreParMatrix *P_mat = NULL;
+         if (assemble_matrix)
+         {
+            if (nlevels == 1)
+            {
+               // Only one level -- no coarsening, finest level
+               ParFiniteElementSpace *pfes
+                  = dynamic_cast<ParFiniteElementSpace*>(&space);
+               if (pfes) { P_mat = pfes->Dof_TrueDof_Matrix(); }
+            }
+            else
+            {
+               ParAlgebraicCoarseSpace *pspace
+                  = dynamic_cast<ParAlgebraicCoarseSpace*>(&space);
+               if (pspace) { P_mat = pspace->GetProlongationHypreParMatrix(); }
+            }
+         }
+         if (P_mat)
+         {
+            if (current_order > 1 && sparsification)
+            {
+               if (print_level >= 1)
+               {
+                  std::cout << "      sparsify AMG." << std::endl;
+               }
+               smoother = new CeedSparsifyAMG(*op, P_mat, Device::Allows(Backend::CUDA));
+            }
+            else
+            {
+               if (print_level >= 1)
+               {
+                  std::cout << "      no-sparsify AMG." << std::endl;
+               }
+               smoother = new CeedAMG(*op, P_mat, Device::Allows(Backend::CUDA),
+                                      amgx_config_file);
+            }
+         }
+         else
+#endif
+         {
+            smoother = BuildSmootherFromCeed(*op, true);
+         }
+      }
+      AddLevel(op, smoother, true, true);
+   }
+}
+
+AlgebraicCeedMultigrid::~AlgebraicCeedMultigrid()
+{
+   for (int i=0; i<ceed_operators.Size(); ++i)
+   {
+      CeedOperatorDestroy(&ceed_operators[i]);
+   }
+}
+
+void AlgebraicSpaceHierarchy::AddCoarseLevel(AlgebraicCoarseSpace* space,
+                                             CeedElemRestriction er)
+{
+   MFEM_VERIFY(meshes.Size() >= 1, "At least one level must exist!");
+   Mesh* finemesh = meshes[0];
+   meshes.Prepend(finemesh); // every entry of meshes points to finest mesh
+   ownedMeshes.Prepend(false);
+   fespaces.Prepend(space);
+   ownedFES.Prepend(true); // owns all but finest
+   ceed_interpolations.Prepend(new MFEMCeedInterpolation(
+                                  internal::ceed,
+                                  space->GetCeedCoarseToFine(),
+                                  space->GetCeedElemRestriction(),
+                                  er)
+      );
+   const SparseMatrix *R = fespaces[1]->GetRestrictionMatrix();
+   if (R)
+   {
+      R->BuildTranspose();
+      R_tr.Prepend(new TransposeOperator(*R));
+   }
+   else
+   {
+      R_tr.Prepend(NULL);
+   }
+   prolongations.Prepend(ceed_interpolations[0]->SetupRAP(
+                            space->GetProlongationMatrix(), R_tr[0]));
+   ownedProlongations.Prepend(prolongations[0] != ceed_interpolations[0]);   
+}
+
+// the ifdefs and dynamic casts are very ugly, but the interface is kinda nice?
+void AlgebraicSpaceHierarchy::PrependPCoarsenedLevel(
+   int current_order, int order_reduction)
+{
+   MFEM_VERIFY(fespaces.Size() >= 1, "At least one level must exist!");
+   int dim = meshes[0]->Dimension();
+
+   AlgebraicCoarseSpace *fine_alg_space =
+      dynamic_cast<AlgebraicCoarseSpace*>(fespaces[0]);
+   CeedElemRestriction current_er;
+#ifdef MFEM_USE_MPI
+   GroupCommunicator *gc = NULL;
+#endif
+   if (fine_alg_space)
+   {
+      current_er = fine_alg_space->GetCeedElemRestriction();
+#ifdef MFEM_USE_MPI
+      ParAlgebraicCoarseSpace *par_alg_space =
+         dynamic_cast<ParAlgebraicCoarseSpace*>(fine_alg_space);
+      if (par_alg_space) { gc = par_alg_space->GetGroupCommunicator(); }
+#endif
+   }
+   else
+   {
+      current_er = fine_er;
+#ifdef MFEM_USE_MPI
+      ParFiniteElementSpace *pfes =
+         dynamic_cast<ParFiniteElementSpace*>(fespaces[0]);
+      if (pfes) { gc = &pfes->GroupComm(); }
+#endif
+   }
+   AlgebraicCoarseSpace *space;
+#ifdef MFEM_USE_MPI
+   if (gc)
+   {
+      space = new ParAlgebraicCoarseSpace(
+         *fespaces[0], current_er, current_order, dim, order_reduction, gc);
+   }
+   else
+#endif
+   {
+      space = new AlgebraicCoarseSpace(
+         *fespaces[0], current_er, current_order, dim, order_reduction);
+   }
+   AddCoarseLevel(space, current_er);
+}
+
+AlgebraicSpaceHierarchy::AlgebraicSpaceHierarchy(FiniteElementSpace &fes)
+{
+   int order = fes.GetOrder(0);
+
+   meshes.Prepend(fes.GetMesh());
+   ownedMeshes.Prepend(false);
+   fespaces.Prepend(&fes);
+   ownedFES.Prepend(false);
+
+   Ceed ceed = internal::ceed;
+   InitCeedTensorRestriction(fes, ceed, &fine_er);
+}
+
+AlgebraicCoarseSpace::AlgebraicCoarseSpace(
+   FiniteElementSpace &fine_fes,
+   CeedElemRestriction fine_er,
+   int order,
+   int dim,
+   int order_reduction_
+) : order_reduction(order_reduction_)
+{
+   order_reduction = order_reduction_;
+
+   CeedATPMGElemRestriction(order, order_reduction, fine_er,
+                            &ceed_elem_restriction, dof_map );
+   CeedBasisATPMGCoarseToFine(internal::ceed, order+1, dim,
+                              order_reduction, &coarse_to_fine );
+   CeedElemRestrictionGetLVectorSize(ceed_elem_restriction, &ndofs);
+   mesh = fine_fes.GetMesh();
+}
+
+AlgebraicCoarseSpace::~AlgebraicCoarseSpace()
+{
+   free(dof_map);
+   CeedBasisDestroy(&coarse_to_fine);
+   CeedElemRestrictionDestroy(&ceed_elem_restriction);
+}
+
+#ifdef MFEM_USE_MPI
+
+ParAlgebraicCoarseSpace::ParAlgebraicCoarseSpace(
+   FiniteElementSpace &fine_fes,
+   CeedElemRestriction fine_er,
+   int order,
+   int dim,
+   int order_reduction_,
+   GroupCommunicator *gc_fine)
+   : AlgebraicCoarseSpace(fine_fes, fine_er, order, dim, order_reduction_)
+{
+   int lsize;
+   CeedElemRestrictionGetLVectorSize(ceed_elem_restriction, &lsize);
+   const Table &group_ldof_fine = gc_fine->GroupLDofTable();
+
+   ldof_group.SetSize(lsize);
+   ldof_group = 0;
+
+   GroupTopology &group_topo = gc_fine->GetGroupTopology();
+   gc = new GroupCommunicator(group_topo);
+   Table &group_ldof = gc->GroupLDofTable();
+   group_ldof.MakeI(group_ldof_fine.Size());
+   for (int g=1; g<group_ldof_fine.Size(); ++g)
+   {
+      int nldof_fine_g = group_ldof_fine.RowSize(g);
+      const int *ldof_fine_g = group_ldof_fine.GetRow(g);
+      for (int i=0; i<nldof_fine_g; ++i)
+      {
+         int icoarse = dof_map[ldof_fine_g[i]];
+         if (icoarse >= 0)
+         {
+            group_ldof.AddAColumnInRow(g);
+            ldof_group[icoarse] = g;
+         }
+      }
+   }
+   group_ldof.MakeJ();
+   for (int g=1; g<group_ldof_fine.Size(); ++g)
+   {
+      int nldof_fine_g = group_ldof_fine.RowSize(g);
+      const int *ldof_fine_g = group_ldof_fine.GetRow(g);
+      for (int i=0; i<nldof_fine_g; ++i)
+      {
+         int icoarse = dof_map[ldof_fine_g[i]];
+         if (icoarse >= 0)
+         {
+            group_ldof.AddConnection(g, icoarse);
+         }
+      }
+   }
+   group_ldof.ShiftUpI();
+   gc->Finalize();
+   ldof_ltdof.SetSize(lsize);
+   ldof_ltdof = -2;
+   int ltsize = 0;
+   for (int i=0; i<lsize; ++i)
+   {
+      int g = ldof_group[i];
+      if (group_topo.IAmMaster(g))
+      {
+         ldof_ltdof[i] = ltsize;
+         ++ltsize;
+      }
+   }
+   gc->SetLTDofTable(ldof_ltdof);
+   gc->Bcast(ldof_ltdof);
+
+   R_mat = new SparseMatrix(ltsize, lsize);
+   for (int j=0; j<lsize; ++j)
+   {
+      if (group_topo.IAmMaster(ldof_group[j]))
+      {
+         int i = ldof_ltdof[j];
+         R_mat->Set(i,j,1.0);
+      }
+   }
+   R_mat->Finalize();
+
+   if (Device::Allows(Backend::DEVICE_MASK))
+   {
+      P = new DeviceConformingProlongationOperator(*gc, R_mat);
+   }
+   else
+   {
+      P = new ConformingProlongationOperator(lsize, *gc);
+   }
+   P_mat = NULL;
+}
+
+HypreParMatrix *ParAlgebraicCoarseSpace::GetProlongationHypreParMatrix()
+{
+   if (P_mat) { return P_mat; }
+
+   ParMesh *pmesh = dynamic_cast<ParMesh*>(mesh);
+   MFEM_VERIFY(pmesh != NULL, "");
+   Array<HYPRE_Int> dof_offsets, tdof_offsets, tdof_nb_offsets;
+   Array<HYPRE_Int> *offsets[2] = {&dof_offsets, &tdof_offsets};
+   int lsize = P->Height();
+   int ltsize = P->Width();
+   HYPRE_Int loc_sizes[2] = {lsize, ltsize};
+   pmesh->GenerateOffsets(2, loc_sizes, offsets);
+
+   MPI_Comm comm = pmesh->GetComm();
+
+   const GroupTopology &group_topo = gc->GetGroupTopology();
+
+   if (HYPRE_AssumedPartitionCheck())
+   {
+      // communicate the neighbor offsets in tdof_nb_offsets
+      int nsize = group_topo.GetNumNeighbors()-1;
+      MPI_Request *requests = new MPI_Request[2*nsize];
+      MPI_Status  *statuses = new MPI_Status[2*nsize];
+      tdof_nb_offsets.SetSize(nsize+1);
+      tdof_nb_offsets[0] = tdof_offsets[0];
+
+      // send and receive neighbors' local tdof offsets
+      int request_counter = 0;
+      for (int i = 1; i <= nsize; i++)
+      {
+         MPI_Irecv(&tdof_nb_offsets[i], 1, HYPRE_MPI_INT,
+                   group_topo.GetNeighborRank(i), 5365, comm,
+                   &requests[request_counter++]);
+      }
+      for (int i = 1; i <= nsize; i++)
+      {
+         MPI_Isend(&tdof_nb_offsets[0], 1, HYPRE_MPI_INT,
+                   group_topo.GetNeighborRank(i), 5365, comm,
+                   &requests[request_counter++]);
+      }
+      MPI_Waitall(request_counter, requests, statuses);
+
+      delete [] statuses;
+      delete [] requests;
+   }
+
+   HYPRE_Int *i_diag = Memory<HYPRE_Int>(lsize+1);
+   HYPRE_Int *j_diag = Memory<HYPRE_Int>(ltsize);
+   int diag_counter;
+
+   HYPRE_Int *i_offd = Memory<HYPRE_Int>(lsize+1);
+   HYPRE_Int *j_offd = Memory<HYPRE_Int>(lsize-ltsize);
+   int offd_counter;
+
+   HYPRE_Int *cmap   = Memory<HYPRE_Int>(lsize-ltsize);
+
+   HYPRE_Int *col_starts = tdof_offsets;
+   HYPRE_Int *row_starts = dof_offsets;
+
+   Array<Pair<HYPRE_Int, int> > cmap_j_offd(lsize-ltsize);
+
+   i_diag[0] = i_offd[0] = 0;
+   diag_counter = offd_counter = 0;
+   for (int i_ldof = 0; i_ldof < lsize; i_ldof++)
+   {
+      int g = ldof_group[i_ldof];
+      int i_ltdof = ldof_ltdof[i_ldof];
+      if (group_topo.IAmMaster(g))
+      {
+         j_diag[diag_counter++] = i_ltdof;
+      }
+      else
+      {
+         HYPRE_Int global_tdof_number;
+         int g = ldof_group[i_ldof];
+         if (HYPRE_AssumedPartitionCheck())
+         {
+            global_tdof_number
+               = i_ltdof + tdof_nb_offsets[group_topo.GetGroupMaster(g)];
+         }
+         else
+         {
+            global_tdof_number
+               = i_ltdof + tdof_offsets[group_topo.GetGroupMasterRank(g)];
+         }
+
+         cmap_j_offd[offd_counter].one = global_tdof_number;
+         cmap_j_offd[offd_counter].two = offd_counter;
+         offd_counter++;
+      }
+      i_diag[i_ldof+1] = diag_counter;
+      i_offd[i_ldof+1] = offd_counter;
+   }
+
+   SortPairs<HYPRE_Int, int>(cmap_j_offd, offd_counter);
+
+   for (int i = 0; i < offd_counter; i++)
+   {
+      cmap[i] = cmap_j_offd[i].one;
+      j_offd[cmap_j_offd[i].two] = i;
+   }
+
+   P_mat = new HypreParMatrix(
+      comm, pmesh->GetMyRank(), pmesh->GetNRanks(),
+      row_starts, col_starts,
+      i_diag, j_diag, i_offd, j_offd,
+      cmap, offd_counter
+   );
+
+   P_mat->CopyRowStarts();
+   P_mat->CopyColStarts();
+
+   return P_mat;
+}
+
+ParAlgebraicCoarseSpace::~ParAlgebraicCoarseSpace()
+{
+   delete P;
+   delete R_mat;
+   delete P_mat;
+   delete gc;
+}
+
+#endif
+
+
+} // namespace mfem
+#endif // MFEM_USE_CEED

fem/ceed_algebraic.hpp

@@ -0,0 +1,198 @@
+// Copyright (c) 2010-2020, 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_CEED_ALGEBRAIC_HPP
+#define MFEM_CEED_ALGEBRAIC_HPP
+
+#include "../config/config.hpp"
+
+#ifdef MFEM_USE_CEED
+#include "fespacehierarchy.hpp"
+#include "multigrid.hpp"
+#include "libceed/ceedsolvers-utility.h"
+#include "libceed/ceed-wrappers.hpp"
+
+namespace mfem
+{
+
+/** @brief A way to use algebraic levels in a Multigrid object
+
+    This is analogous to a FiniteElementSpace but with no Mesh information,
+    constructed in a semi-algebraic way. */
+class AlgebraicCoarseSpace : public FiniteElementSpace
+{
+public:
+   AlgebraicCoarseSpace(FiniteElementSpace &fine_fes, CeedElemRestriction fine_er,
+                        int order, int dim, int order_reduction_);
+   int GetOrderReduction() const { return order_reduction; }
+   CeedElemRestriction GetCeedElemRestriction() const { return ceed_elem_restriction; }
+   CeedBasis GetCeedCoarseToFine() const { return coarse_to_fine; }
+   virtual const Operator *GetProlongationMatrix() const override { return NULL; }
+   virtual const SparseMatrix *GetRestrictionMatrix() const override { return NULL; }
+   ~AlgebraicCoarseSpace();
+
+protected:
+   int *dof_map;
+   int order_reduction;
+   CeedElemRestriction ceed_elem_restriction;
+   CeedBasis coarse_to_fine;
+};
+
+#ifdef MFEM_USE_MPI
+
+/** @brief Parallel version of AlgebraicCoarseSpace
+
+    This provides prolongation and restriction matrices for RAP-type
+    parallel operators and potential explicit assembly. */
+class ParAlgebraicCoarseSpace : public AlgebraicCoarseSpace
+{
+public:
+   ParAlgebraicCoarseSpace(
+      FiniteElementSpace &fine_fes,
+      CeedElemRestriction fine_er,
+      int order,
+      int dim,
+      int order_reduction_,
+      GroupCommunicator *gc_fine
+   );
+   virtual const Operator *GetProlongationMatrix() const override { return P; }
+   virtual const SparseMatrix *GetRestrictionMatrix() const override { return R_mat; }
+   GroupCommunicator *GetGroupCommunicator() const { return gc; }
+   HypreParMatrix *GetProlongationHypreParMatrix();
+   ~ParAlgebraicCoarseSpace();
+
+private:
+   SparseMatrix *R_mat;
+   GroupCommunicator *gc;
+   ConformingProlongationOperator *P;
+   HypreParMatrix *P_mat;
+   Array<int> ldof_group, ldof_ltdof;
+};
+
+#endif
+
+/** @brief Hierarchy of AlgebraicCoarseSpace objects for use in Multigrid object */
+class AlgebraicSpaceHierarchy : public FiniteElementSpaceHierarchy
+{
+public:
+   /** @brief Construct hierarchy based on finest FiniteElementSpace
+
+       The given space is a real (geometric) space, but the coarse spaces
+       are constructed semi-algebraically with no mesh information. */
+   AlgebraicSpaceHierarchy(FiniteElementSpace &fespace);
+
+   AlgebraicCoarseSpace& GetAlgebraicCoarseSpace(int level)
+   {
+      MFEM_ASSERT(level < GetNumLevels() - 1, "");
+      return static_cast<AlgebraicCoarseSpace&>(*fespaces[level]);
+   }
+
+   ~AlgebraicSpaceHierarchy()
+   {
+      CeedElemRestrictionDestroy(&fine_er);
+      for (int i=0; i<R_tr.Size(); ++i)
+      {
+         delete R_tr[i];
+      }
+      for (int i=0; i<ceed_interpolations.Size(); ++i)
+      {
+         delete ceed_interpolations[i];
+      }
+   }
+
+   /** Prepend an already constructed coarse space to the hierarchy,
+       managing meshes, fespaces, and other arrays appropriately.
+
+       Analogous to FESpaceHierarchy::AddLevel() */
+   void AddCoarseLevel(AlgebraicCoarseSpace* space,
+                       CeedElemRestriction er);
+
+   /// Analogous to FiniteElementSpaceHierarchy::AddOrderRefinedLevel()
+   /// could probably make this happen with just order_reduction if you
+   /// want to save more info in the hierarchy
+   void PrependPCoarsenedLevel(int current_order,
+                               int order_reduction);
+
+private:
+   CeedElemRestriction fine_er;
+   Array<MFEMCeedInterpolation*> ceed_interpolations;
+   Array<TransposeOperator*> R_tr;
+};
+
+/** @brief Extension of Multigrid object to algebraically generated coarse spaces */
+class AlgebraicCeedMultigrid : public GeometricMultigrid
+{
+public:
+   /** @brief Constructs multigrid solver based on existing space hierarchy
+
+       This only works if the Ceed device backend is enabled.
+
+       @param[in] hierarchy   Hierarchy of (algebraic) spaces
+       @param[in] form        partially assembled BilinearForm on finest level
+       @param[in] ess_tdofs   List of essential true dofs on finest level
+       @param[in] print_level 0 is silent
+       @param[in] contrast_threshold Threshold to control p-coarsening
+       @param[in] switch_amg_order Controls when to switch from p-coarsening to AMG
+       @param[in] sparsification controls whether the coarsest grid is "sparsified"
+                                 if it is not already lowest-order
+   */
+   AlgebraicCeedMultigrid(
+      AlgebraicSpaceHierarchy &hierarchy,
+      BilinearForm &form,
+      const Array<int> &ess_tdofs,
+      int print_level=1,
+      double contrast_threshold=1000.0,
+      int switch_amg_order=2,
+      bool collocate_coarse=true,
+      bool sparsification=true,
+      const std::string amgx_config_file=""
+   );
+   virtual void SetOperator(const Operator &op) override { }
+   ~AlgebraicCeedMultigrid();
+
+private:
+   OperatorHandle fine_operator;
+   Array<CeedOperator> ceed_operators;
+};
+
+/** @brief Wrapper for AlgebraicCeedMultigrid object
+
+    This exists so that the algebraic Ceed-based idea has the simplest
+    possible one-line interface. Finer control (choosing smoothers, w-cycle,
+    parameters) can be exercised with the AlgebraicCeedMultigrid object. */
+class AlgebraicCeedSolver : public Solver
+{
+private:
+   AlgebraicSpaceHierarchy fespaces;
+   AlgebraicCeedMultigrid multigrid;
+
+public:
+   /** @brief Constructs algebraic multigrid hierarchy and solver.
+
+       This only works if the Ceed device backend is enabled.
+
+       @param[in] form      partially assembled BilinearForm on finest level
+       @param[in] ess_tdofs List of essential true dofs on finest level
+   */
+   AlgebraicCeedSolver(BilinearForm &form, const Array<int>& ess_tdofs) :
+      fespaces(*form.FESpace()),
+      multigrid(fespaces, form, ess_tdofs)
+   { }
+
+   void Mult(const Vector& x, Vector& y) const { multigrid.Mult(x, y); }
+   void SetOperator(const Operator& op) { multigrid.SetOperator(op); }
+};
+
+} // namespace mfem
+
+#endif // MFEM_USE_CEED
+
+#endif // MFEM_CEED_ALGEBRAIC_HPP

fem/fem.hpp

@@ -62,4 +62,8 @@
 #include "adios2datacollection.hpp"
 #endif
 
+#ifdef MFEM_USE_CEED
+#include "ceed_algebraic.hpp"
+#endif
+
 #endif

fem/fespacehierarchy.hpp

@@ -35,6 +35,11 @@ protected:
 
 public:
 
+   /// Construct an empty finite element space hierarchy. This is useful if the
+   /// hierarchy is constructed by coarsening a fine space, rather than refining
+   /// a coarse space.
+   FiniteElementSpaceHierarchy() { }
+
    /// @brief Constructs a space hierarchy with the given mesh and space on the
    /// coarsest level.
    /** The ownership of the mesh and space may be transferred to the

fem/libceed/ceed-assemble.cpp

@@ -0,0 +1,319 @@
+// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "ceed-assemble.hpp"
+
+#ifdef MFEM_USE_CEED
+#include "ceedsolvers-utility.h"
+
+int CeedHackReallocArray(size_t n, size_t unit, void *p)
+{
+   *(void **)p = realloc(*(void **)p, n*unit);
+   if (n && unit && !*(void **)p)
+      return CeedError(NULL, 1, "realloc failed to allocate %zd members of size "
+                       "%zd\n", n, unit);
+   return 0;
+}
+
+#define CeedHackRealloc(n, p) CeedHackReallocArray((n), sizeof(**(p)), p)
+
+int CeedHackFree(void *p)
+{
+   free(*(void **)p);
+   *(void **)p = NULL;
+   return 0;
+}
+
+namespace mfem
+{
+
+int CeedSingleOperatorFullAssemble(CeedOperator op, SparseMatrix *out)
+{
+   int ierr;
+   Ceed ceed;
+   ierr = CeedOperatorGetCeed(op, &ceed); CeedChk(ierr);
+
+   // Assemble QFunction
+   CeedQFunction qf;
+   ierr = CeedOperatorGetQFunction(op, &qf); CeedChk(ierr);
+   CeedInt numinputfields, numoutputfields;
+   ierr= CeedQFunctionGetNumArgs(qf, &numinputfields, &numoutputfields);
+   CeedChk(ierr);
+   CeedVector assembledqf;
+   CeedElemRestriction rstr_q;
+   ierr = CeedOperatorLinearAssembleQFunction(
+             op, &assembledqf, &rstr_q, CEED_REQUEST_IMMEDIATE); CeedChk(ierr);
+
+   CeedInt qflength;
+   ierr = CeedVectorGetLength(assembledqf, &qflength); CeedChk(ierr);
+
+   CeedOperatorField *input_fields;
+   CeedOperatorField *output_fields;
+   ierr = CeedOperatorGetFields(op, &input_fields, &output_fields); CeedChk(ierr);
+
+   // Determine active input basis
+   CeedQFunctionField *qffields;
+   ierr = CeedQFunctionGetFields(qf, &qffields, NULL); CeedChk(ierr);
+   CeedInt numemodein = 0, ncomp, dim = 1;
+   CeedEvalMode *emodein = NULL;
+   CeedBasis basisin = NULL;
+   CeedElemRestriction rstrin = NULL;
+   for (CeedInt i=0; i<numinputfields; i++)
+   {
+      CeedVector vec;
+      ierr = CeedOperatorFieldGetVector(input_fields[i], &vec); CeedChk(ierr);
+      if (vec == CEED_VECTOR_ACTIVE)
+      {
+         ierr = CeedOperatorFieldGetBasis(input_fields[i], &basisin);
+         CeedChk(ierr);
+         ierr = CeedBasisGetNumComponents(basisin, &ncomp); CeedChk(ierr);
+         ierr = CeedBasisGetDimension(basisin, &dim); CeedChk(ierr);
+         ierr = CeedOperatorFieldGetElemRestriction(input_fields[i], &rstrin);
+         CeedChk(ierr);
+         CeedEvalMode emode;
+         ierr = CeedQFunctionFieldGetEvalMode(qffields[i], &emode);
+         CeedChk(ierr);
+         switch (emode)
+         {
+            case CEED_EVAL_NONE:
+            case CEED_EVAL_INTERP:
+               ierr = CeedHackRealloc(numemodein + 1, &emodein); CeedChk(ierr);
+               emodein[numemodein] = emode;
+               numemodein += 1;
+               break;
+            case CEED_EVAL_GRAD:
+               ierr = CeedHackRealloc(numemodein + dim, &emodein); CeedChk(ierr);
+               for (CeedInt d=0; d<dim; d++)
+               {
+                  emodein[numemodein+d] = emode;
+               }
+               numemodein += dim;
+               break;
+            case CEED_EVAL_WEIGHT:
+            case CEED_EVAL_DIV:
+            case CEED_EVAL_CURL:
+               break; // Caught by QF Assembly
+         }
+      }
+   }
+
+   // Determine active output basis
+   ierr = CeedQFunctionGetFields(qf, NULL, &qffields); CeedChk(ierr);
+   CeedInt numemodeout = 0;
+   CeedEvalMode *emodeout = NULL;
+   CeedBasis basisout = NULL;
+   CeedElemRestriction rstrout = NULL;
+   for (CeedInt i=0; i<numoutputfields; i++)
+   {
+      CeedVector vec;
+      ierr = CeedOperatorFieldGetVector(output_fields[i], &vec); CeedChk(ierr);
+      if (vec == CEED_VECTOR_ACTIVE)
+      {
+         ierr = CeedOperatorFieldGetBasis(output_fields[i], &basisout);
+         CeedChk(ierr);
+         ierr = CeedOperatorFieldGetElemRestriction(output_fields[i], &rstrout);
+         CeedChk(ierr);
+         CeedChk(ierr);
+         CeedEvalMode emode;
+         ierr = CeedQFunctionFieldGetEvalMode(qffields[i], &emode);
+         CeedChk(ierr);
+         switch (emode)
+         {
+            case CEED_EVAL_NONE:
+            case CEED_EVAL_INTERP:
+               ierr = CeedHackRealloc(numemodeout + 1, &emodeout); CeedChk(ierr);
+               emodeout[numemodeout] = emode;
+               numemodeout += 1;
+               break;
+            case CEED_EVAL_GRAD:
+               ierr = CeedHackRealloc(numemodeout + dim, &emodeout); CeedChk(ierr);
+               for (CeedInt d=0; d<dim; d++)
+               {
+                  emodeout[numemodeout+d] = emode;
+               }
+               numemodeout += dim;
+               break;
+            case CEED_EVAL_WEIGHT:
+            case CEED_EVAL_DIV:
+            case CEED_EVAL_CURL:
+               break; // Caught by QF Assembly
+         }
+      }
+   }
+
+   CeedInt nnodes, nelem, elemsize, nqpts;
+   ierr = CeedElemRestrictionGetNumElements(rstrin, &nelem); CeedChk(ierr);
+   ierr = CeedElemRestrictionGetElementSize(rstrin, &elemsize); CeedChk(ierr);
+   ierr = CeedElemRestrictionGetLVectorSize(rstrin, &nnodes); CeedChk(ierr);
+   ierr = CeedBasisGetNumQuadraturePoints(basisin, &nqpts); CeedChk(ierr);
+
+   // Determine elem_dof relation
+   CeedVector index_vec;
+   ierr = CeedVectorCreate(ceed, nnodes, &index_vec); CeedChk(ierr);
+   CeedScalar *array;
+   ierr = CeedVectorGetArray(index_vec, CEED_MEM_HOST, &array); CeedChk(ierr);
+   for (CeedInt i = 0; i < nnodes; ++i)
+   {
+      array[i] = i;
+   }
+   ierr = CeedVectorRestoreArray(index_vec, &array); CeedChk(ierr);
+   CeedVector elem_dof;
+   ierr = CeedVectorCreate(ceed, nelem * elemsize, &elem_dof); CeedChk(ierr);
+   ierr = CeedVectorSetValue(elem_dof, 0.0); CeedChk(ierr);
+   CeedElemRestrictionApply(rstrin, CEED_NOTRANSPOSE, index_vec,
+                            elem_dof, CEED_REQUEST_IMMEDIATE); CeedChk(ierr);
+   const CeedScalar * elem_dof_a;
+   ierr = CeedVectorGetArrayRead(elem_dof, CEED_MEM_HOST, &elem_dof_a);
+   CeedChk(ierr);
+   ierr = CeedVectorDestroy(&index_vec); CeedChk(ierr);
+
+   /// loop over elements and put in SparseMatrix
+   // SparseMatrix * out = new SparseMatrix(nnodes, nnodes);
+   MFEM_ASSERT(out->Height() == nnodes, "Sizes don't match!");
+   MFEM_ASSERT(out->Width() == nnodes, "Sizes don't match!");
+   const CeedScalar *interpin, *gradin;
+   ierr = CeedBasisGetInterp(basisin, &interpin); CeedChk(ierr);
+   ierr = CeedBasisGetGrad(basisin, &gradin); CeedChk(ierr);
+
+   const CeedScalar * assembledqfarray;
+   ierr = CeedVectorGetArrayRead(assembledqf, CEED_MEM_HOST, &assembledqfarray);
+   CeedChk(ierr);
+
+   CeedInt layout[3];
+   ierr = CeedElemRestrictionGetELayout(rstr_q, &layout); CeedChk(ierr);
+   ierr = CeedElemRestrictionDestroy(&rstr_q); CeedChk(ierr);
+
+   // enforce structurally symmetric for later elimination
+   const int skip_zeros = 0;
+   MFEM_ASSERT(numemodein == numemodeout, "My undestanding fails in this case.");
+   for (int e = 0; e < nelem; ++e)
+   {
+      /// get Array<int> for use in SparseMatrix::AddSubMatrix()
+      Array<int> rows(elemsize);
+      for (int i = 0; i < elemsize; ++i)
+      {
+         rows[i] = elem_dof_a[e * elemsize + i];
+      }
+
+      // form element matrix itself
+      DenseMatrix Bmat(nqpts * numemodein, elemsize);
+      Bmat = 0.0;
+      // Store block-diagonal D matrix as collection of small dense blocks
+      DenseTensor Dmat(numemodeout, numemodein, nqpts);
+      Dmat = 0.0;
+      DenseMatrix elem_mat(elemsize, elemsize);
+      elem_mat = 0.0;
+      for (int q = 0; q < nqpts; ++q)
+      {
+         for (int n = 0; n < elemsize; ++n)
+         {
+            CeedInt din = -1;
+            for (int ein = 0; ein < numemodein; ++ein)
+            {
+               if (emodein[ein] == CEED_EVAL_INTERP)
+               {
+                  Bmat(numemodein * q + ein, n) += interpin[q * elemsize + n];
+               }
+               else if (emodein[ein] == CEED_EVAL_GRAD)
+               {
+                  din += 1;
+                  Bmat(numemodein * q + ein, n) += gradin[(din*nqpts+q) * elemsize + n];
+               }
+               else
+               {
+                  MFEM_ASSERT(false, "Not implemented!");
+               }
+            }
+         }
+         for (int ei = 0; ei < numemodein; ++ei)
+         {
+            for (int ej = 0; ej < numemodein; ++ej)
+            {
+               const int comp = ei * numemodein + ej;
+               const int index = q*layout[0] + comp*layout[1] + e*layout[2];
+               Dmat(ei, ej, q) += assembledqfarray[index];
+            }
+         }
+      }
+      DenseMatrix BTD(elemsize, nqpts*numemodein);
+      // Compute B^T*D
+      BTD = 0.0;
+      for (int j=0; j<elemsize; ++j)
+      {
+         for (int q=0; q<nqpts; ++q)
+         {
+            int qq = numemodein*q;
+            for (int ei = 0; ei < numemodein; ++ei)
+            {
+               for (int ej = 0; ej < numemodein; ++ej)
+               {
+                  BTD(j,qq+ei) += Bmat(qq+ej,j)*Dmat(ej,ei,q);
+               }
+            }
+         }
+      }
+
+      Mult(BTD, Bmat, elem_mat);
+
+      /// put element matrix in sparsemat
+      out->AddSubMatrix(rows, rows, elem_mat, skip_zeros);
+   }
+
+   ierr = CeedVectorRestoreArrayRead(elem_dof, &elem_dof_a); CeedChk(ierr);
+   ierr = CeedVectorDestroy(&elem_dof); CeedChk(ierr);
+   ierr = CeedVectorRestoreArrayRead(assembledqf, &assembledqfarray);
+   CeedChk(ierr);
+   ierr = CeedVectorDestroy(&assembledqf); CeedChk(ierr);
+   ierr = CeedHackFree(&emodein); CeedChk(ierr);
+   ierr = CeedHackFree(&emodeout); CeedChk(ierr);
+
+   return 0;
+}
+
+int CeedOperatorFullAssemble(CeedOperator op, SparseMatrix **mat)
+{
+   int ierr;
+
+   CeedElemRestriction er;
+   ierr = CeedOperatorGetActiveElemRestriction(op, &er); CeedChk(ierr);
+   CeedInt nnodes;
+   ierr = CeedElemRestrictionGetLVectorSize(er, &nnodes); CeedChk(ierr);
+
+   SparseMatrix *out = new SparseMatrix(nnodes, nnodes);
+
+   bool isComposite;
+   ierr = CeedOperatorIsComposite(op, &isComposite); CeedChk(ierr);
+   if (isComposite)
+   {
+      CeedInt numsub;
+      CeedOperator *subops;
+      CeedOperatorGetNumSub(op, &numsub);
+      ierr = CeedOperatorGetSubList(op, &subops); CeedChk(ierr);
+      for (int i = 0; i < numsub; ++i)
+      {
+         ierr = CeedSingleOperatorFullAssemble(subops[i], out); CeedChk(ierr);
+      }
+   }
+   else
+   {
+      ierr = CeedSingleOperatorFullAssemble(op, out); CeedChk(ierr);
+   }
+   // enforce structurally symmetric for later elimination
+   const int skip_zeros = 0;
+   out->Finalize(skip_zeros);
+   *mat = out;
+
+   return 0;
+}
+
+} // namespace mfem
+
+#endif

fem/libceed/ceed-assemble.hpp

@@ -0,0 +1,32 @@
+// Copyright (c) 2010-2020, 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_CEED_ASSEMBLE_HPP
+#define MFEM_CEED_ASSEMBLE_HPP
+
+#include "../../config/config.hpp"
+
+#ifdef MFEM_USE_CEED
+#include <ceed.h>
+#include "ceedsolvers-interpolation.h"
+#include "../../linalg/sparsemat.hpp"
+
+namespace mfem
+{
+
+/// Assembles (single processor) a CeedOperator as an mfem::SparseMatrix
+int CeedOperatorFullAssemble(CeedOperator op, SparseMatrix **mat);
+
+} // namespace mfem
+
+#endif
+
+#endif

fem/libceed/ceed-wrappers.cpp

@@ -0,0 +1,221 @@
+// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "ceed-wrappers.hpp"
+
+#ifdef MFEM_USE_CEED
+#include <ceed-backend.h>
+#include "ceed.hpp"
+#include "ceedsolvers-utility.h"
+
+namespace mfem
+{
+
+/** Manages memory for using mfem::Vector s for Ceed operations */
+class MFEMCeedVectorContext
+{
+public:
+   MFEMCeedVectorContext(const mfem::Vector& in, mfem::Vector& out,
+                         CeedVector ceed_in_, CeedVector ceed_out_)
+      :
+      ceed_in(ceed_in_), ceed_out(ceed_out_)
+   {
+      CeedGetPreferredMemType(internal::ceed, &mem);
+      if ( Device::Allows(Backend::DEVICE_MASK) && mem==CEED_MEM_DEVICE )
+      {
+         in_ptr = in.Read();
+         out_ptr = out.ReadWrite();
+      }
+      else
+      {
+         in_ptr = in.HostRead();
+         out_ptr = out.HostReadWrite();
+         mem = CEED_MEM_HOST;
+      }
+
+      CeedVectorSetArray(ceed_in, mem, CEED_USE_POINTER,
+                         const_cast<CeedScalar*>(in_ptr));
+      CeedVectorSetArray(ceed_out, mem, CEED_USE_POINTER, out_ptr);
+   }
+
+   ~MFEMCeedVectorContext()
+   {
+      CeedVectorTakeArray(ceed_in, mem, const_cast<CeedScalar**>(&in_ptr));
+      CeedVectorTakeArray(ceed_out, mem, &out_ptr);
+   }
+
+private:
+   CeedVector ceed_in, ceed_out;
+   const CeedScalar *in_ptr;
+   CeedScalar *out_ptr;
+   CeedMemType mem;
+};
+
+class UnconstrainedMFEMCeedOperator : public Operator
+{
+public:
+   UnconstrainedMFEMCeedOperator(CeedOperator oper);
+   ~UnconstrainedMFEMCeedOperator();
+   virtual void Mult(const Vector& x, Vector& y) const;
+   CeedOperator GetCeedOperator() const { return oper_; }
+   using Operator::SetupRAP;
+private:
+   CeedOperator oper_;
+   CeedVector u_, v_;
+};
+
+UnconstrainedMFEMCeedOperator::UnconstrainedMFEMCeedOperator(CeedOperator oper)
+   : oper_(oper)
+{
+   int ierr = 0;
+   Ceed ceed;
+   ierr += CeedOperatorGetCeed(oper, &ceed);
+   CeedElemRestriction er;
+   ierr += CeedOperatorGetActiveElemRestriction(oper, &er);
+   int s;
+   ierr += CeedElemRestrictionGetLVectorSize(er, &s);
+   height = width = s;
+   ierr += CeedVectorCreate(ceed, height, &v_);
+   ierr += CeedVectorCreate(ceed, width, &u_);
+   MFEM_ASSERT(ierr == 0, "CEED error");
+}
+
+UnconstrainedMFEMCeedOperator::~UnconstrainedMFEMCeedOperator()
+{
+   int ierr = 0;
+   ierr += CeedVectorDestroy(&v_);
+   ierr += CeedVectorDestroy(&u_);
+   MFEM_ASSERT(ierr == 0, "CEED error");
+}
+
+void UnconstrainedMFEMCeedOperator::Mult(const Vector& x, Vector& y) const
+{
+   // would like to use MFEMCeedVectorContext here, does not seem to work
+
+   y = 0.0;
+
+   // I specifically do not want to call the constructor or destructor
+   // of CeedData, this is kind of a hack.
+   CeedData * data = (CeedData*) malloc(sizeof(CeedData));
+   data->u = u_;
+   data->v = v_;
+   data->oper = oper_;
+
+   CeedAddMult(data, x, y);
+
+   free(data);
+}
+
+MFEMCeedOperator::MFEMCeedOperator(
+   CeedOperator oper,
+   const Array<int> &ess_tdofs_,
+   const Operator *P_)
+   : ess_tdofs(ess_tdofs_), P(P_)
+{
+   unconstrained_op = new UnconstrainedMFEMCeedOperator(oper);
+   Operator *rap = unconstrained_op->SetupRAP(P, P);
+   height = width = rap->Height();
+   bool own_rap = (rap != unconstrained_op);
+   constrained_op = new ConstrainedOperator(rap, ess_tdofs, own_rap);
+}
+
+MFEMCeedOperator::MFEMCeedOperator(CeedOperator oper, const Operator *P_)
+   : MFEMCeedOperator(oper, Array<int>(), P_)
+{ }
+
+MFEMCeedOperator::~MFEMCeedOperator()
+{
+   delete constrained_op;
+   delete unconstrained_op;
+}
+
+void MFEMCeedOperator::Mult(const Vector& x, Vector& y) const
+{
+   constrained_op->Mult(x, y);
+}
+
+CeedOperator MFEMCeedOperator::GetCeedOperator() const
+{
+   return unconstrained_op->GetCeedOperator();
+}
+
+const Array<int> &MFEMCeedOperator::GetEssentialTrueDofs() const
+{
+   return ess_tdofs;
+}
+
+const Operator *MFEMCeedOperator::GetProlongation() const
+{
+   return P;
+}
+
+int MFEMCeedInterpolation::Initialize(
+   Ceed ceed, CeedBasis basisctof,
+   CeedElemRestriction erestrictu_coarse, CeedElemRestriction erestrictu_fine)
+{
+   int ierr = 0;
+
+   ierr = CeedInterpolationCreate(ceed, basisctof, erestrictu_coarse,
+                                  erestrictu_fine, &ceed_interp_); CeedChk(ierr);
+
+   ierr = CeedVectorCreate(ceed, height, &v_); CeedChk(ierr);
+   ierr = CeedVectorCreate(ceed, width, &u_); CeedChk(ierr);
+
+   return 0;
+}
+
+MFEMCeedInterpolation::MFEMCeedInterpolation(
+   Ceed ceed, CeedBasis basisctof,
+   CeedElemRestriction erestrictu_coarse,
+   CeedElemRestriction erestrictu_fine)
+{
+   int lo_nldofs, ho_nldofs;
+   CeedElemRestrictionGetLVectorSize(erestrictu_coarse, &lo_nldofs);
+   CeedElemRestrictionGetLVectorSize(erestrictu_fine, &ho_nldofs);
+   height = ho_nldofs;
+   width = lo_nldofs;
+   owns_basis_ = false;
+   Initialize(ceed, basisctof, erestrictu_coarse, erestrictu_fine);
+}
+
+MFEMCeedInterpolation::~MFEMCeedInterpolation()
+{
+   CeedVectorDestroy(&v_);
+   CeedVectorDestroy(&u_);
+   if (owns_basis_)
+   {
+      CeedBasisDestroy(&basisctof_);
+   }
+   CeedInterpolationDestroy(&ceed_interp_);
+}
+
+void MFEMCeedInterpolation::Mult(const mfem::Vector& x, mfem::Vector& y) const
+{
+   int ierr = 0;
+   MFEMCeedVectorContext context(x, y, u_, v_);
+
+   ierr += CeedInterpolationInterpolate(ceed_interp_, u_, v_);
+   MFEM_ASSERT(ierr == 0, "CEED error");
+}
+
+void MFEMCeedInterpolation::MultTranspose(const mfem::Vector& x,
+                                          mfem::Vector& y) const
+{
+   int ierr = 0;
+   MFEMCeedVectorContext context(x, y, v_, u_);
+
+   ierr += CeedInterpolationRestrict(ceed_interp_, v_, u_);
+   MFEM_ASSERT(ierr == 0, "CEED error");
+}
+
+}
+
+#endif

fem/libceed/ceed-wrappers.hpp

@@ -0,0 +1,77 @@
+// Copyright (c) 2010-2020, 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_CEED_WRAPPERS_HPP
+#define MFEM_CEED_WRAPPERS_HPP
+
+#include "../../config/config.hpp"
+
+#ifdef MFEM_USE_CEED
+#include <ceed.h>
+#include "ceedsolvers-interpolation.h"
+#include "../../linalg/operator.hpp"
+
+namespace mfem
+{
+
+/// Wraps a CeedOperator in an mfem::Operator
+class MFEMCeedOperator : public Operator
+{
+public:
+   MFEMCeedOperator(CeedOperator oper, const Array<int> &ess_tdofs_,
+                    const Operator *P_);
+   MFEMCeedOperator(CeedOperator oper, const Operator *P_);
+   ~MFEMCeedOperator();
+   void Mult(const Vector& x, Vector& y) const;
+   CeedOperator GetCeedOperator() const;
+   const Array<int> &GetEssentialTrueDofs() const;
+   const Operator *GetProlongation() const;
+private:
+   Array<int> ess_tdofs;
+   const Operator *P;
+   class UnconstrainedMFEMCeedOperator *unconstrained_op;
+   ConstrainedOperator *constrained_op;
+};
+
+/// Wraps CeedInterpolation object in an mfem::Operator
+class MFEMCeedInterpolation : public mfem::Operator
+{
+public:
+   MFEMCeedInterpolation(
+      Ceed ceed, CeedBasis basisctof,
+      CeedElemRestriction erestrictu_coarse,
+      CeedElemRestriction erestrictu_fine);
+
+   ~MFEMCeedInterpolation();
+
+   virtual void Mult(const mfem::Vector& x, mfem::Vector& y) const;
+
+   virtual void MultTranspose(const mfem::Vector& x, mfem::Vector& y) const;
+
+   using Operator::SetupRAP;
+private:
+   int Initialize(Ceed ceed, CeedBasis basisctof,
+                  CeedElemRestriction erestrictu_coarse,
+                  CeedElemRestriction erestrictu_fine);
+
+   CeedBasis basisctof_;
+   CeedVector u_, v_;
+
+   CeedInterpolation ceed_interp_;
+
+   bool owns_basis_;
+};
+
+}
+
+#endif
+
+#endif

fem/libceed/ceed.cpp

@@ -590,8 +590,8 @@ static void InitCeedTensorBasis(const FiniteElementSpace &fes,
                            qW.GetData(), basis);
 }
 
-static void InitCeedTensorRestriction(const FiniteElementSpace &fes,
-                                      Ceed ceed, CeedElemRestriction *restr)
+void InitCeedTensorRestriction(const FiniteElementSpace &fes,
+                               Ceed ceed, CeedElemRestriction *restr)
 {
    Mesh *mesh = fes.GetMesh();
    const FiniteElement *fe = fes.GetFE(0);

fem/libceed/ceed.hpp

@@ -232,6 +232,15 @@ void InitCeedBasisAndRestriction(const FiniteElementSpace &fes,
 /// Return the path to the libCEED q-function headers.
 const std::string &GetCeedPath();
 
+/// Initialize a CeedBasis and a CeedElemRestriction
+void InitCeedBasisAndRestriction(const FiniteElementSpace &fes,
+                                 const IntegrationRule &ir,
+                                 Ceed ceed, CeedBasis *basis,
+                                 CeedElemRestriction *restr);
+
+void InitCeedTensorRestriction(const FiniteElementSpace &fes,
+                               Ceed ceed, CeedElemRestriction *restr);
+
 // Hash table for CeedBasis
 using CeedBasisKey =
    std::tuple<const FiniteElementSpace*, const IntegrationRule*, int, int, int>;

fem/libceed/ceedsolvers-atpmg.cpp

@@ -0,0 +1,728 @@
+// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "ceedsolvers-atpmg.h"
+#include "ceedsolvers-utility.h"
+
+#ifdef MFEM_USE_CEED
+#include <ceed-backend.h>
+
+#include <math.h>
+// todo: should probably use Ceed memory wrappers instead of calloc/free?
+#include <stdlib.h>
+
+int coarse_1d_edof(int i, int P1d, int coarse_P1d)
+{
+   int coarse_i = (i < coarse_P1d - 1) ? i : -1;
+   if (i == P1d - 1)
+   {
+      coarse_i = coarse_P1d - 1;
+   }
+   return coarse_i;
+}
+
+int reverse_coarse_1d_edof(int i, int P1d, int coarse_P1d)
+{
+   int coarse_i;
+   if (i > P1d - coarse_P1d)
+   {
+      coarse_i = i - (P1d - coarse_P1d);
+   }
+   else
+   {
+      coarse_i = -1;
+   }
+   if (i == 0)
+   {
+      coarse_i = 0;
+   }
+   return coarse_i;
+}
+
+int min4(int a, int b, int c, int d)
+{
+   if (a <= b && a <= c && a <= d)
+   {
+      return a;
+   }
+   else if (b <= a && b <= c && b <= d)
+   {
+      return b;
+   }
+   else if (c <= a && c <= b && c <= d)
+   {
+      return c;
+   }
+   else
+   {
+      return d;
+   }
+}
+
+int CeedATPMGElemRestriction(int order,
+                             int order_reduction,
+                             CeedElemRestriction er_in,
+                             CeedElemRestriction* er_out,
+                             CeedInt *&dof_map)
+{
+   int ierr;
+   Ceed ceed;
+   ierr = CeedElemRestrictionGetCeed(er_in, &ceed); CeedChk(ierr);
+
+   CeedInt numelem, numnodes, numcomp, elemsize;
+   ierr = CeedElemRestrictionGetNumElements(er_in, &numelem); CeedChk(ierr);
+   ierr = CeedElemRestrictionGetLVectorSize(er_in, &numnodes); CeedChk(ierr);
+   ierr = CeedElemRestrictionGetElementSize(er_in, &elemsize); CeedChk(ierr);
+   ierr = CeedElemRestrictionGetNumComponents(er_in, &numcomp); CeedChk(ierr);
+   if (numcomp != 1)
+   {
+      // todo: this will require more thought
+      return CeedError(ceed, 1, "Not implemented!");
+   }
+
+   int P1d = order + 1;
+   int coarse_P1d = P1d - order_reduction;
+   int dim = (log((double) elemsize) / log((double) P1d)) + 1.e-3;
+
+   CeedVector in_lvec, in_evec;
+   ierr = CeedElemRestrictionCreateVector(er_in, &in_lvec, &in_evec);
+   CeedChk(ierr);
+
+   CeedScalar * lvec_data;
+   ierr = CeedVectorGetArray(in_lvec, CEED_MEM_HOST, &lvec_data); CeedChk(ierr);
+   for (int i = 0; i < numnodes; ++i)
+   {
+      lvec_data[i] = (CeedScalar) i;
+   }
+   ierr = CeedVectorRestoreArray(in_lvec, &lvec_data); CeedChk(ierr);
+
+   // todo: I am making assumptions about the ordering of the evec that in
+   // principle are decided by the backend, which I do not control
+   ierr = CeedElemRestrictionApply(er_in, CEED_NOTRANSPOSE, in_lvec, in_evec,
+                                   CEED_REQUEST_IMMEDIATE); CeedChk(ierr);
+   ierr = CeedVectorDestroy(&in_lvec); CeedChk(ierr);
+   const CeedScalar * in_elem_dof;
+   ierr = CeedVectorGetArrayRead(in_evec, CEED_MEM_HOST, &in_elem_dof);
+   CeedChk(ierr);
+
+   // map high-order ldof to low-order ldof
+   // !! caller's responsibility to free
+   dof_map = (CeedInt*) calloc(numnodes, sizeof(CeedInt));
+   for (int i = 0; i < numnodes; ++i)
+   {
+      dof_map[i] = -1;
+   }
+
+   const bool debug = false;
+   CeedInt coarse_elemsize = pow(coarse_P1d, dim);
+   if (debug)
+      printf("coarse_P1d = %d, numelem = %d, coarse_elemsize = %d\n",
+             coarse_P1d, numelem, coarse_elemsize);
+   CeedInt * out_elem_dof = (CeedInt*) calloc(coarse_elemsize * numelem,
+                                              sizeof(CeedInt));
+
+   int running_out_ldof_count = 0;
+   if (dim == 2)
+   {
+      for (int e = 0; e < numelem; ++e)
+      {
+         for (int i = 0; i < P1d; ++i)
+         {
+            for (int j = 0; j < P1d; ++j)
+            {
+               int in_edof = i*P1d + j;
+               int in_ldof = in_elem_dof[e*elemsize + in_edof] + 1.e-10;
+               bool i_edge = (i == 0 || i == P1d - 1);
+               bool j_edge = (j == 0 || j == P1d - 1);
+               int coarse_i, coarse_j;
+               if (i_edge == j_edge)   // vertices and interiors
+               {
+                  // note that interiors could be done with elements in parallel
+                  // (you'd have to rethink numbering but it could be done in advance)
+                  coarse_i = coarse_1d_edof(i, P1d, coarse_P1d);
+                  coarse_j = coarse_1d_edof(j, P1d, coarse_P1d);
+               }
+               else     // edges (without vertices)
+               {
+                  int left_in_edof, left_in_ldof, right_in_edof, right_in_ldof;
+                  if (i_edge)
+                  {
+                     left_in_edof = i*P1d + 0;
+                     right_in_edof = i*P1d + (P1d - 1);
+                     left_in_ldof = in_elem_dof[e*elemsize + left_in_edof] + 1.e-10;
+                     right_in_ldof = in_elem_dof[e*elemsize + right_in_edof] + 1.e-10;
+                     coarse_i = coarse_1d_edof(i, P1d, coarse_P1d);
+                     coarse_j = (left_in_ldof < right_in_ldof) ?
+                                coarse_1d_edof(j, P1d, coarse_P1d) : reverse_coarse_1d_edof(j, P1d, coarse_P1d);
+                  }
+                  else
+                  {
+                     left_in_edof = 0*P1d + j;
+                     right_in_edof = (P1d - 1)*P1d + j;
+                     left_in_ldof = in_elem_dof[e*elemsize + left_in_edof] + 1.e-10;
+                     right_in_ldof = in_elem_dof[e*elemsize + right_in_edof] + 1.e-10;
+                     coarse_i = (left_in_ldof < right_in_ldof) ?
+                                coarse_1d_edof(i, P1d, coarse_P1d) : reverse_coarse_1d_edof(i, P1d, coarse_P1d);
+                     coarse_j = coarse_1d_edof(j, P1d, coarse_P1d);
+                  }
+               }
+               if (coarse_i >= 0 && coarse_j >= 0)
+               {
+                  int out_edof = coarse_i*coarse_P1d + coarse_j;
+                  if (dof_map[in_ldof] >= 0)
+                  {
+                     out_elem_dof[e*coarse_elemsize + out_edof] = dof_map[in_ldof];
+                  }
+                  else
+                  {
+                     out_elem_dof[e*coarse_elemsize + out_edof] = running_out_ldof_count;
+                     dof_map[in_ldof] = running_out_ldof_count;
+                     running_out_ldof_count++;
+                  }
+               }
+            }
+         }
+      }
+   }
+   else if (dim == 3)
+   {
+      // this code is a disaster TODO
+      for (int e = 0; e < numelem; ++e)
+      {
+         if (debug) { printf("  e %d\n", e); }
+         for (int i = 0; i < P1d; ++i)
+         {
+            for (int j = 0; j < P1d; ++j)
+            {
+               for (int k = 0; k < P1d; ++k)
+               {
+                  int in_edof = i*P1d*P1d + j*P1d + k;
+                  int in_ldof = in_elem_dof[e*elemsize + in_edof] + 1.e-10;
+                  int coarse_i, coarse_j, coarse_k;
+                  bool i_edge = (i == 0 || i == P1d - 1);
+                  bool j_edge = (j == 0 || j == P1d - 1);
+                  bool k_edge = (k == 0 || k == P1d - 1);
+                  int topo = 0;
+                  if (i_edge) { topo++; }
+                  if (j_edge) { topo++; }
+                  if (k_edge) { topo++; }
+                  if (topo == 0 || topo == 3)
+                  {
+                     // vertices and interiors
+                     coarse_i = coarse_1d_edof(i, P1d, coarse_P1d);
+                     coarse_j = coarse_1d_edof(j, P1d, coarse_P1d);
+                     coarse_k = coarse_1d_edof(k, P1d, coarse_P1d);
+                  }
+                  else if (topo == 2)
+                  {
+                     // edge
+                     int left_in_edof, left_in_ldof, right_in_edof, right_in_ldof;
+                     if (!i_edge)
+                     {
+                        left_in_edof = 0*P1d*P1d + j*P1d + k;
+                        right_in_edof = (P1d - 1)*P1d*P1d + j*P1d + k;
+                        left_in_ldof = in_elem_dof[e*elemsize + left_in_edof] + 1.e-10;
+                        right_in_ldof = in_elem_dof[e*elemsize + right_in_edof] + 1.e-10;
+                        coarse_i = (left_in_ldof < right_in_ldof) ?
+                                   coarse_1d_edof(i, P1d, coarse_P1d) : reverse_coarse_1d_edof(i, P1d, coarse_P1d);
+                        coarse_j = coarse_1d_edof(j, P1d, coarse_P1d);
+                        coarse_k = coarse_1d_edof(k, P1d, coarse_P1d);
+                     }
+                     else if (!j_edge)
+                     {
+                        left_in_edof = i*P1d*P1d + 0*P1d + k;
+                        right_in_edof = i*P1d*P1d + (P1d - 1)*P1d + k;
+                        left_in_ldof = in_elem_dof[e*elemsize + left_in_edof] + 1.e-10;
+                        right_in_ldof = in_elem_dof[e*elemsize + right_in_edof] + 1.e-10;
+                        coarse_i = coarse_1d_edof(i, P1d, coarse_P1d);
+                        coarse_j = (left_in_ldof < right_in_ldof) ?
+                                   coarse_1d_edof(j, P1d, coarse_P1d) : reverse_coarse_1d_edof(j, P1d, coarse_P1d);
+                        coarse_k = coarse_1d_edof(k, P1d, coarse_P1d);
+                     }
+                     else
+                     {
+                        if (k_edge) { return CeedError(ceed, 1, "Nonsense!"); }
+                        left_in_edof = i*P1d*P1d + j*P1d + 0;
+                        right_in_edof = i*P1d*P1d + j*P1d + (P1d - 1);
+                        left_in_ldof = in_elem_dof[e*elemsize + left_in_edof] + 1.e-10;
+                        right_in_ldof = in_elem_dof[e*elemsize + right_in_edof] + 1.e-10;
+                        coarse_i = coarse_1d_edof(i, P1d, coarse_P1d);
+                        coarse_j = coarse_1d_edof(j, P1d, coarse_P1d);
+                        coarse_k = (left_in_ldof < right_in_ldof) ?
+                                   coarse_1d_edof(k, P1d, coarse_P1d) : reverse_coarse_1d_edof(k, P1d, coarse_P1d);
+                     }
+                  }
+                  else
+                  {
+                     if (topo != 1) { return CeedError(ceed, 1, "Nonsense!"); }
+                     // face
+                     int bottom_left_edof, bottom_right_edof, top_left_edof, top_right_edof;
+                     int bottom_left_ldof, bottom_right_ldof, top_left_ldof, top_right_ldof;
+                     if (i_edge)
+                     {
+                        bottom_left_edof = i*P1d*P1d + 0*P1d + 0;
+                        bottom_right_edof = i*P1d*P1d + 0*P1d + (P1d - 1);
+                        top_right_edof = i*P1d*P1d + (P1d - 1)*P1d + (P1d - 1);
+                        top_left_edof = i*P1d*P1d + (P1d - 1)*P1d + 0;
+                        bottom_left_ldof = in_elem_dof[e*elemsize + bottom_left_edof] + 1.e-10;
+                        bottom_right_ldof = in_elem_dof[e*elemsize + bottom_right_edof] + 1.e-10;
+                        top_right_ldof = in_elem_dof[e*elemsize + top_right_edof] + 1.e-10;
+                        top_left_ldof = in_elem_dof[e*elemsize + top_left_edof] + 1.e-10;
+                        int m = min4(bottom_left_ldof, bottom_right_ldof, top_right_ldof,
+                                     top_left_ldof);
+                        coarse_i = coarse_1d_edof(i, P1d, coarse_P1d);
+                        if (m == bottom_left_ldof)
+                        {
+                           coarse_j = coarse_1d_edof(j, P1d, coarse_P1d);
+                           coarse_k = coarse_1d_edof(k, P1d, coarse_P1d);
+                        }
+                        else if (m == bottom_right_ldof)     // j=0, k=P1d-1
+                        {
+                           coarse_j = coarse_1d_edof(j, P1d, coarse_P1d);
+                           coarse_k = reverse_coarse_1d_edof(k, P1d, coarse_P1d);
+                        }
+                        else if (m == top_right_ldof)
+                        {
+                           coarse_j = reverse_coarse_1d_edof(j, P1d, coarse_P1d);
+                           coarse_k = reverse_coarse_1d_edof(k, P1d, coarse_P1d);
+                        }
+                        else     // j=P1d-1, k=0
+                        {
+                           coarse_j = reverse_coarse_1d_edof(j, P1d, coarse_P1d);
+                           coarse_k = coarse_1d_edof(k, P1d, coarse_P1d);
+                        }
+                     }
+                     else if (j_edge)
+                     {
+                        bottom_left_edof = 0*P1d*P1d + j*P1d + 0;
+                        bottom_right_edof = 0*P1d*P1d + j*P1d + (P1d - 1);
+                        top_right_edof = (P1d - 1)*P1d*P1d + j*P1d + (P1d - 1);
+                        top_left_edof = (P1d - 1)*P1d*P1d + j*P1d + 0;
+                        bottom_left_ldof = in_elem_dof[e*elemsize + bottom_left_edof] + 1.e-10;
+                        bottom_right_ldof = in_elem_dof[e*elemsize + bottom_right_edof] + 1.e-10;
+                        top_right_ldof = in_elem_dof[e*elemsize + top_right_edof] + 1.e-10;
+                        top_left_ldof = in_elem_dof[e*elemsize + top_left_edof] + 1.e-10;
+                        int m = min4(bottom_left_ldof, bottom_right_ldof, top_right_ldof,
+                                     top_left_ldof);
+                        coarse_j = coarse_1d_edof(j, P1d, coarse_P1d);
+                        if (m == bottom_left_ldof)
+                        {
+                           coarse_i = coarse_1d_edof(i, P1d, coarse_P1d);
+                           coarse_k = coarse_1d_edof(k, P1d, coarse_P1d);
+                        }
+                        else if (m == bottom_right_ldof)     // i=0, k=P1d-1
+                        {
+                           coarse_i = coarse_1d_edof(i, P1d, coarse_P1d);
+                           coarse_k = reverse_coarse_1d_edof(k, P1d, coarse_P1d);
+                        }
+                        else if (m == top_right_ldof)
+                        {
+                           coarse_i = reverse_coarse_1d_edof(i, P1d, coarse_P1d);
+                           coarse_k = reverse_coarse_1d_edof(k, P1d, coarse_P1d);
+                        }
+                        else     // i=P1d-1, k=0
+                        {
+                           coarse_i = reverse_coarse_1d_edof(i, P1d, coarse_P1d);
+                           coarse_k = coarse_1d_edof(k, P1d, coarse_P1d);
+                        }
+                     }
+                     else
+                     {
+                        if (!k_edge) { return CeedError(ceed, 1, "Nonsense!"); }
+                        bottom_left_edof = 0*P1d*P1d + 0*P1d + k;
+                        bottom_right_edof = 0*P1d*P1d + (P1d - 1)*P1d + k;
+                        top_right_edof = (P1d - 1)*P1d*P1d + (P1d - 1)*P1d + k;
+                        top_left_edof = (P1d - 1)*P1d*P1d + 0*P1d + k;
+                        bottom_left_ldof = in_elem_dof[e*elemsize + bottom_left_edof] + 1.e-10;
+                        bottom_right_ldof = in_elem_dof[e*elemsize + bottom_right_edof] + 1.e-10;
+                        top_right_ldof = in_elem_dof[e*elemsize + top_right_edof] + 1.e-10;
+                        top_left_ldof = in_elem_dof[e*elemsize + top_left_edof] + 1.e-10;
+                        int m = min4(bottom_left_ldof, bottom_right_ldof, top_right_ldof,
+                                     top_left_ldof);
+                        coarse_k = coarse_1d_edof(k, P1d, coarse_P1d);
+                        if (m == bottom_left_ldof)
+                        {
+                           coarse_i = coarse_1d_edof(i, P1d, coarse_P1d);
+                           coarse_j = coarse_1d_edof(j, P1d, coarse_P1d);
+                        }
+                        else if (m == bottom_right_ldof)     // i=0, j=P1d-1
+                        {
+                           coarse_i = coarse_1d_edof(i, P1d, coarse_P1d);
+                           coarse_j = reverse_coarse_1d_edof(j, P1d, coarse_P1d);
+                        }
+                        else if (m == top_right_ldof)
+                        {
+                           coarse_i = reverse_coarse_1d_edof(i, P1d, coarse_P1d);
+                           coarse_j = reverse_coarse_1d_edof(j, P1d, coarse_P1d);
+                        }
+                        else     // i=P1d-1, j=0
+                        {
+                           coarse_i = reverse_coarse_1d_edof(i, P1d, coarse_P1d);
+                           coarse_j = coarse_1d_edof(j, P1d, coarse_P1d);
+                        }
+                     }
+                  }
+                  if (coarse_i >= 0 && coarse_j >= 0 && coarse_k >= 0)
+                  {
+                     int out_edof = coarse_i*coarse_P1d*coarse_P1d + coarse_j*coarse_P1d + coarse_k;
+                     if (dof_map[in_ldof] >= 0)
+                     {
+                        out_elem_dof[e*coarse_elemsize + out_edof] = dof_map[in_ldof];
+                     }
+                     else
+                     {
+                        out_elem_dof[e*coarse_elemsize + out_edof] = running_out_ldof_count;
+                        dof_map[in_ldof] = running_out_ldof_count;
+                        running_out_ldof_count++;
+                     }
+                     if (debug)
+                        printf("      ci=%d,cj=%d,ck=%d, out_edof=%d, out_ldof=%d\n",
+                               coarse_i, coarse_j, coarse_k, out_edof,
+                               out_elem_dof[e*coarse_elemsize + out_edof]);
+                  }
+               }
+            }
+         }
+      }
+
+   }
+   else
+   {
+      return CeedError(ceed, 1, "Bad dimension!");
+   }
+
+   ierr = CeedVectorRestoreArrayRead(in_evec, &in_elem_dof); CeedChk(ierr);
+   ierr = CeedVectorDestroy(&in_evec); CeedChk(ierr);
+
+   ierr = CeedElemRestrictionCreate(ceed, numelem, coarse_elemsize, numcomp,
+                                    0, running_out_ldof_count,
+                                    CEED_MEM_HOST, CEED_COPY_VALUES, out_elem_dof,
+                                    er_out); CeedChk(ierr);
+
+   free(out_elem_dof);
+
+   return 0;
+}
+
+
+int CeedBasisATPMGCoarseToFine(Ceed ceed, int P1d, int dim, int order_reduction,
+                               CeedBasis *basisc2f)
+{
+   // this assumes Lobatto nodes on fine and coarse again
+   // (not so hard to generalize, but we would have to write it ourselves instead of
+   // calling the following Ceed function)
+   int ierr;
+   ierr = CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, P1d - order_reduction, P1d,
+                                          CEED_GAUSS_LOBATTO, basisc2f); CeedChk(ierr);
+   return 0;
+}
+
+int CeedBasisATPMGCoarseToFine(CeedBasis basisin,
+                               CeedBasis *basisc2f,
+                               int order_reduction)
+{
+   int ierr;
+   Ceed ceed;
+   ierr = CeedBasisGetCeed(basisin, &ceed); CeedChk(ierr);
+
+   CeedInt dim, P1d;
+   ierr = CeedBasisGetDimension(basisin, &dim); CeedChk(ierr);
+   ierr = CeedBasisGetNumNodes1D(basisin, &P1d); CeedChk(ierr);
+   ierr = CeedBasisATPMGCoarseToFine(ceed, P1d, dim, order_reduction,
+                                     basisc2f); CeedChk(ierr);
+   return 0;
+}
+
+int CeedBasisATPMGCoarsen(CeedBasis basisin,
+                          CeedBasis basisc2f,
+                          CeedBasis* basisout,
+                          int order_reduction)
+{
+   int ierr;
+   Ceed ceed;
+   ierr = CeedBasisGetCeed(basisin, &ceed); CeedChk(ierr);
+
+   CeedInt dim, ncomp, P1d, Q1d;
+   ierr = CeedBasisGetDimension(basisin, &dim); CeedChk(ierr);
+   ierr = CeedBasisGetNumComponents(basisin, &ncomp); CeedChk(ierr);
+   ierr = CeedBasisGetNumNodes1D(basisin, &P1d); CeedChk(ierr);
+   ierr = CeedBasisGetNumQuadraturePoints1D(basisin, &Q1d); CeedChk(ierr);
+
+   CeedInt coarse_P1d = P1d - order_reduction;
+
+   const CeedScalar *interp1d;
+   ierr = CeedBasisGetInterp1D(basisin, &interp1d); CeedChk(ierr);
+   const CeedScalar * grad1d;
+   ierr = CeedBasisGetGrad1D(basisin, &grad1d); CeedChk(ierr);
+
+   CeedScalar * coarse_interp1d = (CeedScalar*) calloc(
+                                     coarse_P1d * Q1d, sizeof(CeedScalar));
+   CeedScalar * coarse_grad1d = (CeedScalar*) calloc(
+                                   coarse_P1d * Q1d, sizeof(CeedScalar));
+
+   CeedScalar * fine_nodal_points = (CeedScalar*) calloc(
+                                       P1d, sizeof(CeedScalar));
+
+   // these things are in [-1, 1], not [0, 1], which matters
+   // (todo: how can we determine this or something related, algebraically?)
+   /* one way you might be able to tell is to just run this algorithm
+      with coarse_P1d = 2 (ie, linear) and look for symmetry in the coarse
+      basis matrix? */
+   ierr = CeedLobattoQuadrature(P1d, fine_nodal_points, NULL); CeedChk(ierr);
+   for (int i = 0; i < P1d; ++i)
+   {
+      fine_nodal_points[i] = 0.5 * fine_nodal_points[i] + 0.5; // cheating
+   }
+
+   const CeedScalar *interp_ctof;
+   ierr = CeedBasisGetInterp1D(basisc2f, &interp_ctof); CeedChk(ierr);
+
+   for (int i = 0; i < Q1d; ++i)
+   {
+      for (int j = 0; j < coarse_P1d; ++j)
+      {
+         coarse_interp1d[i * coarse_P1d + j] = 0.0;
+         coarse_grad1d[i * coarse_P1d + j] = 0.0;
+         for (int k = 0; k < P1d; ++k)
+         {
+            coarse_interp1d[i * coarse_P1d + j] += interp_ctof[k * coarse_P1d + j] *
+                                                   interp1d[i * P1d + k];
+
+            // below: you could (a) interpolate the derivative, or (b) take derivative of the interpolant
+            // (b) matches p-multigrid and behaves better, but (a) actually behaves reasonably well
+            // and might be useful in some situations?
+
+            // coarse_grad1d[i * coarse_P1d + j] += grad_ctof[k * coarse_P1d + j] * interp1d[i * P1d + k];
+            coarse_grad1d[i * coarse_P1d + j] += interp_ctof[k * coarse_P1d + j] * grad1d[i
+                                                                                          * P1d + k];
+         }
+      }
+   }
+
+   const CeedScalar * qref1d;
+   ierr = CeedBasisGetQRef(basisin, &qref1d); CeedChk(ierr);
+   const CeedScalar * qweight1d;
+   ierr = CeedBasisGetQWeights(basisin, &qweight1d); CeedChk(ierr);
+   ierr = CeedBasisCreateTensorH1(ceed, dim, ncomp,
+                                  coarse_P1d, Q1d, coarse_interp1d, coarse_grad1d,
+                                  qref1d, qweight1d, basisout); CeedChk(ierr);
+
+   const bool debug = false;
+   if (debug)
+   {
+      for (int q = 0; q < Q1d; q++)
+      {
+         for (int i = 0; i < P1d; i++)
+         {
+            printf("  interp1d[%d,%d] = %f\n",
+                   q, i, interp1d[q * P1d + i]);
+         }
+      }
+      for (int q = 0; q < Q1d; q++)
+      {
+         for (int i = 0; i < coarse_P1d; i++)
+         {
+            printf("  coarse_interp1d[%d,%d] = %f\n",
+                   q, i, coarse_interp1d[q * coarse_P1d + i]);
+         }
+      }
+      for (int q = 0; q < Q1d; q++)
+      {
+         for (int i = 0; i < P1d; i++)
+         {
+            printf("  grad1d[%d,%d] = %f\n",
+                   q, i, grad1d[q * P1d + i]);
+         }
+      }
+      for (int q = 0; q < Q1d; q++)
+      {
+         for (int i = 0; i < coarse_P1d; i++)
+         {
+            printf("  coarse_grad1d[%d,%d] = %f\n",
+                   q, i, coarse_grad1d[q * coarse_P1d + i]);
+         }
+      }
+   }
+
+   free(fine_nodal_points);
+   free(coarse_interp1d);
+   free(coarse_grad1d);
+
+   return 0;
+}
+
+int CeedATPMGOperator(CeedOperator oper, int order_reduction,
+                      CeedElemRestriction coarse_er,
+                      CeedBasis coarse_basis_in,
+                      CeedBasis basis_ctof_in,
+                      CeedOperator* out)
+{
+   int ierr;
+   Ceed ceed;
+   ierr = CeedOperatorGetCeed(oper, &ceed); CeedChk(ierr);
+
+   CeedQFunction qf;
+   ierr = CeedOperatorGetQFunction(oper, &qf); CeedChk(ierr);
+   CeedInt numinputfields, numoutputfields;
+   ierr = CeedQFunctionGetNumArgs(qf, &numinputfields, &numoutputfields);
+   CeedQFunctionField *inputqfields, *outputqfields;
+   ierr = CeedQFunctionGetFields(qf, &inputqfields, &outputqfields); CeedChk(ierr);
+   CeedOperatorField *inputfields, *outputfields;
+   ierr = CeedOperatorGetFields(oper, &inputfields, &outputfields); CeedChk(ierr);
+
+   CeedElemRestriction * er_input =
+      (CeedElemRestriction*) calloc(numinputfields, sizeof(CeedElemRestriction));
+   CeedElemRestriction * er_output =
+      (CeedElemRestriction*) calloc(numoutputfields, sizeof(CeedElemRestriction));
+   CeedVector * if_vector =
+      (CeedVector*) calloc(numinputfields, sizeof(CeedVector));
+   CeedVector * of_vector =
+      (CeedVector*) calloc(numoutputfields, sizeof(CeedVector));
+   CeedBasis * basis_input =
+      (CeedBasis*) calloc(numinputfields, sizeof(CeedBasis));
+   CeedBasis * basis_output =
+      (CeedBasis*) calloc(numoutputfields, sizeof(CeedBasis));
+   CeedBasis cbasis = coarse_basis_in;
+
+   int active_input_basis = -1;
+   for (int i = 0; i < numinputfields; ++i)
+   {
+      ierr = CeedOperatorFieldGetElemRestriction(inputfields[i],
+                                                 &er_input[i]); CeedChk(ierr);
+      ierr = CeedOperatorFieldGetVector(inputfields[i], &if_vector[i]); CeedChk(ierr);
+      ierr = CeedOperatorFieldGetBasis(inputfields[i], &basis_input[i]);
+      CeedChk(ierr);
+      if (if_vector[i] == CEED_VECTOR_ACTIVE)
+      {
+         if (active_input_basis < 0)
+         {
+            active_input_basis = i;
+         }
+         else if (basis_input[i] != basis_input[active_input_basis])
+         {
+            return CeedError(ceed, 1, "Two different active input basis!");
+         }
+      }
+   }
+   for (int i = 0; i < numoutputfields; ++i)
+   {
+      ierr = CeedOperatorFieldGetElemRestriction(outputfields[i],
+                                                 &er_output[i]); CeedChk(ierr);
+      ierr = CeedOperatorFieldGetVector(outputfields[i], &of_vector[i]);
+      CeedChk(ierr);
+      ierr = CeedOperatorFieldGetBasis(outputfields[i], &basis_output[i]);
+      CeedChk(ierr);
+      if (of_vector[i] == CEED_VECTOR_ACTIVE)
+      {
+         // should already be coarsened
+         if (basis_output[i] != basis_input[active_input_basis])
+         {
+            return CeedError(ceed, 1, "Input and output basis do not match!");
+         }
+         if (er_output[i] != er_input[active_input_basis])
+         {
+            return CeedError(ceed, 1, "Input and output elem-restriction do not match!");
+         }
+      }
+   }
+
+   CeedOperator coper;
+   ierr = CeedOperatorCreate(ceed, qf, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE,
+                             &coper); CeedChk(ierr);
+
+   for (int i = 0; i < numinputfields; ++i)
+   {
+      char * fieldname;
+      ierr = CeedQFunctionFieldGetName(inputqfields[i], &fieldname); CeedChk(ierr);
+      if (if_vector[i] == CEED_VECTOR_ACTIVE)
+      {
+         ierr = CeedOperatorSetField(coper, fieldname, coarse_er, cbasis,
+                                     if_vector[i]); CeedChk(ierr);
+      }
+      else
+      {
+         ierr = CeedOperatorSetField(coper, fieldname, er_input[i], basis_input[i],
+                                     if_vector[i]); CeedChk(ierr);
+      }
+   }
+   for (int i = 0; i < numoutputfields; ++i)
+   {
+      char * fieldname;
+      ierr = CeedQFunctionFieldGetName(outputqfields[i], &fieldname); CeedChk(ierr);
+      if (of_vector[i] == CEED_VECTOR_ACTIVE)
+      {
+         ierr = CeedOperatorSetField(coper, fieldname, coarse_er, cbasis,
+                                     of_vector[i]); CeedChk(ierr);
+      }
+      else
+      {
+         ierr = CeedOperatorSetField(coper, fieldname, er_output[i], basis_output[i],
+                                     of_vector[i]); CeedChk(ierr);
+      }
+   }
+   free(er_input);
+   free(er_output);
+   free(if_vector);
+   free(of_vector);
+   free(basis_input);
+   free(basis_output);
+
+   *out = coper;
+   return 0;
+}
+
+
+int CeedATPMGOperator(CeedOperator oper, int order_reduction,
+                      CeedElemRestriction coarse_er,
+                      CeedBasis *coarse_basis_out,
+                      CeedBasis *basis_ctof_out,
+                      CeedOperator *out)
+{
+   int ierr;
+
+   CeedQFunction qf;
+   ierr = CeedOperatorGetQFunction(oper, &qf); CeedChk(ierr);
+   CeedInt numinputfields, numoutputfields;
+   ierr = CeedQFunctionGetNumArgs(qf, &numinputfields, &numoutputfields);
+   CeedOperatorField *inputfields;
+   ierr = CeedOperatorGetFields(oper, &inputfields, NULL); CeedChk(ierr);
+
+   CeedBasis basis;
+   ierr = CeedOperatorGetActiveBasis(oper, &basis); CeedChk(ierr);
+   ierr = CeedBasisATPMGCoarseToFine(basis, basis_ctof_out, order_reduction);
+   CeedChk(ierr);
+   ierr = CeedBasisATPMGCoarsen(basis, *basis_ctof_out, coarse_basis_out,
+                                order_reduction); CeedChk(ierr);
+   ierr = CeedATPMGOperator(oper, order_reduction, coarse_er, *coarse_basis_out,
+                            *basis_ctof_out, out); CeedChk(ierr);
+   return 0;
+}
+
+int CeedATPMGBundle(CeedOperator oper, int order_reduction,
+                    CeedBasis* coarse_basis_out,
+                    CeedBasis* basis_ctof_out,
+                    CeedElemRestriction* er_out,
+                    CeedOperator* coarse_oper,
+                    CeedInt *&dof_map)
+{
+   int ierr;
+   CeedInt order;
+   ierr = CeedOperatorGetOrder(oper, &order); CeedChk(ierr);
+   CeedElemRestriction ho_er;
+   ierr = CeedOperatorGetActiveElemRestriction(oper, &ho_er); CeedChk(ierr);
+   ierr = CeedATPMGElemRestriction(order, order_reduction, ho_er, er_out, dof_map);
+   CeedChk(ierr);
+   ierr = CeedATPMGOperator(oper, order_reduction, *er_out, coarse_basis_out,
+                            basis_ctof_out, coarse_oper); CeedChk(ierr);
+   return 0;
+}
+
+#endif // MFEM_USE_CEED

fem/libceed/ceedsolvers-atpmg.h

@@ -0,0 +1,91 @@
+// Copyright (c) 2010-2020, 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_CEEDSOLVERS_ATPMG_H
+#define MFEM_CEEDSOLVERS_ATPMG_H
+
+#include "../../config/config.hpp"
+
+#ifdef MFEM_USE_CEED
+
+#include <ceed.h>
+
+/** @brief Take given (high-order) CeedElemRestriction and make a new
+    CeedElemRestriction, which corresponds to a lower-order problem.
+
+    Assumes a Gauss-Lobatto basis and tensor product elements.
+
+    This is a setup routine that operates on the host. */
+int CeedATPMGElemRestriction(int order,
+                             int order_reduction,
+                             CeedElemRestriction er_in,
+                             CeedElemRestriction* er_out,
+                             CeedInt *&dof_map);
+
+/** @brief Create coarse-to-fine basis, given number of input nodes and order reduction.
+
+    Assumes Gauss-Lobatto basis. This is useful because it does not require an
+    input CeedBasis object, which depends on choice of quadrature rule, whereas
+    the coarse-to-fine operator is independent of quadrature. */
+int CeedBasisATPMGCoarseToFine(Ceed ceed, int P1d, int dim, int order_reduction,
+                               CeedBasis *basisc2f);
+
+/** @brief Given basis basisin, reduces its order by order_reduction and
+    return basisout (which has the same height (Q1d) but is narrower
+    (smaller P1d))
+
+    The algorithm we describe in the writeup takes the locations of
+    the fine nodes as input, but this particular implementation simply
+    assumes Gauss-Lobatto, and furthermore assumes the MFEM [0, 1]
+    reference element (rather than the Ceed/Petsc [-1, 1] element) */
+int CeedBasisATPMGCoarsen(CeedBasis basisin, CeedBasis* basisout,
+                          CeedBasis* basis_ctof,
+                          int order_reduction);
+
+/** @brief Coarsen a CeedOperator using semi-algebraic p-multigrid
+
+    This implementation does not coarsen the integration points at all.
+
+    @param[in] oper              the operator to coarsen
+    @param[in] order_reduction   how much to coarsen (order p)
+    @param[in] coarse_er         CeedElemRestriction for coarse operator
+                                 (see CeedATPMGElemRestriction)
+    @param[out] coarse_basis_out CeedBasis for coarser operator
+    @param[out] out              coarsened CeedOperator
+*/
+int CeedATPMGOperator(CeedOperator oper, int order_reduction,
+                      CeedElemRestriction coarse_er,
+                      CeedBasis* coarse_basis_out,
+                      CeedBasis* basis_ctof_out,
+                      CeedOperator* out);
+
+/** @brief Given (fine) CeedOperator, produces everything you need for a coarse
+    level (operator and interpolation.
+
+    @param[in]  oper             Fine CeedOperator to coarsen
+    @param[in]  order_reduction  Amount to reduce the order (p) of the operator
+    @param[out] coarse_basis_out CeedBasis for coarse operator
+    @param[out] basis_ctof_out   CeedBasis describing interpolation from coarse to fine
+    @param[out] er_out           CeedElemRestriction for coarse operator
+    @param[out] coarse_oper      coarse operator itself
+    @param[out] dof_map          maps high-order ldof to low-order ldof, needed for
+                                 further coarsening
+*/
+int CeedATPMGBundle(CeedOperator oper, int order_reduction,
+                    CeedBasis* coarse_basis_out,
+                    CeedBasis* basis_ctof_out,
+                    CeedElemRestriction* er_out,
+                    CeedOperator* coarse_oper,
+                    CeedInt *&dof_map);
+
+#endif // MFEM_USE_CEED
+
+#endif // include guard

fem/libceed/ceedsolvers-interpolation.cpp

@@ -0,0 +1,167 @@
+// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "ceedsolvers-interpolation.h"
+#include "ceedsolvers-utility.h"
+
+#include "../../general/forall.hpp"
+using namespace mfem;
+
+
+#ifdef MFEM_USE_CEED
+#include <stdlib.h>
+
+int CeedInterpolationCreate(Ceed ceed, CeedBasis basisctof,
+                            CeedElemRestriction erestrictu_coarse,
+                            CeedElemRestriction erestrictu_fine,
+                            CeedInterpolation *interp)
+{
+   int ierr;
+
+   int height, width;
+   ierr = CeedElemRestrictionGetLVectorSize(erestrictu_coarse, &width);
+   CeedChk(ierr);
+   CeedElemRestrictionGetLVectorSize(erestrictu_fine, &height); CeedChk(ierr);
+
+   // interpolation qfunction
+   const int bp3_ncompu = 1;
+   CeedQFunction qf_restrict, qf_prolong;
+   ierr = CeedQFunctionCreateIdentity(ceed, bp3_ncompu, CEED_EVAL_NONE,
+                                      CEED_EVAL_INTERP, &qf_restrict); CeedChk(ierr);
+   ierr = CeedQFunctionCreateIdentity(ceed, bp3_ncompu, CEED_EVAL_INTERP,
+                                      CEED_EVAL_NONE, &qf_prolong); CeedChk(ierr);
+
+   CeedVector c_fine_multiplicity;
+   ierr = CeedVectorCreate(ceed, height, &c_fine_multiplicity); CeedChk(ierr);
+   ierr = CeedVectorSetValue(c_fine_multiplicity, 0.0); CeedChk(ierr);
+
+   // Create the restriction operator
+   // Restriction - Fine to coarse
+   CeedOperator op_interp, op_restrict;
+   ierr = CeedOperatorCreate(ceed, qf_restrict, CEED_QFUNCTION_NONE,
+                             CEED_QFUNCTION_NONE, &op_restrict); CeedChk(ierr);
+   ierr = CeedOperatorSetField(op_restrict, "input", erestrictu_fine,
+                               CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); CeedChk(ierr);
+   ierr = CeedOperatorSetField(op_restrict, "output", erestrictu_coarse,
+                               basisctof, CEED_VECTOR_ACTIVE); CeedChk(ierr);
+
+   // Interpolation - Coarse to fine
+   // Create the prolongation operator
+   ierr = CeedOperatorCreate(ceed, qf_prolong, CEED_QFUNCTION_NONE,
+                             CEED_QFUNCTION_NONE, &op_interp); CeedChk(ierr);
+   CeedOperatorSetField(op_interp, "input", erestrictu_coarse,
+                        basisctof, CEED_VECTOR_ACTIVE);
+   CeedOperatorSetField(op_interp, "output", erestrictu_fine,
+                        CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE);
+
+   CeedVector fine_multiplicity_r;
+   ierr = CeedElemRestrictionGetMultiplicity(
+             erestrictu_fine, c_fine_multiplicity); CeedChk(ierr);
+   ierr = CeedVectorCreate(ceed, height, &fine_multiplicity_r); CeedChk(ierr);
+
+   CeedScalar* fine_r_data;
+   const CeedScalar* fine_data;
+   ierr = CeedVectorGetArray(fine_multiplicity_r, CEED_MEM_HOST,
+                             &fine_r_data); CeedChk(ierr);
+   ierr = CeedVectorGetArrayRead(c_fine_multiplicity, CEED_MEM_HOST,
+                                 &fine_data); CeedChk(ierr);
+   MFEM_FORALL(i, height,
+   {fine_r_data[i] = 1.0 / fine_data[i];});
+
+   ierr = CeedVectorRestoreArray(fine_multiplicity_r, &fine_r_data); CeedChk(ierr);
+   ierr = CeedVectorRestoreArrayRead(c_fine_multiplicity, &fine_data);
+   CeedChk(ierr);
+   ierr = CeedVectorDestroy(&c_fine_multiplicity); CeedChk(ierr);
+
+   CeedVector fine_work;
+   ierr = CeedVectorCreate(ceed, height, &fine_work); CeedChk(ierr);
+
+   *interp = (CeedInterpolation) calloc(1,
+                                        sizeof(struct CeedInterpolation_private));
+   (*interp)->ceed = ceed;
+   (*interp)->qf_restrict = qf_restrict;
+   (*interp)->qf_prolong = qf_prolong;
+   (*interp)->op_interp = op_interp;
+   (*interp)->op_restrict = op_restrict;
+   (*interp)->fine_multiplicity_r = fine_multiplicity_r;
+   (*interp)->fine_work = fine_work;
+
+   return 0;
+}
+
+/// this is not implemented with reference counting etc.
+int CeedInterpolationDestroy(CeedInterpolation *interp)
+{
+   int ierr;
+
+   ierr = CeedQFunctionDestroy(&(*interp)->qf_restrict); CeedChk(ierr);
+   ierr = CeedQFunctionDestroy(&(*interp)->qf_prolong); CeedChk(ierr);
+   ierr = CeedOperatorDestroy(&(*interp)->op_interp); CeedChk(ierr);
+   ierr = CeedOperatorDestroy(&(*interp)->op_restrict); CeedChk(ierr);
+   ierr = CeedVectorDestroy(&(*interp)->fine_multiplicity_r); CeedChk(ierr);
+   ierr = CeedVectorDestroy(&(*interp)->fine_work); CeedChk(ierr);
+
+   free(*interp);
+   return 0;
+}
+
+/// @todo could use a CEED_REQUEST here
+int CeedInterpolationInterpolate(CeedInterpolation interp,
+                                 CeedVector in, CeedVector out)
+{
+   int ierr;
+
+   ierr = CeedOperatorApply(interp->op_interp, in, out,
+                            CEED_REQUEST_IMMEDIATE); CeedChk(ierr);
+   ierr = CeedVectorPointwiseMult(out, interp->fine_multiplicity_r); CeedChk(ierr);
+
+   return 0;
+}
+
+/// @todo could use a CEED_REQUEST here
+/// @todo using MFEM_FORALL in this Ceed-like function is ugly
+int CeedInterpolationRestrict(CeedInterpolation interp,
+                              CeedVector in, CeedVector out)
+{
+   int ierr;
+
+   int length;
+   ierr = CeedVectorGetLength(in, &length); CeedChk(ierr);
+
+   const CeedScalar *multiplicitydata, *indata;
+   CeedScalar *workdata;
+   CeedMemType mem;
+   if (Device::Allows(Backend::DEVICE_MASK))
+   {
+      mem = CEED_MEM_DEVICE;
+   }
+   else
+   {
+      mem = CEED_MEM_HOST;
+   }
+   ierr = CeedVectorGetArrayRead(in, mem, &indata); CeedChk(ierr);
+   ierr = CeedVectorGetArrayRead(interp->fine_multiplicity_r, mem,
+                                 &multiplicitydata); CeedChk(ierr);
+   ierr = CeedVectorGetArray(interp->fine_work, mem, &workdata); CeedChk(ierr);
+   MFEM_FORALL(i, length,
+   {workdata[i] = indata[i] * multiplicitydata[i];});
+   ierr = CeedVectorRestoreArrayRead(in, &indata); CeedChk(ierr);
+   ierr = CeedVectorRestoreArrayRead(interp->fine_multiplicity_r,
+                                     &multiplicitydata); CeedChk(ierr);
+   ierr = CeedVectorRestoreArray(interp->fine_work, &workdata); CeedChk(ierr);
+
+   ierr = CeedOperatorApply(interp->op_restrict, interp->fine_work, out,
+                            CEED_REQUEST_IMMEDIATE); CeedChk(ierr);
+
+   return 0;
+}
+
+#endif // MFEM_USE_CEED

fem/libceed/ceedsolvers-interpolation.h

@@ -0,0 +1,48 @@
+// Copyright (c) 2010-2020, 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_CEEDSOLVERS_INTERPOLATION_H
+#define MFEM_CEEDSOLVERS_INTERPOLATION_H
+
+#include "../../config/config.hpp"
+
+#ifdef MFEM_USE_CEED
+#include <ceed.h>
+
+typedef struct CeedInterpolation_private *CeedInterpolation;
+
+/** @brief Multigrid interpolation operator in Ceed framework
+
+    Interpolation/restriction has two components, an element-wise
+    interpolation and then a scaling to correct multiplicity
+    on shared ldofs. This encapsulates those two in one object. */
+struct CeedInterpolation_private {
+  Ceed ceed;
+  CeedQFunction qf_restrict, qf_prolong;
+  CeedOperator op_interp, op_restrict;
+  CeedVector fine_multiplicity_r;
+  CeedVector fine_work;
+};
+
+int CeedInterpolationCreate(Ceed ceed, CeedBasis basisctof,
+                            CeedElemRestriction erestrictu_coarse,
+                            CeedElemRestriction erestrictu_fine,
+                            CeedInterpolation *interp);
+int CeedInterpolationDestroy(CeedInterpolation *interp);
+
+int CeedInterpolationInterpolate(CeedInterpolation interp,
+                                 CeedVector in, CeedVector out);
+int CeedInterpolationRestrict(CeedInterpolation interp,
+                              CeedVector in, CeedVector out);
+
+#endif // MFEM_USE_CEED
+
+#endif // include guard

fem/libceed/ceedsolvers-qcoarsen.cpp

@@ -0,0 +1,553 @@
+// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "ceedsolvers-qcoarsen.h"
+#include "ceedsolvers-utility.h"
+
+#ifdef MFEM_USE_CEED
+#include <ceed-backend.h>
+
+// todo: should probably use Ceed memory wrappers instead of calloc/free?
+#include <stdlib.h>
+#include "linear.h"
+
+#include <math.h>  // for fabs(), which I actualy don't like
+
+#include <fstream>
+
+/** @brief Coarsen the rows (integration points) of a CeedBasis
+
+    Originally thought this would be like CeedBasisATPMGCoarsen(),
+    but the "interpolation" in quadrature points is not quite
+    analogous, so I think I'm going to do something simpler but less
+    algebraic.
+
+    We will need something like qbasisctof in order to coarsen the
+    linearassembled vector, in fact qbasisctof is essentially P_Q
+    and is essential to Pazner's theorem.
+
+    For the first coarsening, fine_quadmode = CEED_GAUSS (almost certainly)
+    You can decide the coarse_quadmode, CEED_GAUSS_LOBATTO leads to collocated B operators
+    For later coarsenings, fine_quadmode has to match what you did before for coarse
+
+    @param[in] basisin   the CeedBasis, already p-coarsened
+    @param[out] basisout new CeedBasis, same width, but shorter
+    @param[out] qbasisctof describes interpolation of integration points (P_Q)
+    @param[in] order_reduction  amount to coarsen
+    @param[in] collocated_coarse  whether to use collocated quadrature on coarser level
+    @param[in] fine_quadmode  points for fine quadrature rule
+*/
+int CeedBasisQCoarsen(CeedBasis basisin, CeedBasis* basisout,
+                      CeedBasis* qbasisctof,
+                      int order_reduction,
+                      CeedQuadMode fine_quadmode,
+                      CeedQuadMode coarse_quadmode)
+{
+   int ierr;
+   Ceed ceed;
+   ierr = CeedBasisGetCeed(basisin, &ceed); CeedChk(ierr);
+
+   CeedInt dim, ncomp, P1d, Q1d;
+   ierr = CeedBasisGetDimension(basisin, &dim); CeedChk(ierr);
+   ierr = CeedBasisGetNumComponents(basisin, &ncomp); CeedChk(ierr);
+   ierr = CeedBasisGetNumNodes1D(basisin, &P1d); CeedChk(ierr);
+   ierr = CeedBasisGetNumQuadraturePoints1D(basisin, &Q1d); CeedChk(ierr);
+
+   CeedInt coarse_Q1d = Q1d - order_reduction;
+
+   // the "grad" part of qbasisctof will be meaningless, we only use the
+   // "interp" part.
+   if (coarse_quadmode == CEED_GAUSS)
+   {
+      ierr = CeedBasisCreateTensorH1Gauss(ceed, dim, ncomp, coarse_Q1d, Q1d,
+                                          fine_quadmode, qbasisctof); CeedChk(ierr);
+   }
+   else if (coarse_quadmode == CEED_GAUSS_LOBATTO)
+   {
+      ierr = CeedBasisCreateTensorH1Lagrange(ceed, dim, ncomp, coarse_Q1d, Q1d,
+                                             fine_quadmode, qbasisctof); CeedChk(ierr);
+   }
+   else
+   {
+      return CeedError(ceed, 1, "Bad quadrature mode!");
+   }
+
+   // the Ceed reference element is [-1, 1], while the MFEM element is [0, 1]
+   // which means with order_reduction=0 we get different gradients in the
+   // bases; should actually interpolate or something...
+   ierr = CeedBasisCreateMFEMTensorH1Lagrange(ceed, dim, ncomp, P1d, coarse_Q1d,
+                                              coarse_quadmode, basisout); CeedChk(ierr);
+
+   return 0;
+}
+
+/** Given the (CeedVector) output of CeedOperatorLinearAssembleQFunction,
+    coarsen it according to the local integration-interpolation qbasisctof
+    and the numbering encoded in rstr_q, returning a (smaller)
+    CeedVector with the coarsened D operator. */
+int CeedQFunctionCoarsenAssembledVector(CeedVector assembledqf,
+                                        CeedElemRestriction rstr_q,
+                                        CeedBasis qbasisctof,
+                                        CeedVector* coarse_assembledqf_out)
+{
+   int ierr;
+   Ceed ceed;
+   ierr = CeedVectorGetCeed(assembledqf, &ceed); CeedChk(ierr);
+
+   CeedInt layout[3];
+   ierr = CeedElemRestrictionGetELayout(rstr_q, &layout); CeedChk(ierr);
+   CeedInt qflength;
+   ierr = CeedVectorGetLength(assembledqf, &qflength); CeedChk(ierr);
+
+   // may want to ignore layout *entirely* here
+   int elemsize, lsize, ncomp, nelem;
+   ierr = CeedElemRestrictionGetLVectorSize(rstr_q, &lsize); CeedChk(ierr);
+   ierr = CeedElemRestrictionGetNumElements(rstr_q, &nelem); CeedChk(ierr);
+   ierr = CeedElemRestrictionGetElementSize(rstr_q, &elemsize); CeedChk(ierr);
+   ierr = CeedElemRestrictionGetNumComponents(rstr_q, &ncomp); CeedChk(ierr);  
+
+   CeedInt coarse_qflength;
+   CeedInt P, Q, basis_ncomp;
+   ierr = CeedBasisGetNumNodes(qbasisctof, &P); CeedChk(ierr);
+   ierr = CeedBasisGetNumQuadraturePoints(qbasisctof, &Q); CeedChk(ierr);
+   ierr = CeedBasisGetNumComponents(qbasisctof, &basis_ncomp); CeedChk(ierr);
+
+   if (qflength != ncomp * nelem * Q)
+   {
+      return CeedError(ceed, 1, "original qfunction vector does not match rstr_q!");
+   }
+
+   // note well ncomp (not basis_ncomp) on line below (they are different)
+   coarse_qflength = ncomp * nelem * P;
+   if (Q != elemsize)
+   {
+      return CeedError(ceed, 1, "qbasisctof does not match rstr_q!");
+   }
+
+   CeedVector coarse_assembledqf;
+   ierr = CeedVectorCreate(ceed, coarse_qflength, &coarse_assembledqf); CeedChk(ierr);
+   ierr = CeedVectorSetValue(coarse_assembledqf, 0.0); CeedChk(ierr);
+
+   const CeedScalar* finedata;
+   CeedScalar* coarsedata;
+   ierr = CeedVectorGetArrayRead(assembledqf, CEED_MEM_HOST, &finedata); CeedChk(ierr);
+   ierr = CeedVectorGetArray(coarse_assembledqf, CEED_MEM_HOST, &coarsedata); CeedChk(ierr);
+
+   // rows associated with fine, cols associated with coarse
+   // they're both quadpoints, but interface thinks rows(fine) are quad, cols(coarse) are basis
+   // also note I am applying the *transpose*
+   const CeedScalar* ctof_interp;
+   ierr = CeedBasisGetInterp(qbasisctof, &ctof_interp); CeedChk(ierr);
+   for (int k = 0; k < coarse_qflength; ++k)
+   {
+      coarsedata[k] = 0.0;
+   }
+   const int d_per_fineelem = ncomp * Q;
+   const int d_per_coarseelem = ncomp * P;
+
+   for (int e = 0; e < nelem; ++e)
+   {
+      for (int j = 0; j < P; ++j)  // associated with coarse
+      {
+         for (int c = 0; c < ncomp; ++c)
+         {
+            const int output_index = e*d_per_coarseelem + c*P + j;
+            for (int i = 0; i < Q; ++i) // associated with fine
+            {
+               const int input_index = e*d_per_fineelem + c*Q + i;
+               coarsedata[output_index] +=
+                  ctof_interp[i*P + j] * finedata[input_index];
+            }
+         }
+      }
+   }
+
+   ierr = CeedVectorRestoreArrayRead(assembledqf, &finedata); CeedChk(ierr);
+   ierr = CeedVectorRestoreArray(coarse_assembledqf, &coarsedata); CeedChk(ierr);
+
+   *coarse_assembledqf_out = coarse_assembledqf;
+
+   return 0;
+}
+
+/** Given an ElemRestriction rstr_q and a local integration-interpolation,
+    return a CeedElemRestriction with the integration points coarsened.
+
+    @param[in] rstr_q
+    @param[in] qbasisctof (only used for dimensions/sizes!)
+    @param[out] coarse_rstr_q
+    @param[out] ncomp
+*/
+int CeedElementRestrictionQCoarsen(CeedElemRestriction rstr_q,
+                                   CeedBasis qbasisctof,
+                                   CeedElemRestriction* coarse_rstr_q,
+                                   CeedInt* ncomp)
+{
+   int ierr;
+   Ceed ceed;
+   ierr = CeedElemRestrictionGetCeed(rstr_q, &ceed); CeedChk(ierr);
+
+   // layout and strides are different; we may only care about strides
+   CeedInt layout[3];
+   ierr = CeedElemRestrictionGetELayout(rstr_q, &layout); CeedChk(ierr);
+   CeedInt strides[3];
+   ierr = CeedElemRestrictionGetStrides(rstr_q, &strides); CeedChk(ierr);
+
+   CeedInt coarse_lsize;
+   CeedInt coarse_strides[3];
+
+   // some of these are only used for sanity checking
+   CeedInt q_nelem, q_elemsize, q_lsize, q_ncomp;
+   ierr = CeedElemRestrictionGetNumElements(rstr_q, &q_nelem); CeedChk(ierr);
+   ierr = CeedElemRestrictionGetElementSize(rstr_q, &q_elemsize); CeedChk(ierr);
+   ierr = CeedElemRestrictionGetLVectorSize(rstr_q, &q_lsize); CeedChk(ierr);
+   ierr = CeedElemRestrictionGetNumComponents(rstr_q, &q_ncomp); CeedChk(ierr);
+   *ncomp = q_ncomp;
+
+   CeedInt coarse_elemsize, Q; // Q is *fine* quadpoints, coarse_elemsize is *coarse* quadpoints
+   ierr = CeedBasisGetNumNodes(qbasisctof, &coarse_elemsize); CeedChk(ierr);
+   ierr = CeedBasisGetNumQuadraturePoints(qbasisctof, &Q); CeedChk(ierr);
+
+   coarse_lsize = q_nelem * coarse_elemsize * q_ncomp;
+
+   coarse_strides[0] = strides[0];  // always 1, as far as I can tell
+   if (strides[1] == q_elemsize && strides[2] == q_elemsize * q_ncomp)
+   {
+      // characteristic of host layout/strides?
+      coarse_strides[1] = coarse_elemsize;
+      coarse_strides[2] = q_ncomp * coarse_elemsize;
+   }
+   else if (strides[1] == q_elemsize * q_nelem && strides[2] == q_elemsize)
+   {
+      // characteristic of device layout/strides?
+      coarse_strides[1] = q_nelem * coarse_elemsize;
+      coarse_strides[2] = coarse_elemsize;
+   }
+   else
+   {
+      return CeedError(ceed, 1, "I do not understand");
+   }
+   ierr = CeedElemRestrictionCreateStrided(ceed, q_nelem, coarse_elemsize, q_ncomp,
+                                           coarse_lsize, coarse_strides,
+                                           coarse_rstr_q); CeedChk(ierr); 
+
+   return 0;
+}
+
+/// Does not work for composite operator
+int CeedSingleOperatorGetHeuristics(CeedOperator oper, CeedScalar* minq,
+                                    CeedScalar* maxq, CeedScalar* absmin)
+{
+   int ierr;
+   CeedQFunction qfin;
+   ierr = CeedOperatorGetQFunction(oper, &qfin); CeedChk(ierr);
+
+   CeedVector assembledqf;
+   CeedElemRestriction rstr_q;
+   ierr = CeedOperatorLinearAssembleQFunction(
+      oper, &assembledqf, &rstr_q, CEED_REQUEST_IMMEDIATE); CeedChk(ierr);
+
+   CeedInt assembledqf_len;
+   ierr = CeedVectorGetLength(assembledqf, &assembledqf_len); CeedChk(ierr);
+   const CeedScalar * tempdata;
+   ierr = CeedVectorGetArrayRead(assembledqf, CEED_MEM_HOST, &tempdata); CeedChk(ierr);
+
+   *minq = 1.e+12;
+   *maxq = -1.e+12;
+   *absmin = 1.e+12;
+   for (int i = 0; i < assembledqf_len; ++i)
+   {
+      *maxq = std::max(*maxq, tempdata[i]);
+      *minq = std::min(*minq, tempdata[i]);
+      *absmin = std::min(*absmin, fabs(tempdata[i]));
+   }
+   ierr = CeedVectorRestoreArrayRead(assembledqf, &tempdata); CeedChk(ierr);
+
+   ierr = CeedElemRestrictionDestroy(&rstr_q); CeedChk(ierr);
+   ierr = CeedVectorDestroy(&assembledqf); CeedChk(ierr);
+   return 0;
+}
+
+int CeedOperatorGetHeuristics(CeedOperator oper, CeedScalar* minq,
+                              CeedScalar* maxq, CeedScalar* absmin)
+{
+   int ierr;
+   bool isComposite;
+   ierr = CeedOperatorIsComposite(oper, &isComposite); CeedChk(ierr);
+   if (!isComposite)
+   {
+      return CeedSingleOperatorGetHeuristics(oper, minq, maxq, absmin);
+   }
+   
+   *minq = 1.e+12;
+   *maxq = -1.e+12;
+   *absmin = 1.e+12;
+
+   int nsub;
+   ierr = CeedOperatorGetNumSub(oper, &nsub); CeedChk(ierr);
+   CeedOperator *subops;
+   ierr = CeedOperatorGetSubList(oper, &subops); CeedChk(ierr);
+   for (int isub=0; isub<nsub; ++isub)
+   {
+      CeedOperator subop = subops[isub];
+      CeedScalar lminq, lmaxq, labsmin;
+      ierr = CeedSingleOperatorGetHeuristics(subop, &lminq, &lmaxq,
+                                             &labsmin); CeedChk(ierr);
+      *minq = std::min(lminq, *minq);
+      *maxq = std::max(lmaxq, *maxq);
+      *absmin = std::min(*absmin, labsmin);
+   }
+
+   return 0;
+}
+
+int CeedQFunctionQCoarsen(CeedOperator oper, CeedInt qorder_reduction,
+                          CeedVector* coarse_assembledqf,
+                          CeedElemRestriction* coarse_rstr_q, CeedBasis* qcoarse_basis,
+                          CeedQFunction* qfout, CeedQFunctionContext* context_ptr,
+                          CeedQuadMode fine_qmode, CeedQuadMode coarse_qmode)
+{
+   int ierr;
+   Ceed ceed;
+   ierr = CeedOperatorGetCeed(oper, &ceed); CeedChk(ierr);
+
+   CeedQFunction qfin;
+   ierr = CeedOperatorGetQFunction(oper, &qfin); CeedChk(ierr);
+   CeedInt vlength;
+   ierr = CeedQFunctionGetVectorLength(qfin, &vlength); CeedChk(ierr);
+
+   CeedVector assembledqf;
+   CeedElemRestriction rstr_q;
+   ierr = CeedOperatorLinearAssembleQFunction(
+      oper, &assembledqf, &rstr_q, CEED_REQUEST_IMMEDIATE); CeedChk(ierr);
+
+   CeedBasis qbasisctof; // P_Q
+   CeedBasis fine_basis;
+   CeedOperatorField active_field;
+   ierr = CeedOperatorGetActiveField(oper, &active_field); CeedChk(ierr);
+   ierr = CeedOperatorFieldGetBasis(active_field, &fine_basis); CeedChk(ierr);
+   ierr = CeedBasisQCoarsen(fine_basis, qcoarse_basis, &qbasisctof,
+                            qorder_reduction, fine_qmode, coarse_qmode); CeedChk(ierr);
+
+   int ncomp_rstr; // components in ElementRestriction rstr_q
+   ierr = CeedElementRestrictionQCoarsen(rstr_q, qbasisctof, coarse_rstr_q,
+                                         &ncomp_rstr); CeedChk(ierr);
+   CeedInt coarse_vlength = vlength; /// this looks wrong but is probably right
+
+   // coarsen: coarsen the vector assembledqf itself (using basisctof == P_Q)
+   ierr = CeedQFunctionCoarsenAssembledVector(assembledqf, rstr_q, qbasisctof,
+                                              coarse_assembledqf); CeedChk(ierr);
+
+   CeedInt coarse_layout[3];
+   ierr = CeedElemRestrictionGetELayout(*coarse_rstr_q, &coarse_layout); CeedChk(ierr);
+
+   /// hack_loc is a hack, in general try qcoarsen_linearfunc_loc
+   /// this runs on tuxbox CPU because the location isn't even queried except
+   /// with nvcc?
+   const char* hack_loc = "/usr/WS1/barker29/ceed-solvers/include/linear.h:qcoarsen_linearfunc";
+   ierr = CeedQFunctionCreateInterior(ceed, coarse_vlength, qcoarsen_linearfunc,
+                                      hack_loc, qfout); CeedChk(ierr);
+
+   struct LinearQFunctionContext * context = 
+      (struct LinearQFunctionContext *) calloc(1, sizeof(struct LinearQFunctionContext));
+   context->dim = -1;
+   context->ncomp = ncomp_rstr;
+   for (int i = 0; i < 3; ++i)
+   {
+      context->layout[i] = coarse_layout[i];
+   }
+   CeedQFunctionContext qf_context;
+   ierr = CeedQFunctionContextCreate(ceed, &qf_context); CeedChk(ierr);
+   ierr = CeedQFunctionContextSetData(qf_context, CEED_MEM_HOST, CEED_COPY_VALUES,
+                                      sizeof(*context), context); CeedChk(ierr);
+   ierr = CeedQFunctionSetContext(*qfout, qf_context); CeedChk(ierr);
+   *context_ptr = qf_context;
+   free(context);
+
+   ierr = CeedElemRestrictionDestroy(&rstr_q); CeedChk(ierr);
+   ierr = CeedBasisDestroy(&qbasisctof); CeedChk(ierr);
+   ierr = CeedVectorDestroy(&assembledqf); CeedChk(ierr);
+
+   return 0;
+}
+
+/** @brief given a CeedOperator, use the "assembled" qfunction to create your
+    own qfunction that has the same action (works for linear operators),
+    and build a new CeedOperator around that.
+
+    oper is in, qorder_reduction is in, everything else is out
+*/
+int CeedOperatorQCoarsen(CeedOperator oper, int qorder_reduction,
+                         CeedOperator* out, CeedVector* coarse_assembledqf,
+                         CeedQFunctionContext* context_ptr,
+                         CeedQuadMode fine_qmode, CeedQuadMode coarse_qmode)
+{
+   int ierr;
+   Ceed ceed;
+   ierr = CeedOperatorGetCeed(oper, &ceed); CeedChk(ierr);
+
+   CeedQFunction qfin;
+   ierr = CeedOperatorGetQFunction(oper, &qfin); CeedChk(ierr);
+   CeedElemRestriction coarse_rstr_q;
+   CeedQFunction qfout;
+   CeedBasis qcoarse_basis;
+   ierr = CeedQFunctionQCoarsen(oper, qorder_reduction, coarse_assembledqf,
+                                &coarse_rstr_q, &qcoarse_basis, &qfout,
+                                context_ptr,
+                                fine_qmode, coarse_qmode); CeedChk(ierr);
+
+   CeedInt numinputfields, numoutputfields;
+   ierr = CeedQFunctionGetNumArgs(qfin, &numinputfields, &numoutputfields); CeedChk(ierr);
+   CeedQFunctionField *inputqfields, *outputqfields;
+   ierr = CeedQFunctionGetFields(qfin, &inputqfields, &outputqfields); CeedChk(ierr);
+   CeedOperatorField *inputfields, *outputfields;
+   ierr = CeedOperatorGetFields(oper, &inputfields, &outputfields); CeedChk(ierr);
+
+   // Determine active input basis, get dimension, numemodein
+   CeedInt size;
+   char * fieldname;
+   CeedInt numemodein = 0, dim = 1;
+   CeedEvalMode emodein;
+   CeedBasis basisin = NULL;
+   CeedVector vec;
+   for (CeedInt i=0; i<numinputfields; i++)
+   {
+      ierr = CeedOperatorFieldGetVector(inputfields[i], &vec); CeedChk(ierr);
+      if (vec == CEED_VECTOR_ACTIVE)
+      {
+         ierr = CeedOperatorFieldGetBasis(inputfields[i], &basisin); CeedChk(ierr);
+         ierr = CeedBasisGetDimension(basisin, &dim); CeedChk(ierr);
+         ierr = CeedQFunctionFieldGetEvalMode(inputqfields[i], &emodein); CeedChk(ierr);
+         switch (emodein)
+         {
+         case CEED_EVAL_NONE:
+         case CEED_EVAL_INTERP:
+            numemodein += 1;
+            break;
+         case CEED_EVAL_GRAD:
+            numemodein += dim;
+            break;
+         case CEED_EVAL_WEIGHT:
+         case CEED_EVAL_DIV:
+         case CEED_EVAL_CURL:
+            break; // Caught by QF Assembly
+         }
+         ierr = CeedQFunctionFieldGetName(inputqfields[i], &fieldname); CeedChk(ierr);
+         ierr = CeedQFunctionFieldGetSize(inputqfields[i], &size); CeedChk(ierr);
+         ierr = CeedQFunctionAddInput(qfout, fieldname, size, emodein); CeedChk(ierr);
+      }
+      else
+      {
+         // don't do anything for inactive fields
+         // maybe count them to make sure there's exactly one?
+      }
+   }
+
+   ierr = CeedQFunctionAddInput(qfout, "assembled", numemodein*numemodein,
+                                CEED_EVAL_NONE); CeedChk(ierr);
+
+   // Determine active output basis, count emodeout
+   CeedInt numemodeout = 0;
+   CeedEvalMode emodeout;
+   for (CeedInt i=0; i<numoutputfields; i++)
+   {
+      ierr = CeedOperatorFieldGetVector(outputfields[i], &vec); CeedChk(ierr);
+      if (vec == CEED_VECTOR_ACTIVE)
+      {
+         ierr = CeedQFunctionFieldGetEvalMode(outputqfields[i], &emodeout);
+         CeedChk(ierr);
+         switch (emodeout)
+         {
+         case CEED_EVAL_NONE:
+         case CEED_EVAL_INTERP:
+            numemodeout += 1;
+            break;
+         case CEED_EVAL_GRAD:
+            numemodeout += dim;
+            break;
+         case CEED_EVAL_WEIGHT:
+         case CEED_EVAL_DIV:
+         case CEED_EVAL_CURL:
+            break; // Caught by QF Assembly
+         }
+         ierr = CeedQFunctionFieldGetName(outputqfields[i], &fieldname); CeedChk(ierr);
+         ierr = CeedQFunctionFieldGetSize(outputqfields[i], &size); CeedChk(ierr);
+         ierr = CeedQFunctionAddOutput(qfout, fieldname, size, emodeout); CeedChk(ierr);
+      }
+      else
+      {
+         // don't do anything for inactive fields
+      }
+   }
+
+   CeedOperator qcoper;
+   ierr = CeedOperatorCreate(ceed, qfout, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE,
+                             &qcoper); CeedChk(ierr);
+
+   bool linear_thing_set = false;
+   CeedBasis basis;
+   CeedElemRestriction er_input;
+   for (int i = 0; i < numinputfields; ++i)
+   {
+      ierr = CeedQFunctionFieldGetName(inputqfields[i], &fieldname); CeedChk(ierr);
+      ierr = CeedOperatorFieldGetVector(inputfields[i], &vec); CeedChk(ierr);
+      ierr = CeedOperatorFieldGetBasis(inputfields[i], &basis); CeedChk(ierr);
+      ierr = CeedOperatorFieldGetElemRestriction(inputfields[i], &er_input); CeedChk(ierr);
+      if (vec == CEED_VECTOR_ACTIVE)
+      {
+         ierr = CeedOperatorSetField(qcoper, fieldname, er_input, qcoarse_basis,
+                                     CEED_VECTOR_ACTIVE); CeedChk(ierr);
+      }
+      else
+      {
+         if (linear_thing_set)
+         {
+            return CeedError(ceed, 1, "Only know how to do one non-active vector!");
+         }
+         ierr = CeedOperatorSetField(qcoper, "assembled", coarse_rstr_q,
+                                     CEED_BASIS_COLLOCATED,
+                                     *coarse_assembledqf); CeedChk(ierr);
+         linear_thing_set = true;
+      }
+   }
+   if (!linear_thing_set)
+   {
+      return CeedError(ceed, 1, "Did not find active vector!");
+   }
+   CeedElemRestriction er_output;
+   for (int i = 0; i < numoutputfields; ++i)
+   {
+      ierr = CeedQFunctionFieldGetName(outputqfields[i], &fieldname); CeedChk(ierr);
+      ierr = CeedOperatorFieldGetVector(outputfields[i], &vec); CeedChk(ierr);
+      ierr = CeedOperatorFieldGetBasis(outputfields[i], &basis); CeedChk(ierr);
+      ierr = CeedOperatorFieldGetElemRestriction(outputfields[i], &er_output); CeedChk(ierr);
+      if (vec == CEED_VECTOR_ACTIVE)
+      {
+         ierr = CeedOperatorSetField(qcoper, fieldname, er_output, qcoarse_basis,
+                                     CEED_VECTOR_ACTIVE); CeedChk(ierr);
+      }
+      else
+      {
+         return CeedError(ceed, 1, "Don't think this should happen!");
+      }
+   }
+
+   // the following probably do not really get destroyed, but their refcounts
+   // get reduced, so now they are owned by qcoper
+   ierr = CeedElemRestrictionDestroy(&coarse_rstr_q); CeedChk(ierr);
+   ierr = CeedQFunctionDestroy(&qfout); CeedChk(ierr);
+   ierr = CeedBasisDestroy(&qcoarse_basis); CeedChk(ierr);
+
+   *out = qcoper;
+   return 0;
+}
+
+#endif // MFEM_USE_CEED

fem/libceed/ceedsolvers-qcoarsen.h

@@ -0,0 +1,33 @@
+// Copyright (c) 2010-2020, 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_CEEDSOLVERS_QCOARSEN_H
+#define MFEM_CEEDSOLVERS_QCOARSEN_H
+
+#include "../../config/config.hpp"
+
+#ifdef MFEM_USE_CEED
+
+#include <ceed.h>
+
+CEED_EXTERN int CeedOperatorQCoarsen(CeedOperator oper, int qorder_reduction,
+                                     CeedOperator* out, CeedVector* coarse_assembledqf,
+                                     CeedQFunctionContext* context_ptr,
+                                     CeedQuadMode fine_qmode, CeedQuadMode coarse_qmode);
+
+/// @todo eventually I want this per element or something?
+/// (also, this might be more expensive than we want)
+CEED_EXTERN int CeedOperatorGetHeuristics(CeedOperator oper, CeedScalar* minq,
+                                          CeedScalar* maxq, CeedScalar* absmin);
+
+#endif
+
+#endif

fem/libceed/ceedsolvers-sparsify.cpp

@@ -0,0 +1,450 @@
+// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "ceedsolvers-sparsify.h"
+#include "ceedsolvers-utility.h"
+
+#ifdef MFEM_USE_CEED
+#include <ceed-backend.h>
+
+#include <math.h>
+// todo: should probably use Ceed memory wrappers instead of calloc/free?
+#include <stdlib.h>
+
+typedef struct {
+  int index;
+  double basis_value;
+  double grad_value;
+} basis_element;
+
+/// for C (not C++) sorting using stdlib.h qsort()
+/// sort by basis_value
+int basis_element_compare_abs(const void *a, const void *b) {
+  basis_element *x = (basis_element *) a;
+  basis_element *y = (basis_element *) b;
+
+  // it is actually not so clear that absolute comparison
+  // is the right thing here
+  double xv = fabs(x->basis_value);
+  double yv = fabs(y->basis_value);
+  if (xv > yv) {
+    return  1;
+  }
+  if (xv < yv) {
+    return -1;
+  }
+  return 0;
+}
+
+/// sort by grad_value
+int basis_element_compare_grad_abs(const void *a, const void *b) {
+  basis_element *x = (basis_element *) a;
+  basis_element *y = (basis_element *) b;
+
+  double xv = fabs(x->grad_value);
+  double yv = fabs(y->grad_value);
+  if (xv > yv) {
+    return  1;
+  }
+  if (xv < yv) {
+    return -1;
+  }
+  return 0;
+}
+
+int basis_element_compare(const void *a, const void *b) {
+  basis_element *x = (basis_element *) a;
+  basis_element *y = (basis_element *) b;
+
+  double xv = x->basis_value;
+  double yv = y->basis_value;
+  if (xv > yv) {
+    return  1;
+  }
+  if (xv < yv) {
+    return -1;
+  }
+  return 0;
+}
+
+/// this guy puts the largest element of row at the end,
+/// and then the ones closest to it in *index* next to it
+/// at the end
+int near_largest(basis_element * row, int size) {
+  if (size < 4) return 0;
+
+  int largest_index = -1;
+  double largest_value = -1.e+10;
+  for (int i = 0; i < size; i++) {
+    double v = row[i].basis_value;
+    if (v > largest_value) {
+      largest_value = v;
+      largest_index = i;
+    }
+  }
+  int near_index_a = largest_index - 1;
+  if (near_index_a < 0) {
+    near_index_a = largest_index + 2;
+  }
+  int near_index_b = largest_index + 1;
+  if (near_index_b > size - 1) {
+    near_index_b = largest_index - 2;
+  }
+
+  basis_element temp = row[size - 1];
+  row[size - 1] = row[largest_index];
+  row[largest_index] = temp;
+
+  temp = row[size - 2];
+  row[size - 2] = row[near_index_a];
+  row[near_index_a] = temp;
+  
+  temp = row[size - 3];
+  row[size - 3] = row[near_index_b];
+  row[near_index_b] = temp;
+
+  return 0;
+}
+
+typedef struct {
+  double dist;
+  int index;
+} distindex;
+
+int compare_distindex(const void *a, const void *b) {
+  distindex *x = (distindex *) a;
+  distindex *y = (distindex *) b;
+
+  if (x->dist > y->dist) {
+    return  1;
+  }
+  if (x->dist < y->dist) {
+    return -1;
+  }
+  return 0;
+}
+
+/**
+   Replacement for near_largest, assume Gauss-Lobatto for nodes and
+   Gauss-Legendre for quad points, sorts by physical proximity
+
+   row has size P1d, qrefx is the coordinate of the quadpoint
+   corresponding to the row.
+*/
+int nearby_physical(basis_element * row, int P1d, double qrefx)
+{
+  int ierr;
+
+  CeedScalar * nodes = (CeedScalar*) calloc(P1d, sizeof(CeedScalar));
+  ierr = CeedLobattoQuadrature(P1d, nodes, NULL); CeedChk(ierr);
+  // so these nodes are in [-1, 1]
+  // but qrefx is in [0, 1], but why?
+
+  // I am just going to copy the row because I don't care about efficiency
+  basis_element old_row[P1d];
+  distindex distance_row[P1d];
+  for (int i = 0; i < P1d; ++i) {
+    const double node_loc = 0.5 * (nodes[i] + 1.0);
+    old_row[i] = row[i];
+    distance_row[i].dist = fabs(node_loc - qrefx);
+    distance_row[i].index = i;
+  }
+  qsort(distance_row, P1d, sizeof(distindex), compare_distindex);
+  for (int i = 0; i < P1d; ++i) {
+    // we want closest ones at the end, so reverse
+    // of the result from qsort
+    row[P1d - i - 1] = old_row[distance_row[i].index];
+  }
+  free(nodes);
+
+  return 0;
+}
+
+int CeedBasisSparsifyScaling(CeedBasis basisin, CeedBasis* basisout,
+                             SparsifySelectionStrategy sel_strategy,
+                             int parameter) {
+  int ierr;
+  Ceed ceed;
+  ierr = CeedBasisGetCeed(basisin, &ceed); CeedChk(ierr);
+  if (parameter < 0) {
+    return CeedError(ceed, 1, "Bad parameter value!");
+  }
+
+  CeedInt dim, ncomp, P1d, Q1d;
+  ierr = CeedBasisGetDimension(basisin, &dim); CeedChk(ierr);
+  ierr = CeedBasisGetNumComponents(basisin, &ncomp); CeedChk(ierr);
+  ierr = CeedBasisGetNumNodes1D(basisin, &P1d); CeedChk(ierr);
+  ierr = CeedBasisGetNumQuadraturePoints1D(basisin, &Q1d); CeedChk(ierr);
+
+  const CeedScalar *interp1d;
+  ierr = CeedBasisGetInterp1D(basisin, &interp1d); CeedChk(ierr);
+  const CeedScalar * grad1d;
+  ierr = CeedBasisGetGrad1D(basisin, &grad1d); CeedChk(ierr);
+  const CeedScalar * qref1d;
+  ierr = CeedBasisGetQRef(basisin, &qref1d); CeedChk(ierr);
+  const CeedScalar * qweight1d;
+  ierr = CeedBasisGetQWeights(basisin, &qweight1d); CeedChk(ierr);
+
+  CeedScalar * sparse_interp1d = (CeedScalar*) calloc(
+    P1d * Q1d, sizeof(CeedScalar));
+  CeedScalar * sparse_grad1d = (CeedScalar*) calloc(
+    P1d * Q1d, sizeof(CeedScalar));
+
+  basis_element * row = (basis_element*) calloc(P1d, sizeof(basis_element));
+  int sc = (P1d < parameter + 1) ? P1d : parameter + 1;
+    
+  for (int q = 0; q < Q1d; ++q) {
+    int basis_nonzeros = 0;
+    for (int i = 0; i < P1d; ++i) {
+      sparse_interp1d[q * P1d + i] = 0.0;
+      sparse_grad1d[q * P1d + i] = 0.0;
+
+      row[i].index = i;
+      row[i].basis_value = interp1d[q * P1d + i];
+      row[i].grad_value = grad1d[q * P1d + i];
+      if (fabs(row[i].basis_value) > 1.e-12) {
+        basis_nonzeros++;
+      }
+    }
+    if (sel_strategy == SPARSIFY_LARGEST_GRAD_ABS ||
+        basis_nonzeros < sc) {
+      qsort(row, P1d, sizeof(basis_element), basis_element_compare_grad_abs);
+    } else if (sel_strategy == SPARSIFY_LARGEST_ABS) {
+      qsort(row, P1d, sizeof(basis_element), basis_element_compare_abs);
+    } else if (sel_strategy == SPARSIFY_LARGEST_POSITIVE) {
+      qsort(row, P1d, sizeof(basis_element), basis_element_compare);
+    } else if (sel_strategy == SPARSIFY_NEARBY) {
+      // near_largest(row, P1d);
+      nearby_physical(row, P1d, qref1d[q]);
+    } else {
+      return CeedError(ceed, 1, "Bad strategy!");
+    }
+
+    double basis_sum = 0.0;
+    double grad_sum_plus = 0.0;
+    double grad_sum_minus = 0.0;
+    for (int k = 0; k < sc; k++) {
+      int reverse_index = P1d - 1 - k;
+      double basis_value = row[reverse_index].basis_value;
+      double grad_value = row[reverse_index].grad_value;
+      basis_sum += basis_value;
+      if (grad_value > 0.0) {
+        grad_sum_plus += grad_value;
+      } else {
+        grad_sum_minus += (-grad_value);
+      }
+    }
+    double grad_scale_plus = (grad_sum_plus + grad_sum_minus) /
+      (2.0 * grad_sum_plus);
+    if (grad_sum_plus == 0.0) {
+      grad_scale_plus = 1.0;
+    }
+    double grad_scale_minus = (grad_sum_plus + grad_sum_minus) /
+      (2.0 * grad_sum_minus);
+    if (grad_sum_minus == 0.0) {
+      grad_scale_minus = 1.0;
+    }
+    for (int k = 0; k < sc; k++) {
+      int reverse_index = P1d - 1 - k;
+      int i = row[reverse_index].index;
+      double basis_value = row[reverse_index].basis_value / basis_sum;
+      double grad_value = row[reverse_index].grad_value;
+      sparse_interp1d[q * P1d + i] = basis_value;
+      if (grad_value > 0.0) {
+        sparse_grad1d[q * P1d + i] = grad_value * grad_scale_plus;
+      } else {
+        sparse_grad1d[q * P1d + i] = grad_value * grad_scale_minus;
+      }
+    }
+  }
+  free(row);
+
+  ierr = CeedBasisCreateTensorH1(ceed, dim, ncomp,
+                                 P1d, Q1d, sparse_interp1d, sparse_grad1d,
+                                 qref1d, qweight1d, basisout); CeedChk(ierr);
+
+  const bool debug = false;
+  if (debug) {
+    double frob = 0.0;
+    for (int q = 0; q < Q1d; q++) {
+      for (int i = 0; i < P1d; i++) {
+        printf("  interp1d[%d,%d] = %f, sparse = %f\n",
+               q, i, interp1d[q * P1d + i], sparse_interp1d[q * P1d + i]);
+        frob += (interp1d[q * P1d + i] -  sparse_interp1d[q * P1d + i]) *
+          (interp1d[q * P1d + i] -  sparse_interp1d[q * P1d + i]);
+      }
+    }
+    printf("interp1d frob diff = %e\n", sqrt(frob));
+    frob = 0.0;
+    for (int q = 0; q < Q1d; q++) {
+      for (int i = 0; i < P1d; i++) {
+        printf("  grad1d[%d,%d] = %f, sparse = %f\n",
+               q, i, grad1d[q * P1d + i], sparse_grad1d[q * P1d + i]);
+        frob += (grad1d[q * P1d + i] -  sparse_grad1d[q * P1d + i]) *
+          (grad1d[q * P1d + i] -  sparse_grad1d[q * P1d + i]);
+      }
+    }
+    printf("grad1d frob diff = %e\n", sqrt(frob));
+  }
+
+  free(sparse_interp1d);
+  free(sparse_grad1d);
+
+  return 0;
+}
+
+int CeedSparsifySingleOperator(CeedOperator oper, int sparse_parameter,
+                               SparsifySelectionStrategy sel_strategy,
+                               int (*basis_sparsify)(CeedBasis, CeedBasis*, SparsifySelectionStrategy, int),
+                               CeedBasis* sparse_basis_out,
+                               CeedOperator* out) {
+  int ierr;
+  Ceed ceed;
+  ierr = CeedOperatorGetCeed(oper, &ceed); CeedChk(ierr);
+
+  CeedQFunction qf;
+  bool isComposite;
+  ierr = CeedOperatorIsComposite(oper, &isComposite); CeedChk(ierr);
+  if (isComposite) {
+     return CeedError(ceed, 1, "This function does not make sense for composite operator!");
+  }
+  ierr = CeedOperatorGetQFunction(oper, &qf); CeedChk(ierr);
+  CeedInt numinputfields, numoutputfields;
+  ierr = CeedQFunctionGetNumArgs(qf, &numinputfields, &numoutputfields);
+  CeedQFunctionField *inputqfields, *outputqfields;
+  ierr = CeedQFunctionGetFields(qf, &inputqfields, &outputqfields); CeedChk(ierr);
+  CeedOperatorField *inputfields, *outputfields;
+  ierr = CeedOperatorGetFields(oper, &inputfields, &outputfields); CeedChk(ierr);
+
+  CeedElemRestriction * er_input =
+    (CeedElemRestriction*) calloc(numinputfields, sizeof(CeedElemRestriction));
+  CeedElemRestriction * er_output =
+    (CeedElemRestriction*) calloc(numoutputfields, sizeof(CeedElemRestriction));
+  CeedVector * if_vector = 
+    (CeedVector*) calloc(numinputfields, sizeof(CeedVector));
+  CeedVector * of_vector = 
+    (CeedVector*) calloc(numoutputfields, sizeof(CeedVector));
+  CeedBasis * basis_input =
+    (CeedBasis*) calloc(numinputfields, sizeof(CeedBasis));
+  CeedBasis * basis_output = 
+    (CeedBasis*) calloc(numoutputfields, sizeof(CeedBasis));
+  CeedBasis sbasis;
+  int active_input_basis = -1;
+  for (int i = 0; i < numinputfields; ++i) {
+    ierr = CeedOperatorFieldGetElemRestriction(inputfields[i],
+                                               &er_input[i]); CeedChk(ierr);
+    ierr = CeedOperatorFieldGetVector(inputfields[i], &if_vector[i]); CeedChk(ierr);
+    ierr = CeedOperatorFieldGetBasis(inputfields[i], &basis_input[i]); CeedChk(ierr);
+    if (if_vector[i] == CEED_VECTOR_ACTIVE) {
+      if (active_input_basis < 0) {
+        // sparsify!
+        ierr = (*basis_sparsify)(basis_input[i], &sbasis, sel_strategy, sparse_parameter); CeedChk(ierr);
+        active_input_basis = i;
+      } else {
+        // should already be sparsified
+        if (basis_input[i] != basis_input[active_input_basis]) {
+          return CeedError(ceed, 1, "Two different active input basis!");
+        }
+      }
+    }
+  }
+  for (int i = 0; i < numoutputfields; ++i) {
+    ierr = CeedOperatorFieldGetElemRestriction(outputfields[i],
+                                               &er_output[i]); CeedChk(ierr);
+    ierr = CeedOperatorFieldGetVector(outputfields[i], &of_vector[i]); CeedChk(ierr);
+    ierr = CeedOperatorFieldGetBasis(outputfields[i], &basis_output[i]); CeedChk(ierr);
+    if (of_vector[i] == CEED_VECTOR_ACTIVE) {
+      // should already be sparsified
+      if (basis_output[i] != basis_input[active_input_basis]) {
+        return CeedError(ceed, 1, "Input and output basis do not match!");
+      }
+    }
+  }
+
+  CeedOperator soper;
+  ierr = CeedOperatorCreate(ceed, qf, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE,
+                            &soper); CeedChk(ierr);
+
+  for (int i = 0; i < numinputfields; ++i) {
+    char * fieldname;
+    ierr = CeedQFunctionFieldGetName(inputqfields[i], &fieldname); CeedChk(ierr);
+    if (if_vector[i] == CEED_VECTOR_ACTIVE) {
+      ierr = CeedOperatorSetField(soper, fieldname, er_input[i], sbasis, 
+                                  if_vector[i]); CeedChk(ierr);
+    } else {
+      ierr = CeedOperatorSetField(soper, fieldname, er_input[i], basis_input[i], 
+                                  if_vector[i]); CeedChk(ierr);
+    }
+  }
+  for (int i = 0; i < numoutputfields; ++i) {
+    char * fieldname;
+    ierr = CeedQFunctionFieldGetName(outputqfields[i], &fieldname); CeedChk(ierr);
+    if (of_vector[i] == CEED_VECTOR_ACTIVE) {
+      ierr = CeedOperatorSetField(soper, fieldname, er_output[i], sbasis,
+                                  of_vector[i]); CeedChk(ierr);
+    } else {
+      ierr = CeedOperatorSetField(soper, fieldname, er_output[i], basis_output[i],
+                                  of_vector[i]); CeedChk(ierr);
+    }
+  }
+  free(er_input);
+  free(er_output);
+  free(if_vector);
+  free(of_vector);
+  free(basis_input);
+  free(basis_output);
+
+  *sparse_basis_out = sbasis;
+  *out = soper;
+  return 0;
+}
+
+int CeedSparsifyH1Operator(CeedOperator oper, int sparse_parameter,
+                           SparsifySelectionStrategy sel_strategy,
+                           int (*basis_sparsify)(CeedBasis, CeedBasis*, SparsifySelectionStrategy, int),
+                           CeedBasis* sparse_basis_out,
+                           CeedOperator* out) {
+   int ierr;
+   Ceed ceed;
+   ierr = CeedOperatorGetCeed(oper, &ceed); CeedChk(ierr);
+
+   CeedInt numsub;
+   CeedOperator *subops;
+   bool isComposite;
+   ierr = CeedOperatorIsComposite(oper, &isComposite); CeedChk(ierr);
+   if (isComposite) {
+      ierr = CeedOperatorGetNumSub(oper, &numsub); CeedChk(ierr);
+      ierr = CeedOperatorGetSubList(oper, &subops); CeedChk(ierr);
+      if (numsub != 1) {
+         // implementing this case is pretty easy but I am lazy
+         return CeedError(ceed, 1, "Not implemented for multiple suboperators!");
+      }
+      ierr = CeedSparsifySingleOperator(subops[0], sparse_parameter,
+                                        sel_strategy, basis_sparsify,
+                                        sparse_basis_out, out);
+   } else {
+      ierr = CeedSparsifySingleOperator(oper,  sparse_parameter,
+                                        sel_strategy, basis_sparsify,
+                                        sparse_basis_out, out);
+   }      
+   return 0;
+}
+
+int CeedSparsifySimple(CeedOperator oper, CeedBasis* sparse_basis_out,
+                       CeedOperator* out) {
+   int ierr;
+   ierr = CeedSparsifyH1Operator(oper, 1, SPARSIFY_LARGEST_POSITIVE,
+                                 CeedBasisSparsifyScaling, sparse_basis_out,
+                                 out); CeedChk(ierr);
+   return 0;
+}
+
+#endif

fem/libceed/ceedsolvers-sparsify.h

@@ -0,0 +1,70 @@
+// Copyright (c) 2010-2020, 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_CEEDSOLVERS_SPARSIFY_H
+#define MFEM_CEEDSOLVERS_SPARSIFY_H
+
+#include "../../config/config.hpp"
+
+#ifdef MFEM_USE_CEED
+
+#include <ceed.h>
+
+typedef enum {
+  SPARSIFY_LARGEST_ABS = 0,
+  SPARSIFY_LARGEST_POSITIVE = 1,
+  SPARSIFY_NEARBY = 2,
+  SPARSIFY_LARGEST_GRAD_ABS = 3,
+} SparsifySelectionStrategy;
+
+/**
+   For each (fixed) quadrature point, you sort the dof indices
+   based on the value of the basis function corresponding to
+   that dof. Then you zero the basis function and the gradient
+   for all but the largest <parameter> values.
+
+   basisin must be a tensor basis, ie, have 1d tensors available
+
+   scaling is to ensure that rows of interp1d have rowsum == 1,
+   and to ensure that rows of grad1d have rowsum == 0
+*/
+int CeedBasisSparsifyScaling(CeedBasis basisin, CeedBasis* basisout,
+                             SparsifySelectionStrategy sel_strategy,
+                             int parameter);
+
+/**
+   This is all modeled after H1 operators, do not expect it to work
+   on anything else. You also must have some symmetry in the basis
+   (input/output)
+
+   The only original thing happens in (basis_sparsify); this tries to
+   just copy the pointers/data/etc for everything else.
+
+   Note well that if you ask me to sparsify a different operator, mass
+   or something else, this will all fall apart.
+
+   Caller is responsible for deleting sparse_basis_out and out
+   out is what you actually want to use, I only return sparse_basis_out
+   so you can delete it cleanly.
+*/
+int CeedSparsifyH1Operator(CeedOperator oper, int sparse_parameter,
+                           SparsifySelectionStrategy sel_strategy,
+                           int (*basis_sparsify)(CeedBasis, CeedBasis*, SparsifySelectionStrategy, int),
+                           CeedBasis* sparse_basis_out,
+                           CeedOperator* out);
+
+int CeedSparsifySimple(CeedOperator oper, CeedBasis* sparse_basis_out,
+                       CeedOperator* out);
+
+#endif // MFEM_USE_CEED
+
+#endif // include guard

fem/libceed/ceedsolvers-utility.cpp

@@ -0,0 +1,379 @@
+// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "ceedsolvers-utility.h"
+
+#include "../../general/forall.hpp"
+using namespace mfem;
+
+#ifdef MFEM_USE_CEED
+
+/// a = a (pointwise*) b
+/// @todo: using MPI_FORALL in this Ceed-like function is ugly
+int CeedVectorPointwiseMult(CeedVector a, const CeedVector b)
+{
+   int ierr;
+   Ceed ceed;
+   CeedVectorGetCeed(a, &ceed);
+
+   int length, length2;
+   ierr = CeedVectorGetLength(a, &length); CeedChk(ierr);
+   ierr = CeedVectorGetLength(b, &length2); CeedChk(ierr);
+   if (length != length2)
+   {
+      return CeedError(ceed, 1, "Vector sizes don't match");
+   }
+
+   CeedMemType mem;
+   if (Device::Allows(Backend::DEVICE_MASK))
+   {
+      mem = CEED_MEM_DEVICE;
+   }
+   else
+   {
+      mem = CEED_MEM_HOST;
+   }
+   CeedScalar *a_data;
+   const CeedScalar *b_data;
+   ierr = CeedVectorGetArray(a, mem, &a_data); CeedChk(ierr);
+   ierr = CeedVectorGetArrayRead(b, mem, &b_data); CeedChk(ierr);
+   MFEM_FORALL(i, length,
+   {a_data[i] *= b_data[i];});
+
+   ierr = CeedVectorRestoreArray(a, &a_data); CeedChk(ierr);
+   ierr = CeedVectorRestoreArrayRead(b, &b_data); CeedChk(ierr);
+
+   return 0;
+}
+
+/// assumes a tensor-product operator with one active field
+int CeedOperatorGetActiveField(CeedOperator oper, CeedOperatorField *field)
+{
+   int ierr;
+   Ceed ceed;
+   ierr = CeedOperatorGetCeed(oper, &ceed); CeedChk(ierr);
+
+   CeedQFunction qf;
+   bool isComposite;
+   ierr = CeedOperatorIsComposite(oper, &isComposite); CeedChk(ierr);
+   CeedOperator *subops;
+   if (isComposite)
+   {
+      ierr = CeedOperatorGetSubList(oper, &subops); CeedChk(ierr);
+      ierr = CeedOperatorGetQFunction(subops[0], &qf); CeedChk(ierr);
+   }
+   else
+   {
+      ierr = CeedOperatorGetQFunction(oper, &qf); CeedChk(ierr);
+   }
+   CeedInt numinputfields, numoutputfields;
+   ierr = CeedQFunctionGetNumArgs(qf, &numinputfields, &numoutputfields);
+   CeedOperatorField *inputfields;
+   if (isComposite)
+   {
+      ierr = CeedOperatorGetFields(subops[0], &inputfields, NULL); CeedChk(ierr);
+   }
+   else
+   {
+      ierr = CeedOperatorGetFields(oper, &inputfields, NULL); CeedChk(ierr);
+   }
+
+   CeedVector if_vector;
+   bool found = false;
+   int found_index = -1;
+   for (int i = 0; i < numinputfields; ++i)
+   {
+      ierr = CeedOperatorFieldGetVector(inputfields[i], &if_vector); CeedChk(ierr);
+      if (if_vector == CEED_VECTOR_ACTIVE)
+      {
+         if (found)
+         {
+            return CeedError(ceed, 1, "Multiple active vectors in CeedOperator!");
+         }
+         found = true;
+         found_index = i;
+      }
+   }
+   if (!found)
+   {
+      return CeedError(ceed, 1, "No active vector in CeedOperator!");
+   }
+   *field = inputfields[found_index];
+
+   return 0;
+}
+
+/// (a better design splits this into CeedOperatorGetActiveBasis() and
+/// CeedOperatorGetOrder, which calls the basis one)
+/// TODO: unit test
+int CeedOperatorGetOrder(CeedOperator oper, CeedInt * order)
+{
+   int ierr;
+
+   CeedOperatorField active_field;
+   ierr = CeedOperatorGetActiveField(oper, &active_field); CeedChk(ierr);
+   CeedBasis basis;
+   ierr = CeedOperatorFieldGetBasis(active_field, &basis); CeedChk(ierr);
+   int P1d;
+   ierr = CeedBasisGetNumNodes1D(basis, &P1d); CeedChk(ierr);
+   *order = P1d - 1;
+
+   return 0;
+}
+
+int CeedOperatorGetActiveElemRestriction(CeedOperator oper,
+                                         CeedElemRestriction* restr_out)
+{
+   int ierr;
+
+   CeedOperatorField active_field;
+   ierr = CeedOperatorGetActiveField(oper, &active_field); CeedChk(ierr);
+   CeedElemRestriction er;
+   ierr = CeedOperatorFieldGetElemRestriction(active_field, &er); CeedChk(ierr);
+   *restr_out = er;
+
+   return 0;
+}
+
+/// assumes a square operator (you could do rectangular, you'd have
+/// to find separate active input and output fields/restrictions)
+int CeedOperatorGetSize(CeedOperator oper, CeedInt * size)
+{
+   int ierr;
+   CeedElemRestriction er;
+   ierr = CeedOperatorGetActiveElemRestriction(oper, &er); CeedChk(ierr);
+   ierr = CeedElemRestrictionGetLVectorSize(er, size); CeedChk(ierr);
+   return 0;
+}
+
+/** Just like CeedBasisCreateTensorH1Lagrange but with Legendre dofs
+    (basically only one line changed, which seems wasteful) */
+int CeedBasisCreateTensorH1Gauss(Ceed ceed, CeedInt dim, CeedInt ncomp,
+                                 CeedInt P, CeedInt Q, CeedQuadMode qmode,
+                                 CeedBasis *basis) 
+{
+   // Allocate
+   int ierr, i, j, k;
+   CeedScalar c1, c2, c3, c4, dx, *nodes, *interp1d, *grad1d, 
+      *qref1d, *qweight1d, *dummyweights;
+
+   if (dim<1)
+   {
+      return CeedError(ceed, 1, "Basis dimension must be a positive value");
+   }
+
+   // CeedCalloc replaced below
+   interp1d = (CeedScalar*) calloc(P * Q, sizeof(CeedScalar));
+   grad1d = (CeedScalar*) calloc(P * Q, sizeof(CeedScalar));
+   nodes = (CeedScalar*) calloc(P, sizeof(CeedScalar));
+   dummyweights = (CeedScalar*) calloc(P, sizeof(CeedScalar));
+   qref1d = (CeedScalar*) calloc(Q, sizeof(CeedScalar));
+   qweight1d = (CeedScalar*) calloc(Q, sizeof(CeedScalar));
+
+   // Get Nodes and Weights
+   ierr = CeedGaussQuadrature(P, nodes, dummyweights); CeedChk(ierr);
+   switch (qmode)
+   {
+   case CEED_GAUSS:
+      ierr = CeedGaussQuadrature(Q, qref1d, qweight1d); CeedChk(ierr);
+      break;
+   case CEED_GAUSS_LOBATTO:
+      ierr = CeedLobattoQuadrature(Q, qref1d, qweight1d); CeedChk(ierr);
+      break;
+   }
+   // Build B, D matrix
+   // Fornberg, 1998
+   for (i = 0; i  < Q; i++)
+   {
+      c1 = 1.0;
+      c3 = nodes[0] - qref1d[i];
+      interp1d[i*P+0] = 1.0;
+      for (j = 1; j < P; j++)
+      {
+         c2 = 1.0;
+         c4 = c3;
+         c3 = nodes[j] - qref1d[i];
+         for (k = 0; k < j; k++)
+         {
+            dx = nodes[j] - nodes[k];
+            c2 *= dx;
+            if (k == j - 1)
+            {
+               grad1d[i*P + j] = c1*(interp1d[i*P + k] - c4*grad1d[i*P + k]) / c2;
+               interp1d[i*P + j] = - c1*c4*interp1d[i*P + k] / c2;
+            }
+            grad1d[i*P + k] = (c3*grad1d[i*P + k] - interp1d[i*P + k]) / dx;
+            interp1d[i*P + k] = c3*interp1d[i*P + k] / dx;
+         }
+         c1 = c2;
+      }
+   }
+   //  // Pass to CeedBasisCreateTensorH1
+   ierr = CeedBasisCreateTensorH1(ceed, dim, ncomp, P, Q, interp1d, grad1d, qref1d,
+                                  qweight1d, basis); CeedChk(ierr);
+   /*
+     ierr = CeedFree(&interp1d); CeedChk(ierr);
+     ierr = CeedFree(&grad1d); CeedChk(ierr);
+     ierr = CeedFree(&nodes); CeedChk(ierr);
+     ierr = CeedFree(&dummyweights); CeedChk(ierr);
+     ierr = CeedFree(&qref1d); CeedChk(ierr);
+     ierr = CeedFree(&qweight1d); CeedChk(ierr);
+   */
+   free(interp1d);
+   free(grad1d);
+   free(nodes);
+   free(dummyweights);
+   free(qref1d);
+   free(qweight1d);
+   return 0;
+}
+
+/** Ugly hacky copy/paste from CeedBasisCreateTensorH1Lagrange to deal with
+    different reference elements for MFEM/Ceed.
+
+    The correct way to do this is to actually interpolate, but I am lazy. */
+int CeedBasisCreateMFEMTensorH1Lagrange(Ceed ceed, CeedInt dim, CeedInt ncomp,
+                                        CeedInt P, CeedInt Q, CeedQuadMode qmode,
+                                        CeedBasis *basis)
+{
+   // Allocate
+   int ierr, i, j, k;
+   CeedScalar c1, c2, c3, c4, dx, *nodes, *interp1d, *grad1d, *qref1d, *qweight1d;
+
+   if (dim<1)
+   {
+      return CeedError(ceed, 1, "Basis dimension must be a positive value");
+   }
+
+   interp1d = (CeedScalar*) calloc(P*Q, sizeof(CeedScalar));
+   grad1d = (CeedScalar*) calloc(P*Q, sizeof(CeedScalar));
+   nodes = (CeedScalar*) calloc(P, sizeof(CeedScalar));
+   qref1d = (CeedScalar*) calloc(Q, sizeof(CeedScalar));
+   qweight1d = (CeedScalar*) calloc(Q, sizeof(CeedScalar));
+
+   // Get Nodes and Weights
+   ierr = CeedLobattoQuadrature(P, nodes, NULL); CeedChk(ierr);
+   switch (qmode)
+   {
+   case CEED_GAUSS:
+      ierr = CeedGaussQuadrature(Q, qref1d, qweight1d); CeedChk(ierr);
+      break;
+   case CEED_GAUSS_LOBATTO:
+      ierr = CeedLobattoQuadrature(Q, qref1d, qweight1d); CeedChk(ierr);
+      break;
+   }
+
+   /// modification for MFEM reference element
+   for (int j = 0; j < P; ++j)
+   {
+      nodes[j] = 0.5 + 0.5*nodes[j];
+   }
+   for (int q = 0; q < Q; ++q)
+   {
+      qref1d[q] = 0.5 + 0.5*qref1d[q];
+      qweight1d[q] *= 0.5;
+   }
+
+   // Build B, D matrix
+   // Fornberg, 1998
+   for (i = 0; i < Q; i++)
+   {
+      c1 = 1.0;
+      c3 = nodes[0] - qref1d[i];
+      interp1d[i*P+0] = 1.0;
+      for (j = 1; j < P; j++)
+      {
+         c2 = 1.0;
+         c4 = c3;
+         c3 = nodes[j] - qref1d[i];
+         for (k = 0; k < j; k++)
+         {
+            dx = nodes[j] - nodes[k];
+            c2 *= dx;
+            if (k == j - 1)
+            {
+               grad1d[i*P + j] = c1*(interp1d[i*P + k] - c4*grad1d[i*P + k]) / c2;
+               interp1d[i*P + j] = - c1*c4*interp1d[i*P + k] / c2;
+            }
+            grad1d[i*P + k] = (c3*grad1d[i*P + k] - interp1d[i*P + k]) / dx;
+            interp1d[i*P + k] = c3*interp1d[i*P + k] / dx;
+         }
+         c1 = c2;
+      }
+   }
+   //  // Pass to CeedBasisCreateTensorH1
+   ierr = CeedBasisCreateTensorH1(ceed, dim, ncomp, P, Q, interp1d, grad1d, qref1d,
+                                  qweight1d, basis); CeedChk(ierr);
+  
+   /*
+     ierr = CeedFree(&interp1d); CeedChk(ierr);
+     ierr = CeedFree(&grad1d); CeedChk(ierr);
+     ierr = CeedFree(&nodes); CeedChk(ierr);
+     ierr = CeedFree(&qref1d); CeedChk(ierr);
+     ierr = CeedFree(&qweight1d); CeedChk(ierr);
+   */
+   free(interp1d);
+   free(grad1d);
+   free(nodes);
+   free(qref1d);
+   free(qweight1d);
+   return 0;
+}
+
+int CeedOperatorGetActiveBasis(CeedOperator oper, CeedBasis *basis)
+{
+   int ierr;
+   Ceed ceed;
+   ierr = CeedOperatorGetCeed(oper, &ceed); CeedChk(ierr);
+   CeedQFunction qf;
+   ierr = CeedOperatorGetQFunction(oper, &qf); CeedChk(ierr);
+   CeedInt numinputfields, numoutputfields;
+   ierr = CeedQFunctionGetNumArgs(qf, &numinputfields, &numoutputfields);
+   CeedOperatorField *inputfields, *outputfields;
+   ierr = CeedOperatorGetFields(oper, &inputfields, &outputfields); CeedChk(ierr);
+
+   *basis = NULL;
+   for (int i = 0; i < numinputfields; ++i)
+   {
+      CeedVector if_vector;
+      CeedBasis basis_in;
+      ierr = CeedOperatorFieldGetVector(inputfields[i], &if_vector); CeedChk(ierr);
+      ierr = CeedOperatorFieldGetBasis(inputfields[i], &basis_in); CeedChk(ierr);
+      if (if_vector == CEED_VECTOR_ACTIVE)
+      {
+         if (*basis == NULL)
+         {
+            *basis = basis_in;
+         }
+         else if (*basis != basis_in)
+         {
+            return CeedError(ceed, 1, "Two different active input basis!");
+         }
+      }
+   }
+   for (int i = 0; i < numoutputfields; ++i)
+   {
+      CeedVector of_vector;
+      CeedBasis basis_out;
+      ierr = CeedOperatorFieldGetVector(outputfields[i], &of_vector); CeedChk(ierr);
+      ierr = CeedOperatorFieldGetBasis(outputfields[i], &basis_out); CeedChk(ierr);
+      if (of_vector == CEED_VECTOR_ACTIVE)
+      {
+         if (*basis != basis_out)
+         {
+            return CeedError(ceed, 1, "Input and output basis do not match!");
+         }
+      }
+   }
+   return 0;
+}
+
+#endif // MFEM_USE_CEED

fem/libceed/ceedsolvers-utility.h

@@ -0,0 +1,38 @@
+// Copyright (c) 2010-2020, 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 __CEEDUTILITY_H
+#define __CEEDUTILITY_H
+
+#include "../../config/config.hpp"
+
+#ifdef MFEM_USE_CEED
+
+#include <ceed.h>
+#include <ceed-backend.h>
+
+int CeedVectorPointwiseMult(CeedVector a, const CeedVector b);
+int CeedOperatorGetActiveField(CeedOperator oper, CeedOperatorField *field);
+int CeedOperatorGetOrder(CeedOperator oper, CeedInt * order);
+int CeedOperatorGetActiveElemRestriction(CeedOperator oper,
+                                         CeedElemRestriction* restr_out);
+int CeedOperatorGetSize(CeedOperator oper, CeedInt * size);
+int CeedBasisCreateTensorH1Gauss(Ceed ceed, CeedInt dim, CeedInt ncomp,
+                                 CeedInt P, CeedInt Q, CeedQuadMode qmode,
+                                 CeedBasis *basis);
+int CeedBasisCreateMFEMTensorH1Lagrange(Ceed ceed, CeedInt dim, CeedInt ncomp,
+                                        CeedInt P, CeedInt Q, CeedQuadMode qmode,
+                                        CeedBasis *basis);
+int CeedOperatorGetActiveBasis(CeedOperator oper, CeedBasis *basis);
+
+#endif // MFEM_USE_CEED
+
+#endif // include guard

fem/libceed/linear.h

@@ -0,0 +1,56 @@
+
+#ifndef QCOARSEN_LINEAR_H
+#define QCOARSEN_LINEAR_H
+
+struct LinearQFunctionContext { 
+  CeedInt dim;
+  CeedInt ncomp; // beginnning to think this is only thing that matters
+  CeedInt layout[3];
+};
+
+// the idea here is the context has whatever comes out
+// of CeedOperatorLinearAssembleQFunction, and then we can
+// automatically apply it here
+CEED_QFUNCTION(qcoarsen_linearfunc)(void *ctx, const CeedInt Q,
+                                    const CeedScalar *const *in,
+                                    CeedScalar *const *out) {
+  const CeedScalar *ug = in[0];
+  const CeedScalar *linear_data = in[1];
+  CeedScalar *vg = out[0];
+
+  struct LinearQFunctionContext *context = (struct LinearQFunctionContext *) ctx;
+
+  const int ncomp = context->ncomp;
+  // const int elemsize = context->layout[1];
+  // const int lsize = context->layout[2];
+
+  switch (ncomp) {
+  case 1:
+    for (CeedInt i=0; i<Q; i++) {
+      vg[i] = ug[i] * linear_data[i];
+    }
+    break;
+  case 4:
+    for (CeedInt i=0; i<Q; i++) {
+      const CeedScalar ug0 = ug[i+Q*0];
+      const CeedScalar ug1 = ug[i+Q*1];
+      vg[i+Q*0] = linear_data[i+Q*0]*ug0 + linear_data[i+Q*2]*ug1;
+      vg[i+Q*1] = linear_data[i+Q*1]*ug0 + linear_data[i+Q*3]*ug1;
+    }
+    break;
+  case 9:
+    for (CeedInt i=0; i<Q; i++) {
+      const CeedScalar ug0 = ug[i+Q*0];
+      const CeedScalar ug1 = ug[i+Q*1];
+      const CeedScalar ug2 = ug[i+Q*2];
+      vg[i+Q*0] = linear_data[i+Q*0]*ug0 + linear_data[i+Q*3]*ug1 + linear_data[i+Q*6]*ug2;
+      vg[i+Q*1] = linear_data[i+Q*1]*ug0 + linear_data[i+Q*4]*ug1 + linear_data[i+Q*7]*ug2;
+      vg[i+Q*2] = linear_data[i+Q*2]*ug0 + linear_data[i+Q*5]*ug1 + linear_data[i+Q*8]*ug2;
+    }
+    break;
+  }
+
+  return 0;
+}
+
+#endif

fem/pfespace.cpp

@@ -3010,6 +3010,43 @@ void ParFiniteElementSpace::Update(bool want_transform)
    }
 }
 
+ConformingProlongationOperator::ConformingProlongationOperator(
+   int lsize, const GroupCommunicator &gc_, bool local_)
+   : gc(gc_),
+     local(local_)
+{
+   const Table &group_ldof = gc.GroupLDofTable();
+
+   int n_external = 0;
+   for (int g=1; g<group_ldof.Size(); ++g)
+   {
+      if (!gc.GetGroupTopology().IAmMaster(g))
+      {
+         n_external += group_ldof.RowSize(g);
+      }
+   }
+   int tsize = lsize - n_external;
+
+   height = lsize;
+   width = tsize;
+
+   external_ldofs.Reserve(n_external);
+   for (int gr = 1; gr < group_ldof.Size(); gr++)
+   {
+      if (!gc.GetGroupTopology().IAmMaster(gr))
+      {
+         external_ldofs.Append(group_ldof.GetRow(gr), group_ldof.RowSize(gr));
+      }
+   }
+   external_ldofs.Sort();
+}
+
+const GroupCommunicator &ConformingProlongationOperator::GetGroupCommunicator()
+const
+{
+   return gc;
+}
+
 ConformingProlongationOperator::ConformingProlongationOperator(
    const ParFiniteElementSpace &pfes, bool local_)
    : Operator(pfes.GetVSize(), pfes.GetTrueVSize()),
@@ -3121,17 +3158,13 @@ void ConformingProlongationOperator::MultTranspose(
 }
 
 DeviceConformingProlongationOperator::DeviceConformingProlongationOperator(
-   const ParFiniteElementSpace &pfes,
-   bool local_) :
-   ConformingProlongationOperator(pfes),
-   mpi_gpu_aware(Device::GetGPUAwareMPI()),
-   local(local_)
+   const GroupCommunicator &gc_, const SparseMatrix *R, bool local_)
+   : ConformingProlongationOperator(R->Width(), gc_),
+     mpi_gpu_aware(Device::GetGPUAwareMPI()),
+     local(local_)
 {
-   MFEM_ASSERT(pfes.Conforming(), "internal error");
-   const SparseMatrix *R = pfes.GetRestrictionMatrix();
    MFEM_ASSERT(R->Finalized(), "");
    const int tdofs = R->Height();
-   MFEM_ASSERT(tdofs == pfes.GetTrueVSize(), "");
    MFEM_ASSERT(tdofs == R->HostReadI()[tdofs], "");
    ltdof_ldof = Array<int>(const_cast<int*>(R->HostReadJ()), tdofs);
    ltdof_ldof.UseDevice();
@@ -3191,6 +3224,16 @@ DeviceConformingProlongationOperator::DeviceConformingProlongationOperator(
    requests = new MPI_Request[req_counter];
 }
 
+DeviceConformingProlongationOperator::DeviceConformingProlongationOperator(
+   const ParFiniteElementSpace &pfes, bool local_)
+   : DeviceConformingProlongationOperator(pfes.GroupComm(),
+                                          pfes.GetRestrictionMatrix(),
+                                          local_)
+{
+   MFEM_ASSERT(pfes.Conforming(), "internal error");
+   MFEM_ASSERT(pfes.GetRestrictionMatrix()->Height() == pfes.GetTrueVSize(), "");
+}
+
 static void ExtractSubVector(const int N,
                              const Array<int> &indices,
                              const Vector &in, Vector &out)

fem/pfespace.hpp

@@ -414,9 +414,14 @@ protected:
    bool local;
 
 public:
+   ConformingProlongationOperator(int lsize, const GroupCommunicator &gc_,
+                                  bool local_=false);
+
    ConformingProlongationOperator(const ParFiniteElementSpace &pfes,
                                   bool local_=false);
 
+   const GroupCommunicator &GetGroupCommunicator() const;
+
    virtual void Mult(const Vector &x, Vector &y) const;
 
    virtual void MultTranspose(const Vector &x, Vector &y) const;
@@ -461,6 +466,10 @@ protected:
    void ReduceEndAssemble(Vector &dst) const;
 
 public:
+   DeviceConformingProlongationOperator(
+      const GroupCommunicator &gc_, const SparseMatrix *R,
+      bool local_=false);
+
    DeviceConformingProlongationOperator(const ParFiniteElementSpace &pfes,
                                         bool local_=false);
 

linalg/hypre.cpp

@@ -1340,10 +1340,14 @@ static void get_sorted_rows_cols(const Array<int> &rows_cols,
    rows_cols.HostRead();
    hypre_sorted.SetSize(rows_cols.Size());
    bool sorted = true;
-   for (int i = 0; i < rows_cols.Size(); i++)
    {
-      hypre_sorted[i] = rows_cols[i];
-      if (i && rows_cols[i-1] > rows_cols[i]) { sorted = false; }
+      auto in_data = rows_cols.HostRead();
+      auto out_data = hypre_sorted.HostReadWrite();
+      for (int i = 0; i < rows_cols.Size(); i++)
+      {
+         out_data[i] = in_data[i];
+         if (i && in_data[i-1] > in_data[i]) { sorted = false; }
+      }
    }
    if (!sorted) { hypre_sorted.Sort(); }
 }
@@ -1432,7 +1436,7 @@ void HypreParMatrix::EliminateRowsCols(const Array<int> &rows_cols,
    get_sorted_rows_cols(rows_cols, rc_sorted);
 
    internal::hypre_ParCSRMatrixEliminateAXB(
-      A, rc_sorted.Size(), rc_sorted.GetData(), X, B);
+      A, rc_sorted.Size(), rc_sorted.HostReadWrite(), X, B);
 }
 
 HypreParMatrix* HypreParMatrix::EliminateRowsCols(const Array<int> &rows_cols)
@@ -1442,7 +1446,7 @@ HypreParMatrix* HypreParMatrix::EliminateRowsCols(const Array<int> &rows_cols)
 
    hypre_ParCSRMatrix* Ae;
    internal::hypre_ParCSRMatrixEliminateAAe(
-      A, &Ae, rc_sorted.Size(), rc_sorted.GetData());
+      A, &Ae, rc_sorted.Size(), rc_sorted.HostReadWrite());
 
    return new HypreParMatrix(Ae);
 }
@@ -1454,7 +1458,7 @@ HypreParMatrix* HypreParMatrix::EliminateCols(const Array<int> &cols)
 
    hypre_ParCSRMatrix* Ae;
    internal::hypre_ParCSRMatrixEliminateAAe(
-      A, &Ae, rc_sorted.Size(), rc_sorted.GetData(), 1);
+      A, &Ae, rc_sorted.Size(), rc_sorted.HostReadWrite(), 1);
 
    return new HypreParMatrix(Ae);
 }
@@ -1466,7 +1470,7 @@ void HypreParMatrix::EliminateRows(const Array<int> &rows)
       Array<HYPRE_Int> r_sorted;
       get_sorted_rows_cols(rows, r_sorted);
       internal::hypre_ParCSRMatrixEliminateRows(A, r_sorted.Size(),
-                                                r_sorted.GetData());
+                                                r_sorted.HostReadWrite());
    }
 }