Reducing the two nested GMRES solvers to just one. Also implemented a complex version of AMS.

PA for GradientInterpolator and IdentityInterpolator

CHANGELOG

@@ -11,6 +11,12 @@
 Version 4.2.1 (development)
 ===========================
 
+- Added high-order matrix-free auxiliary Maxwell solver for H(curl) problems,
+  as described in Barker and Kolev 2020 (https://doi.org/10.1002/nla.2348).
+
+- Added matrix-free GPU-enabled implementations of GradientInterpolator and
+  IdentityInterpolator.
+
 - Added interface to MUMPS direct solver. Its usage is demonstrated in ex25p.
   See http://mumps.enseeiht.fr/ for more details. Supported versions >= 5.1.1.
 

examples/ex1p_complex.cpp

@@ -0,0 +1,435 @@
+//                       MFEM Example 3 - Parallel Version
+//
+// Compile with: make ex3p_complex
+//
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+// Exact solution, E, and r.h.s., f. See below for implementation.
+double E_exact(const Vector &);
+void gradE_exact(const Vector &, Vector &);
+double f_exact(const Vector &);
+double freq = 1.0, kappa;
+int dim;
+
+
+#define COMPLEX_VERSION
+#define NEUMANN
+
+const double omega = 1.4;
+const double eps = 1.0e-8;
+
+
+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);
+
+   // 2. Parse command-line options.
+   const char *mesh_file = "../data/beam-tet.mesh";
+   int order = 1;
+   bool static_cond = false;
+   bool pa = false;
+   const char *device_config = "cpu";
+   bool visualization = 1;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree).");
+   args.AddOption(&freq, "-f", "--frequency", "Set the frequency for the exact"
+                  " solution.");
+   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
+                  "--no-static-condensation", "Enable static condensation.");
+   args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
+                  "--no-partial-assembly", "Enable Partial Assembly.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      MPI_Finalize();
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+   kappa = freq * M_PI;
+
+   // 3. Enable hardware devices such as GPUs, and programming models such as
+   //    CUDA, OCCA, RAJA and OpenMP based on command line options.
+   Device device(device_config);
+   if (myid == 0) { device.Print(); }
+
+   // 4. Read the (serial) mesh from the given mesh file on all processors.  We
+   //    can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
+   //    and volume meshes with the same code.
+   Mesh *mesh = new Mesh(mesh_file, 1, 1);
+   dim = mesh->Dimension();
+   int sdim = mesh->SpaceDimension();
+
+   // 5. Refine the serial mesh on all processors to increase the resolution. In
+   //    this example we do 'ref_levels' of uniform refinement. We choose
+   //    'ref_levels' to be the largest number that gives a final mesh with no
+   //    more than 1,000 elements.
+   {
+      int ref_levels = (int)floor(log(1000./mesh->GetNE())/log(2.)/dim);
+      for (int l = 0; l < ref_levels; l++)
+      {
+         mesh->UniformRefinement();
+      }
+   }
+
+   // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
+   //    this mesh further in parallel to increase the resolution. Once the
+   //    parallel mesh is defined, the serial mesh can be deleted. Tetrahedral
+   //    meshes need to be reoriented before we can define high-order Nedelec
+   //    spaces on them.
+   ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
+   delete mesh;
+   {
+      int par_ref_levels = 0;
+      for (int l = 0; l < par_ref_levels; l++)
+      {
+         pmesh->UniformRefinement();
+      }
+   }
+   pmesh->ReorientTetMesh();
+
+   // 7. Define a parallel finite element space on the parallel mesh. Here we
+   //    use the Nedelec finite elements of the specified order.
+   FiniteElementCollection *fec = new H1_FECollection(order, dim);
+   ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec);
+   HYPRE_Int size = fespace->GlobalTrueVSize();
+   if (myid == 0)
+   {
+      cout << "Number of finite element unknowns: " << size << endl;
+   }
+
+   // 8. Determine the list of true (i.e. parallel conforming) essential
+   //    boundary dofs. In this example, the boundary conditions are defined
+   //    by marking all the boundary attributes from the mesh as essential
+   //    (Dirichlet) and converting them to a list of true dofs.
+   Array<int> ess_tdof_list;
+
+#ifndef NEUMANN
+   if (pmesh->bdr_attributes.Size())
+   {
+      Array<int> ess_bdr(pmesh->bdr_attributes.Max());
+      ess_bdr = 1;
+      fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+#endif
+
+   //const double imscale = 0.0;
+   const double imscale = -omega;
+
+   Coefficient *im = new ConstantCoefficient(imscale);  // im part
+   //Coefficient *im = new ConstantCoefficient(0.0);  // im part
+
+   FunctionCoefficient E_coef(E_exact);
+   VectorFunctionCoefficient grad_E(sdim, gradE_exact);
+   ProductCoefficient omegaE(imscale, E_coef);
+
+   // 9. Set up the parallel linear form b(.) which corresponds to the
+   //    right-hand side of the FEM linear system, which in this case is
+   //    (f,phi_i) where f is given by the function f_exact and phi_i are the
+   //    basis functions in the finite element fespace.
+   FunctionCoefficient f(f_exact);
+#ifdef COMPLEX_VERSION
+   ParComplexLinearForm *b = new ParComplexLinearForm(fespace);
+   b->AddDomainIntegrator(new DomainLFIntegrator(f), NULL);
+#ifdef NEUMANN
+   b->AddBoundaryIntegrator(NULL, new BoundaryNormalLFIntegrator(grad_E));
+#endif
+   b->AddBoundaryIntegrator(NULL, new BoundaryLFIntegrator(omegaE)); // im part
+#endif
+
+   b->Assemble();
+
+   // 10. Define the solution vector x as a parallel finite element grid function
+   //     corresponding to fespace. Initialize x by projecting the exact
+   //     solution. Note that only values from the boundary edges will be used
+   //     when eliminating the non-homogeneous boundary condition to modify the
+   //     r.h.s. vector b.
+   /*
+   ParGridFunction x(fespace);
+   VectorFunctionCoefficient E(sdim, E_exact);
+   x.ProjectCoefficient(E);
+   */
+
+#ifdef COMPLEX_VERSION
+   // Complex version
+   ParComplexGridFunction x(fespace);
+   x = 0.0;
+   ConstantCoefficient E_im(0.0);
+   //x.ProjectBdrCoefficientTangent(E_Re, E_Im, ess_bdr);
+   x.ProjectCoefficient(E_coef, E_im);
+#endif
+
+   // 11. Set up the parallel bilinear form corresponding to the EM diffusion
+   //     operator curl muinv curl + sigma I, by adding the curl-curl and the
+   //     mass domain integrators.
+   Coefficient *muinv = new ConstantCoefficient(1.0);
+   Coefficient *epscoef = new ConstantCoefficient(eps);
+   Coefficient *imabs = new ConstantCoefficient(fabs(imscale));  // im part
+
+#ifdef COMPLEX_VERSION
+   // Complex version
+   ParSesquilinearForm *a = new ParSesquilinearForm(fespace);
+   if (pa) { a->SetAssemblyLevel(AssemblyLevel::PARTIAL); }
+   a->AddDomainIntegrator(new DiffusionIntegrator(*muinv), NULL);
+   a->AddDomainIntegrator(new MassIntegrator(*epscoef), NULL);
+   a->AddBoundaryIntegrator(NULL, new MassIntegrator(*im));  // im part
+#endif
+
+   // 12. Assemble the parallel bilinear form and the corresponding linear
+   //     system, applying any necessary transformations such as: parallel
+   //     assembly, eliminating boundary conditions, applying conforming
+   //     constraints for non-conforming AMR, static condensation, etc.
+   //if (static_cond) { a->EnableStaticCondensation(); }
+   a->Assemble();
+
+   OperatorPtr A;
+   Vector B, X;
+   a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
+
+   ParBilinearForm a_Re(fespace);
+   a_Re.AddDomainIntegrator(new DiffusionIntegrator(*muinv));
+   a_Re.AddDomainIntegrator(new MassIntegrator(*epscoef));
+
+   if (pa) { a_Re.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
+   a_Re.Assemble();
+
+   OperatorPtr A_Re;
+   a_Re.FormSystemMatrix(ess_tdof_list, A_Re);
+
+   ParBilinearForm a_Im(fespace);
+   a_Im.AddBoundaryIntegrator(new MassIntegrator(*imabs));
+   a_Im.Assemble();
+
+   OperatorPtr A_Im;
+   a_Im.FormSystemMatrix(ess_tdof_list, A_Im);
+
+   // 13. Solve the system AX=B using PCG with the AMS preconditioner from hypre
+   //     (in the full assembly case) or CG with Jacobi preconditioner (in the
+   //     partial assembly case).
+
+   Array<int> offsets(3);
+   offsets[0] = 0;
+   offsets[1] = fespace->GetTrueVSize();
+   offsets[2] = fespace->GetTrueVSize();
+   offsets.PartialSum();
+
+   //OperatorJacobiSmoother massJacobi(a_Im, ess_tdof_list);
+
+   StopWatch sw;
+   sw.Clear();
+   sw.Start();
+
+   if (pa) // Jacobi preconditioning in partial assembly mode
+   {
+      MFEM_VERIFY(false, "TODO");
+      //OperatorJacobiSmoother Jacobi(*a, ess_tdof_list);
+
+      CGSolver cg(MPI_COMM_WORLD);
+      cg.SetRelTol(1e-12);
+      cg.SetMaxIter(1000);
+      cg.SetPrintLevel(1);
+      cg.SetOperator(*A);
+      //cg.SetPreconditioner(Jacobi);
+      cg.Mult(B, X);
+   }
+   else
+   {
+      if (myid == 0)
+      {
+         cout << "Size of linear system: "
+              << A.As<HypreParMatrix>()->GetGlobalNumRows() << endl;
+      }
+
+      HypreBoomerAMG amg(*A_Re.As<HypreParMatrix>());
+
+#ifdef COMPLEX_VERSION
+      BlockDiagonalPreconditioner BlockDP(offsets);
+      BlockDP.SetDiagonalBlock(0, &amg);
+      BlockDP.SetDiagonalBlock(1, &amg);
+
+      Complex_PMHSS PMHSS(A_Re.Ptr(), A_Im.Ptr(), &BlockDP, NULL, 1.0);
+
+      ComplexOperator AspdComplex(A_Re.Ptr(), A_Im.Ptr(), false, false);
+
+      GMRESSolver PMHSSgmres(MPI_COMM_WORLD);
+      PMHSSgmres.SetPrintLevel(1);
+      PMHSSgmres.SetKDim(100);
+      PMHSSgmres.SetMaxIter(100);
+      PMHSSgmres.SetRelTol(1e-6);
+      PMHSSgmres.SetAbsTol(0.0);
+      PMHSSgmres.SetOperator(AspdComplex);
+      PMHSSgmres.SetPreconditioner(PMHSS);
+
+      GMRESSolver gmres(MPI_COMM_WORLD);
+      gmres.SetPrintLevel(1);
+      gmres.SetKDim(1000);
+      gmres.SetMaxIter(100);
+      gmres.SetRelTol(1e-8);
+      gmres.SetAbsTol(0.0);
+      gmres.SetOperator(*A);
+      //gmres.SetPreconditioner(BlockDP);
+      gmres.SetPreconditioner(PMHSS);
+      //gmres.SetPreconditioner(PMHSSgmres);
+#else
+      GMRESSolver gmres(MPI_COMM_WORLD);
+      gmres.SetPrintLevel(1);
+      gmres.SetKDim(1000);
+      gmres.SetMaxIter(100);
+      gmres.SetRelTol(1e-8);
+      gmres.SetAbsTol(0.0);
+      gmres.SetOperator(*A);
+      gmres.SetPreconditioner(ams);
+#endif
+
+      gmres.Mult(B, X);
+   }
+
+   sw.Stop();
+   mfem::out << "Total solve time " <<sw.RealTime() << endl;
+
+   // 14. Recover the parallel grid function corresponding to X. This is the
+   //     local finite element solution on each processor.
+   a->RecoverFEMSolution(X, *b, x);
+
+   // 15. Compute and print the L^2 norm of the error.
+   {
+#ifdef COMPLEX_VERSION
+      double err = x.real().ComputeL2Error(E_coef);
+#else
+      double err = x.ComputeL2Error(E_coef);
+#endif
+      if (myid == 0)
+      {
+         cout << "\n|| E_h - E ||_{L^2} = " << err << '\n' << endl;
+      }
+   }
+
+   // 16. Save the refined mesh and the solution in parallel. This output can
+   //     be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
+   {
+      ostringstream mesh_name, sol_name;
+      mesh_name << "mesh." << setfill('0') << setw(6) << myid;
+      sol_name << "sol." << setfill('0') << setw(6) << myid;
+
+      ofstream mesh_ofs(mesh_name.str().c_str());
+      mesh_ofs.precision(8);
+      pmesh->Print(mesh_ofs);
+
+      ofstream sol_ofs(sol_name.str().c_str());
+      sol_ofs.precision(8);
+#ifdef COMPLEX_VERSION
+      x.real().Save(sol_ofs);
+#else
+      x.Save(sol_ofs);
+#endif
+   }
+
+   // 17. Send the solution by socket to a GLVis server.
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream sol_sock(vishost, visport);
+      sol_sock << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock.precision(8);
+#ifdef COMPLEX_VERSION
+      sol_sock << "solution\n" << *pmesh << x.real() << flush;
+#else
+      sol_sock << "solution\n" << *pmesh << x << flush;
+#endif
+   }
+
+   // 18. Free the used memory.
+   delete a;
+   delete muinv;
+   delete b;
+   delete fespace;
+   delete fec;
+   delete pmesh;
+
+   MPI_Finalize();
+
+   return 0;
+}
+
+#define VERSION_COS
+
+double E_exact(const Vector &x)
+{
+   if (dim == 3)
+   {
+#ifdef VERSION_COS
+      return cos(kappa * x(0)) * cos(kappa * x(1)) * cos(kappa * x(2));
+#else
+      return sin(kappa * x(0)) * sin(kappa * x(1)) * sin(kappa * x(2));
+#endif
+   }
+   else
+   {
+      return 0.0;
+   }
+}
+
+void gradE_exact(const Vector &x, Vector &grad)
+{
+   if (dim == 3)
+   {
+#ifdef VERSION_COS
+      grad(0) = -kappa * sin(kappa * x(0)) * cos(kappa * x(1)) * cos(kappa * x(2));
+      grad(1) = -kappa * sin(kappa * x(1)) * cos(kappa * x(0)) * cos(kappa * x(2));
+      grad(2) = -kappa * sin(kappa * x(2)) * cos(kappa * x(0)) * cos(kappa * x(1));
+#else
+      grad(0) = kappa * cos(kappa * x(0)) * sin(kappa * x(1)) * sin(kappa * x(2));
+      grad(1) = kappa * cos(kappa * x(1)) * sin(kappa * x(0)) * sin(kappa * x(2));
+      grad(2) = kappa * cos(kappa * x(2)) * sin(kappa * x(0)) * sin(kappa * x(1));
+#endif
+   }
+   else
+   {
+      MFEM_VERIFY(false, "");
+   }
+}
+
+// (grad u, grad v) + eps (u, v) = <grad u . n, v> - (div grad u, v) + eps (u, v)
+double f_exact(const Vector &x)
+{
+   if (dim == 3)
+   {
+      const double c = 3.0 * kappa * kappa;
+#ifdef VERSION_COS
+      return (eps + c) * cos(kappa * x(0)) * cos(kappa * x(1)) * cos(kappa * x(2));
+#else
+      return (eps + c) * sin(kappa * x(0)) * sin(kappa * x(1)) * sin(kappa * x(2));
+#endif
+   }
+   else
+   {
+      return 0.0;
+   }
+}

examples/ex3p.cpp

@@ -69,7 +69,10 @@ int main(int argc, char *argv[])
    bool static_cond = false;
    bool pa = false;
    const char *device_config = "cpu";
-   bool visualization = 1;
+   bool visualization = true;
+#ifdef MFEM_USE_AMGX
+   bool useAmgX = false;
+#endif
 
    OptionsParser args(argc, argv);
    args.AddOption(&mesh_file, "-m", "--mesh",
@@ -87,6 +90,11 @@ int main(int argc, char *argv[])
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                   "--no-visualization",
                   "Enable or disable GLVis visualization.");
+#ifdef MFEM_USE_AMGX
+   args.AddOption(&useAmgX, "-amgx", "--useAmgX", "-no-amgx",
+                  "--no-useAmgX",
+                  "Enable or disable AmgX in MatrixFreeAMS.");
+#endif
 
    args.Parse();
    if (!args.Good())
@@ -159,9 +167,10 @@ int main(int argc, char *argv[])
    //    by marking all the boundary attributes from the mesh as essential
    //    (Dirichlet) and converting them to a list of true dofs.
    Array<int> ess_tdof_list;
+   Array<int> ess_bdr;
    if (pmesh->bdr_attributes.Size())
    {
-      Array<int> ess_bdr(pmesh->bdr_attributes.Max());
+      ess_bdr.SetSize(pmesh->bdr_attributes.Max());
       ess_bdr = 1;
       fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
    }
@@ -205,20 +214,20 @@ int main(int argc, char *argv[])
    Vector B, X;
    a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
 
-   // 13. Solve the system AX=B using PCG with the AMS preconditioner from hypre
-   //     (in the full assembly case) or CG with Jacobi preconditioner (in the
-   //     partial assembly case).
-
-   if (pa) // Jacobi preconditioning in partial assembly mode
+   // 13. Solve the system AX=B using PCG with an AMS preconditioner.
+   if (pa)
    {
-      OperatorJacobiSmoother Jacobi(*a, ess_tdof_list);
-
+#ifdef MFEM_USE_AMGX
+      MatrixFreeAMS ams(*a, *A, *fespace, muinv, sigma, NULL, ess_bdr, useAmgX);
+#else
+      MatrixFreeAMS ams(*a, *A, *fespace, muinv, sigma, NULL, ess_bdr);
+#endif
       CGSolver cg(MPI_COMM_WORLD);
       cg.SetRelTol(1e-12);
       cg.SetMaxIter(1000);
       cg.SetPrintLevel(1);
       cg.SetOperator(*A);
-      cg.SetPreconditioner(Jacobi);
+      cg.SetPreconditioner(ams);
       cg.Mult(B, X);
    }
    else

examples/ex3p_complex.cpp

@@ -0,0 +1,517 @@
+//                       MFEM Example 3 - Parallel Version
+//
+// Compile with: make ex3p_complex
+//
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+// Exact solution, E, and r.h.s., f. See below for implementation.
+void E_exact(const Vector &, Vector &);
+void curlE_exact(const Vector &, Vector &);
+void f_exact(const Vector &, Vector &);
+double freq = 1.0, kappa;
+int dim;
+
+
+#define COMPLEX_VERSION
+#define NEUMANN
+#define INDEFINITE
+
+const double omega = 1.4;
+
+
+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);
+
+   // 2. Parse command-line options.
+   const char *mesh_file = "../data/beam-tet.mesh";
+   int order = 1;
+   bool static_cond = false;
+   bool pa = false;
+   const char *device_config = "cpu";
+   bool visualization = 1;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree).");
+   args.AddOption(&freq, "-f", "--frequency", "Set the frequency for the exact"
+                  " solution.");
+   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
+                  "--no-static-condensation", "Enable static condensation.");
+   args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
+                  "--no-partial-assembly", "Enable Partial Assembly.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      MPI_Finalize();
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+   kappa = freq * M_PI;
+
+   // 3. Enable hardware devices such as GPUs, and programming models such as
+   //    CUDA, OCCA, RAJA and OpenMP based on command line options.
+   Device device(device_config);
+   if (myid == 0) { device.Print(); }
+
+   // 4. Read the (serial) mesh from the given mesh file on all processors.  We
+   //    can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
+   //    and volume meshes with the same code.
+   Mesh *mesh = new Mesh(mesh_file, 1, 1);
+   dim = mesh->Dimension();
+   int sdim = mesh->SpaceDimension();
+
+   // 5. Refine the serial mesh on all processors to increase the resolution. In
+   //    this example we do 'ref_levels' of uniform refinement. We choose
+   //    'ref_levels' to be the largest number that gives a final mesh with no
+   //    more than 1,000 elements.
+   {
+      int ref_levels = (int)floor(log(1000./mesh->GetNE())/log(2.)/dim);
+      for (int l = 0; l < ref_levels; l++)
+      {
+         mesh->UniformRefinement();
+      }
+   }
+
+   // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
+   //    this mesh further in parallel to increase the resolution. Once the
+   //    parallel mesh is defined, the serial mesh can be deleted. Tetrahedral
+   //    meshes need to be reoriented before we can define high-order Nedelec
+   //    spaces on them.
+   ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
+   delete mesh;
+   {
+      int par_ref_levels = 0;
+      for (int l = 0; l < par_ref_levels; l++)
+      {
+         pmesh->UniformRefinement();
+      }
+   }
+   pmesh->ReorientTetMesh();
+
+   // 7. Define a parallel finite element space on the parallel mesh. Here we
+   //    use the Nedelec finite elements of the specified order.
+   FiniteElementCollection *fec = new ND_FECollection(order, dim);
+   ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec);
+   HYPRE_Int size = fespace->GlobalTrueVSize();
+   if (myid == 0)
+   {
+      cout << "Number of finite element unknowns: " << size << endl;
+   }
+
+   // 8. Determine the list of true (i.e. parallel conforming) essential
+   //    boundary dofs. In this example, the boundary conditions are defined
+   //    by marking all the boundary attributes from the mesh as essential
+   //    (Dirichlet) and converting them to a list of true dofs.
+   Array<int> ess_tdof_list;
+   Array<int> ess_bdr;
+   ess_bdr.SetSize(pmesh->bdr_attributes.Max());
+   ess_bdr = 0;
+
+#ifndef NEUMANN
+   if (pmesh->bdr_attributes.Size())
+   {
+      ess_bdr = 1;
+      fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+#endif
+
+   //const double imscale = 0.0;
+   const double imscale = omega;
+
+   Coefficient *im = new ConstantCoefficient(imscale);  // im part
+   //Coefficient *im = new ConstantCoefficient(0.0);  // im part
+
+   VectorFunctionCoefficient E_Re(sdim, E_exact);
+   VectorFunctionCoefficient curlE_Re(sdim, curlE_exact);
+
+   ScalarVectorProductCoefficient omegaE(imscale, E_Re);  // im part
+   //ScalarVectorProductCoefficient omegaE(0.0, E_Re);  // im part
+
+   // 9. Set up the parallel linear form b(.) which corresponds to the
+   //    right-hand side of the FEM linear system, which in this case is
+   //    (f,phi_i) where f is given by the function f_exact and phi_i are the
+   //    basis functions in the finite element fespace.
+   VectorFunctionCoefficient f(sdim, f_exact);
+#ifdef COMPLEX_VERSION
+   ParComplexLinearForm *b = new ParComplexLinearForm(fespace);
+   b->AddDomainIntegrator(new VectorFEDomainLFIntegrator(f), NULL);
+   b->AddBoundaryIntegrator(NULL,
+                            new VectorFEDomainLFIntegrator(omegaE));  // im part
+#else
+   // Real version
+   ParLinearForm *b = new ParLinearForm(fespace);
+   b->AddDomainIntegrator(new VectorFEDomainLFIntegrator(f));
+#endif
+
+#ifdef NEUMANN
+   b->AddBoundaryIntegrator(new VectorFEBoundaryTangentLFIntegrator(curlE_Re),
+                            NULL);
+#endif
+
+   b->Assemble();
+
+   // 10. Define the solution vector x as a parallel finite element grid function
+   //     corresponding to fespace. Initialize x by projecting the exact
+   //     solution. Note that only values from the boundary edges will be used
+   //     when eliminating the non-homogeneous boundary condition to modify the
+   //     r.h.s. vector b.
+   /*
+   ParGridFunction x(fespace);
+   VectorFunctionCoefficient E(sdim, E_exact);
+   x.ProjectCoefficient(E);
+   */
+
+#ifdef COMPLEX_VERSION
+   // Complex version
+   ParComplexGridFunction x(fespace);
+   x = 0.0;
+   Vector zero(sdim);
+   zero = 0.0;
+   VectorConstantCoefficient E_Im(zero);
+   //x.ProjectBdrCoefficientTangent(E_Re, E_Im, ess_bdr);
+   x.ProjectCoefficient(E_Re, E_Im);
+#else
+   ParGridFunction x(fespace);
+   x = 0.0;
+   x.ProjectCoefficient(E_Re);
+#endif
+
+   // 11. Set up the parallel bilinear form corresponding to the EM diffusion
+   //     operator curl muinv curl + sigma I, by adding the curl-curl and the
+   //     mass domain integrators.
+   Coefficient *muinv = new ConstantCoefficient(1.0);
+#ifdef INDEFINITE
+   Coefficient *sigma = new ConstantCoefficient(
+      -omega*omega); // indefinite -, definite +
+#else
+   Coefficient *sigma = new ConstantCoefficient(
+      omega*omega); // indefinite -, definite +
+#endif
+   Coefficient *abssigma = new ConstantCoefficient(omega*omega);
+   Coefficient *imabs = new ConstantCoefficient(imscale);  // im part
+   //Coefficient *imabs = new ConstantCoefficient(0.0);  // im part
+   //Coefficient *im = new ConstantCoefficient(0.0);
+
+#ifdef COMPLEX_VERSION
+   // Complex version
+   ParSesquilinearForm *a = new ParSesquilinearForm(fespace);
+   if (pa) { a->SetAssemblyLevel(AssemblyLevel::PARTIAL); }
+   a->AddDomainIntegrator(new CurlCurlIntegrator(*muinv), NULL);
+   //a->AddDomainIntegrator(new VectorFEMassIntegrator(*sigma), new VectorFEMassIntegrator(*im));
+   a->AddDomainIntegrator(new VectorFEMassIntegrator(*sigma), NULL);
+   a->AddBoundaryIntegrator(NULL, new VectorFEMassIntegrator(*im));  // im part
+#else
+   // Real version
+   ParBilinearForm *a = new ParBilinearForm(fespace);
+   if (pa) { a->SetAssemblyLevel(AssemblyLevel::PARTIAL); }
+   a->AddDomainIntegrator(new CurlCurlIntegrator(*muinv));
+   //a->AddDomainIntegrator(new VectorFEMassIntegrator(*sigma), new VectorFEMassIntegrator(*im));
+   a->AddDomainIntegrator(new VectorFEMassIntegrator(*sigma));
+#endif
+
+   // 12. Assemble the parallel bilinear form and the corresponding linear
+   //     system, applying any necessary transformations such as: parallel
+   //     assembly, eliminating boundary conditions, applying conforming
+   //     constraints for non-conforming AMR, static condensation, etc.
+   //if (static_cond) { a->EnableStaticCondensation(); }
+   a->Assemble();
+
+   OperatorPtr A;
+   Vector B, X;
+   a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B);
+
+   ParBilinearForm a_Re(fespace);
+   a_Re.AddDomainIntegrator(new CurlCurlIntegrator(*muinv));
+   a_Re.AddDomainIntegrator(new VectorFEMassIntegrator(*abssigma));
+
+   //if (pa) { a_Re.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
+   a_Re.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+   a_Re.Assemble();
+
+   OperatorPtr A_Re;
+   a_Re.FormSystemMatrix(ess_tdof_list, A_Re);
+
+   ParBilinearForm a_Im(fespace);
+   a_Im.AddBoundaryIntegrator(new VectorFEMassIntegrator(*imabs));
+   a_Im.Assemble();
+
+   OperatorPtr A_Im;
+   a_Im.FormSystemMatrix(ess_tdof_list, A_Im);
+
+   // 13. Solve the system AX=B using PCG with the AMS preconditioner from hypre
+   //     (in the full assembly case) or CG with Jacobi preconditioner (in the
+   //     partial assembly case).
+
+   Array<int> offsets(3);
+   offsets[0] = 0;
+   offsets[1] = fespace->GetTrueVSize();
+   offsets[2] = fespace->GetTrueVSize();
+   offsets.PartialSum();
+
+   //OperatorJacobiSmoother massJacobi(a_Im, ess_tdof_list);
+
+   StopWatch sw;
+   sw.Clear();
+   sw.Start();
+
+   if (pa) // Jacobi preconditioning in partial assembly mode
+   {
+      MFEM_VERIFY(false, "TODO");
+      //OperatorJacobiSmoother Jacobi(*a, ess_tdof_list);
+
+      CGSolver cg(MPI_COMM_WORLD);
+      cg.SetRelTol(1e-12);
+      cg.SetMaxIter(1000);
+      cg.SetPrintLevel(1);
+      cg.SetOperator(*A);
+      //cg.SetPreconditioner(Jacobi);
+      cg.Mult(B, X);
+   }
+   else
+   {
+      if (myid == 0)
+      {
+         cout << "Size of linear system: "
+              << A.As<HypreParMatrix>()->GetGlobalNumRows() << endl;
+      }
+
+      //HypreAMS ams(*A_Re.As<HypreParMatrix>(), fespace);
+
+      // One option is to use the standard real-valued MatrixFreeAMS to precondition
+      // the real part of the complex system in the PMHSS preconditioner (BlockDiagonalPreconditioner).
+      // Another option is to use complex MatrixFreeAMS to precondition the
+      // complex system without PMHSS and without a BlockDiagonalPreconditioner.
+      //#define COMPLEX_AMS
+
+#ifdef MFEM_USE_AMGX
+      bool useAmgX = false;
+      cout << "Built with AMGX, using AMGX " << useAmgX << endl;
+      MatrixFreeAMS ams(a_Re, *A_Re, *fespace, muinv, abssigma, im, imabs, NULL,
+                        ess_bdr, useAmgX);
+      MatrixFreeAMS ams(a_Re, *A_Re, *fespace, muinv, abssigma, NULL, NULL, ess_bdr,
+                        useAmgX);
+#ifdef COMPLEX_AMS
+      MFEM_VERIFY(false, "TODO");
+#endif
+
+#else
+      cout << "Not built with AMGX" << endl;
+#ifdef COMPLEX_AMS
+      MatrixFreeAMS ams(a_Re, *A_Re, A.Ptr(), *fespace, muinv, abssigma, im, imabs,
+                        NULL, ess_bdr);
+#else
+      MatrixFreeAMS ams(a_Re, *A_Re, NULL, *fespace, muinv, abssigma, NULL, NULL,
+                        NULL, ess_bdr);
+#endif
+#endif
+
+#ifdef COMPLEX_VERSION
+
+#ifdef COMPLEX_AMS
+      //MFEM_VERIFY(false, "TODO");
+#else
+      BlockDiagonalPreconditioner BlockDP(offsets);
+      BlockDP.SetDiagonalBlock(0, &ams);
+      BlockDP.SetDiagonalBlock(1, &ams);
+
+      /*
+      BlockDiagonalPreconditioner BlockDP_Im(offsets);
+      BlockDP_Im.SetDiagonalBlock(0, &massJacobi); // TODO: this won't work if it has zeros on diagonal
+      BlockDP_Im.SetDiagonalBlock(1, &massJacobi);
+      */
+
+      //Complex_PMHSS PMHSS(A_Re, A_Im, &BlockDP, &BlockDP_Im);
+      //Complex_PMHSS PMHSS(A_Re, A_Im, &BlockDP, NULL, 2.0 * omega);
+      //Complex_PMHSS PMHSS(A_Re, A_Im, &BlockDP, NULL, omega);
+      Complex_PMHSS PMHSS(A_Re.Ptr(), A_Im.Ptr(), &BlockDP, NULL, 1.0);
+
+      ComplexOperator AspdComplex(A_Re.Ptr(), A_Im.Ptr(), false, false);
+
+      GMRESSolver PMHSSgmres(MPI_COMM_WORLD);
+      PMHSSgmres.SetPrintLevel(1);
+      PMHSSgmres.SetKDim(100);
+      PMHSSgmres.SetMaxIter(100);
+      PMHSSgmres.SetRelTol(1e-6);
+      PMHSSgmres.SetAbsTol(0.0);
+      PMHSSgmres.SetOperator(AspdComplex);
+      PMHSSgmres.SetPreconditioner(PMHSS);
+#endif
+
+      GMRESSolver gmres(MPI_COMM_WORLD);
+      gmres.SetPrintLevel(1);
+      gmres.SetKDim(1000);
+      gmres.SetMaxIter(100);
+      gmres.SetRelTol(1e-8);
+      gmres.SetAbsTol(0.0);
+      gmres.SetOperator(*A);
+      //gmres.SetPreconditioner(BlockDP);
+#ifdef COMPLEX_AMS
+      //MFEM_VERIFY(false, "TODO");
+      gmres.SetPreconditioner(ams);
+#else
+      gmres.SetPreconditioner(PMHSS);
+      //gmres.SetPreconditioner(PMHSSgmres);
+#endif
+
+#else
+      GMRESSolver gmres(MPI_COMM_WORLD);
+      gmres.SetPrintLevel(1);
+      gmres.SetKDim(1000);
+      gmres.SetMaxIter(100);
+      gmres.SetRelTol(1e-8);
+      gmres.SetAbsTol(0.0);
+      gmres.SetOperator(*A);
+      gmres.SetPreconditioner(ams);
+#endif
+
+      gmres.Mult(B, X);
+   }
+
+   sw.Stop();
+   mfem::out << "Total solve time " <<sw.RealTime() << endl;
+
+   // 14. Recover the parallel grid function corresponding to X. This is the
+   //     local finite element solution on each processor.
+   a->RecoverFEMSolution(X, *b, x);
+
+   // 15. Compute and print the L^2 norm of the error.
+   {
+#ifdef COMPLEX_VERSION
+      double err = x.real().ComputeL2Error(E_Re);
+#else
+      double err = x.ComputeL2Error(E_Re);
+#endif
+      if (myid == 0)
+      {
+         cout << "\n|| E_h - E ||_{L^2} = " << err << '\n' << endl;
+      }
+   }
+
+   // 16. Save the refined mesh and the solution in parallel. This output can
+   //     be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
+   {
+      ostringstream mesh_name, sol_name;
+      mesh_name << "mesh." << setfill('0') << setw(6) << myid;
+      sol_name << "sol." << setfill('0') << setw(6) << myid;
+
+      ofstream mesh_ofs(mesh_name.str().c_str());
+      mesh_ofs.precision(8);
+      pmesh->Print(mesh_ofs);
+
+      ofstream sol_ofs(sol_name.str().c_str());
+      sol_ofs.precision(8);
+#ifdef COMPLEX_VERSION
+      x.real().Save(sol_ofs);
+#else
+      x.Save(sol_ofs);
+#endif
+   }
+
+   // 17. Send the solution by socket to a GLVis server.
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream sol_sock(vishost, visport);
+      sol_sock << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock.precision(8);
+#ifdef COMPLEX_VERSION
+      sol_sock << "solution\n" << *pmesh << x.real() << flush;
+#else
+      sol_sock << "solution\n" << *pmesh << x << flush;
+#endif
+   }
+
+   // 18. Free the used memory.
+   delete a;
+   delete sigma;
+   delete muinv;
+   delete b;
+   delete fespace;
+   delete fec;
+   delete pmesh;
+
+   MPI_Finalize();
+
+   return 0;
+}
+
+
+void E_exact(const Vector &x, Vector &E)
+{
+   if (dim == 3)
+   {
+      E(0) = sin(kappa * x(1));
+      E(1) = sin(kappa * x(2));
+      E(2) = sin(kappa * x(0));
+   }
+   else
+   {
+      E(0) = sin(kappa * x(1));
+      E(1) = sin(kappa * x(0));
+      if (x.Size() == 3) { E(2) = 0.0; }
+   }
+}
+
+void curlE_exact(const Vector &x, Vector &curl)
+{
+   if (dim == 3)
+   {
+      curl(0) = kappa * cos(kappa * x(2));
+      curl(1) = kappa * cos(kappa * x(0));
+      curl(2) = kappa * cos(kappa * x(1));
+   }
+   else
+   {
+      MFEM_VERIFY(false, "");
+   }
+}
+
+void f_exact(const Vector &x, Vector &f)
+{
+   if (dim == 3)
+   {
+      // indefinite -m, definite +m
+      const double c = kappa * kappa;
+#ifdef INDEFINITE
+      const double m = -omega * omega;
+#else
+      const double m = omega * omega;
+#endif
+      f(0) = (c + m) * sin(kappa * x(1));
+      f(1) = (c + m) * sin(kappa * x(2));
+      f(2) = (c + m) * sin(kappa * x(0));
+   }
+   else
+   {
+      f(0) = (1. + kappa * kappa) * sin(kappa * x(1));
+      f(1) = (1. + kappa * kappa) * sin(kappa * x(0));
+      if (x.Size() == 3) { f(2) = 0.0; }
+   }
+}

fem/bilinearform.cpp

@@ -1770,9 +1770,41 @@ MixedBilinearForm::~MixedBilinearForm()
    delete ext;
 }
 
+void DiscreteLinearOperator::SetAssemblyLevel(AssemblyLevel assembly_level)
+{
+   if (ext)
+   {
+      MFEM_ABORT("the assembly level has already been set!");
+   }
+   assembly = assembly_level;
+   switch (assembly)
+   {
+      case AssemblyLevel::LEGACYFULL:
+      case AssemblyLevel::FULL:
+         // Use the original implementation for now
+         break;
+      case AssemblyLevel::ELEMENT:
+         mfem_error("Element assembly not supported yet... stay tuned!");
+         break;
+      case AssemblyLevel::PARTIAL:
+         ext = new PADiscreteLinearOperatorExtension(this);
+         break;
+      case AssemblyLevel::NONE:
+         mfem_error("Matrix-free action not supported yet... stay tuned!");
+         break;
+      default:
+         mfem_error("Unknown assembly level");
+   }
+}
 
 void DiscreteLinearOperator::Assemble(int skip_zeros)
 {
+   if (ext)
+   {
+      ext->Assemble();
+      return;
+   }
+
    Array<int> dom_vdofs, ran_vdofs;
    ElementTransformation *T;
    const FiniteElement *dom_fe, *ran_fe;

fem/bilinearform.hpp

@@ -376,6 +376,13 @@ public:
    /// Get the output finite element space prolongation matrix
    virtual const Operator *GetOutputProlongation() const
    { return GetProlongation(); }
+   /** @brief Returns the output fe space restriction matrix, transposed
+
+       Logically, this is the transpose of GetOutputRestriction, but in
+       practice it is convenient to have it in transposed form for
+       construction of RAP operators in matrix-free methods. */
+   virtual const Operator *GetOutputRestrictionTranspose() const
+   { return GetOutputProlongation(); }
    /// Get the output finite element space restriction matrix
    virtual const Operator *GetOutputRestriction() const
    { return GetRestriction(); }
@@ -977,9 +984,18 @@ public:
    /// Access all interpolators added with AddDomainInterpolator().
    Array<BilinearFormIntegrator*> *GetDI() { return &dbfi; }
 
+   /// Set the desired assembly level. The default is AssemblyLevel::FULL.
+   /** This method must be called before assembly. */
+   void SetAssemblyLevel(AssemblyLevel assembly_level);
+
    /** @brief Construct the internal matrix representation of the discrete
        linear operator. */
    virtual void Assemble(int skip_zeros = 1);
+
+   /** @brief Get the output finite element space restriction matrix in
+       transposed form. */
+   virtual const Operator *GetOutputRestrictionTranspose() const
+   { return test_fes->GetRestrictionTransposeOperator(); }
 };
 
 }

fem/bilinearform_ext.cpp

@@ -1021,7 +1021,6 @@ void PAMixedBilinearFormExtension::Update()
       localTrial.UseDevice(true);
       localTrial.SetSize(elem_restrict_trial->Height(),
                          Device::GetMemoryType());
-
    }
    if (elem_restrict_test)
    {
@@ -1221,4 +1220,134 @@ void PAMixedBilinearFormExtension::AssembleDiagonal_ADAt(const Vector &D,
    }
 }
 
+PADiscreteLinearOperatorExtension::PADiscreteLinearOperatorExtension(
+   DiscreteLinearOperator *linop) :
+   PAMixedBilinearFormExtension(linop)
+{
+}
+
+const
+Operator *PADiscreteLinearOperatorExtension::GetOutputRestrictionTranspose()
+const
+{
+   return a->GetOutputRestrictionTranspose();
+}
+
+void PADiscreteLinearOperatorExtension::Assemble()
+{
+   Array<BilinearFormIntegrator*> &integrators = *a->GetDBFI();
+   const int integratorCount = integrators.Size();
+   for (int i = 0; i < integratorCount; ++i)
+   {
+      integrators[i]->AssemblePA(*trialFes, *testFes);
+   }
+
+   test_multiplicity.UseDevice(true);
+   test_multiplicity.SetSize(elem_restrict_test->Width()); // l-vector
+   Vector ones(elem_restrict_test->Height()); // e-vector
+   ones = 1.0;
+
+   const ElementRestriction* elem_restrict =
+      dynamic_cast<const ElementRestriction*>(elem_restrict_test);
+   if (elem_restrict)
+   {
+      elem_restrict->MultTransposeUnsigned(ones, test_multiplicity);
+   }
+   else
+   {
+      mfem_error("A real ElementRestriction is required in this setting!");
+   }
+
+   auto tm = test_multiplicity.ReadWrite();
+   MFEM_FORALL(i, test_multiplicity.Size(),
+   {
+      tm[i] = 1.0 / tm[i];
+   });
+}
+
+void PADiscreteLinearOperatorExtension::AddMult(
+   const Vector &x, Vector &y, const double c) const
+{
+   Array<BilinearFormIntegrator*> &integrators = *a->GetDBFI();
+   const int iSz = integrators.Size();
+
+   // * G operation
+   SetupMultInputs(elem_restrict_trial, x, localTrial,
+                   elem_restrict_test, y, localTest, c);
+
+   // * B^TDB operation
+   for (int i = 0; i < iSz; ++i)
+   {
+      integrators[i]->AddMultPA(localTrial, localTest);
+   }
+
+   // do a kind of "set" rather than "add" in the below
+   // operation as compared to the BilinearForm case
+   // * G^T operation (kind of...)
+   const ElementRestriction* elem_restrict =
+      dynamic_cast<const ElementRestriction*>(elem_restrict_test);
+   if (elem_restrict)
+   {
+      tempY.SetSize(y.Size());
+      elem_restrict->MultLeftInverse(localTest, tempY);
+      y += tempY;
+   }
+   else
+   {
+      mfem_error("In this setting you need a real ElementRestriction!");
+   }
+}
+
+void PADiscreteLinearOperatorExtension::AddMultTranspose(
+   const Vector &x, Vector &y, const double c) const
+{
+   Array<BilinearFormIntegrator*> &integrators = *a->GetDBFI();
+   const int iSz = integrators.Size();
+
+   // do a kind of "set" rather than "add" in the below
+   // operation as compared to the BilinearForm case
+   // * G operation (kinda)
+   Vector xscaled(x);
+   MFEM_VERIFY(x.Size() == test_multiplicity.Size(), "Input vector of wrong size");
+   auto xs = xscaled.ReadWrite();
+   auto tm = test_multiplicity.Read();
+   MFEM_FORALL(i, x.Size(),
+   {
+      xs[i] *= tm[i];
+   });
+   SetupMultInputs(elem_restrict_test, xscaled, localTest,
+                   elem_restrict_trial, y, localTrial, c);
+
+   // * B^TD^TB operation
+   for (int i = 0; i < iSz; ++i)
+   {
+      integrators[i]->AddMultTransposePA(localTest, localTrial);
+   }
+
+   // * G^T operation
+   if (elem_restrict_trial)
+   {
+      tempY.SetSize(y.Size());
+      elem_restrict_trial->MultTranspose(localTrial, tempY);
+      y += tempY;
+   }
+   else
+   {
+      mfem_error("Trial ElementRestriction not defined");
+   }
+}
+
+void PADiscreteLinearOperatorExtension::FormRectangularSystemOperator(
+   const Array<int>& ess1, const Array<int>& ess2, OperatorHandle &A)
+{
+   const Operator *Pi = this->GetProlongation();
+   const Operator *RoT = this->GetOutputRestrictionTranspose();
+   Operator *rap = SetupRAP(Pi, RoT);
+
+   RectangularConstrainedOperator *Arco
+      = new RectangularConstrainedOperator(rap, ess1, ess2, rap != this);
+
+   A.Reset(Arco);
+}
+
 } // namespace mfem

fem/bilinearform_ext.hpp

@@ -21,6 +21,7 @@ namespace mfem
 
 class BilinearForm;
 class MixedBilinearForm;
+class DiscreteLinearOperator;
 
 /// Class extending the BilinearForm class to support different AssemblyLevels.
 /**  FA - Full Assembly
@@ -212,7 +213,7 @@ protected:
    mutable Vector localTrial, localTest, tempY;
    const Operator *elem_restrict_trial; // Not owned
    const Operator *elem_restrict_test;  // Not owned
-private:
+
    /// Helper function to set up inputs/outputs for Mult or MultTranspose
    void SetupMultInputs(const Operator *elem_restrict_x,
                         const Vector &x, Vector &localX,
@@ -258,6 +259,35 @@ public:
    void Update();
 };
 
+
+/**
+   @brief Partial assembly extension for DiscreteLinearOperator
+
+   This acts very much like PAMixedBilinearFormExtension, but its
+   FormRectangularSystemOperator implementation emulates 'Set' rather than
+   'Add' in the assembly case.
+*/
+class PADiscreteLinearOperatorExtension : public PAMixedBilinearFormExtension
+{
+public:
+   PADiscreteLinearOperatorExtension(DiscreteLinearOperator *linop);
+
+   /// Partial assembly of all internal integrators
+   void Assemble();
+
+   void AddMult(const Vector &x, Vector &y, const double c) const;
+
+   void AddMultTranspose(const Vector &x, Vector &y, const double c=1.0) const;
+
+   void FormRectangularSystemOperator(const Array<int>&, const Array<int>&,
+                                      OperatorHandle& A);
+
+   const Operator * GetOutputRestrictionTranspose() const;
+
+private:
+   Vector test_multiplicity;
+};
+
 }
 
 #endif

fem/bilininteg.cpp

@@ -22,14 +22,14 @@ namespace mfem
 
 void BilinearFormIntegrator::AssemblePA(const FiniteElementSpace&)
 {
-   mfem_error ("BilinearFormIntegrator::AssemblePA(...)\n"
+   mfem_error ("BilinearFormIntegrator::AssemblePA(fes)\n"
                "   is not implemented for this class.");
 }
 
 void BilinearFormIntegrator::AssemblePA(const FiniteElementSpace&,
                                         const FiniteElementSpace&)
 {
-   mfem_error ("BilinearFormIntegrator::AssemblePA(...)\n"
+   mfem_error ("BilinearFormIntegrator::AssemblePA(fes, fes)\n"
                "   is not implemented for this class.");
 }
 
@@ -92,7 +92,7 @@ void BilinearFormIntegrator::AddMultPA(const Vector &, Vector &) const
 
 void BilinearFormIntegrator::AddMultTransposePA(const Vector &, Vector &) const
 {
-   mfem_error ("BilinearFormIntegrator::MultAssembledTranspose(...)\n"
+   mfem_error ("BilinearFormIntegrator::AddMultTransposePA(...)\n"
                "   is not implemented for this class.");
 }
 

fem/bilininteg.hpp

@@ -1844,8 +1844,10 @@ protected:
 };
 
 /** Class for integrating the bilinear form a(u,v) := (Q grad u, v) where Q is a
-    scalar coefficient, and v is a vector with components v_i in the same space
-    as u. */
+    scalar coefficient, and v is a vector with components v_i in the same (H1) space
+    as u.
+
+    See also MixedVectorGradientIntegrator when v is in H(curl). */
 class GradientIntegrator : public BilinearFormIntegrator
 {
 protected:
@@ -2989,11 +2991,36 @@ class DiscreteInterpolator : public BilinearFormIntegrator { };
 class GradientInterpolator : public DiscreteInterpolator
 {
 public:
+   GradientInterpolator() : dofquad_fe(NULL) { }
+   virtual ~GradientInterpolator() { delete dofquad_fe; }
+
    virtual void AssembleElementMatrix2(const FiniteElement &h1_fe,
                                        const FiniteElement &nd_fe,
                                        ElementTransformation &Trans,
                                        DenseMatrix &elmat)
    { nd_fe.ProjectGrad(h1_fe, Trans, elmat); }
+
+   using BilinearFormIntegrator::AssemblePA;
+
+   /** @brief Setup method for PA data.
+
+       @param[in] trial_fes   H1 Lagrange space
+       @param[in] test_fes    H(curl) Nedelec space
+    */
+   virtual void AssemblePA(const FiniteElementSpace &trial_fes,
+                           const FiniteElementSpace &test_fes);
+
+   virtual void AddMultPA(const Vector &x, Vector &y) const;
+   virtual void AddMultTransposePA(const Vector &x, Vector &y) const;
+
+private:
+   /// 1D finite element that generates and owns the 1D DofToQuad maps below
+   FiniteElement * dofquad_fe;
+
+   bool B_id; // is the B basis operator (maps_C_C) the identity?
+   const DofToQuad *maps_C_C; // one-d map with Lobatto rows, Lobatto columns
+   const DofToQuad *maps_O_C; // one-d map with Legendre rows, Lobatto columns
+   int dim, ne, o_dofs1D, c_dofs1D;
 };
 
 
@@ -3008,6 +3035,24 @@ public:
                                        ElementTransformation &Trans,
                                        DenseMatrix &elmat)
    { ran_fe.Project(dom_fe, Trans, elmat); }
+
+   using BilinearFormIntegrator::AssemblePA;
+
+   virtual void AssemblePA(const FiniteElementSpace &trial_fes,
+                           const FiniteElementSpace &test_fes);
+
+   virtual void AddMultPA(const Vector &x, Vector &y) const;
+   virtual void AddMultTransposePA(const Vector &x, Vector &y) const;
+
+private:
+   /// 1D finite element that generates and owns the 1D DofToQuad maps below
+   FiniteElement * dofquad_fe;
+
+   const DofToQuad *maps_C_C; // one-d map with Lobatto rows, Lobatto columns
+   const DofToQuad *maps_O_C; // one-d map with Legendre rows, Lobatto columns
+   int dim, ne, o_dofs1D, c_dofs1D;
+
+   Vector pa_data;
 };
 
 

fem/bilininteg_gradient.cpp

@@ -334,7 +334,7 @@ static void PAGradientApplyTranspose2D(const int NE,
                                        const int q1d = 0)
 {
    // TODO
-   MFEM_ASSERT(false, "GradientPAApplyTranspose 3D not implemented.");
+   MFEM_ASSERT(false, "PAGradientApplyTranspose2D not implemented.");
 }
 
 // PA Gradient Apply 3D kernel

fem/fespace.cpp

@@ -2608,9 +2608,9 @@ const Operator &GridTransfer::MakeTrueOperator(
    else // Parallel() == true
    {
 #ifdef MFEM_USE_MPI
-      const SparseMatrix *out_R = fes_out.GetRestrictionMatrix();
       if (oper_type == Operator::Hypre_ParCSR)
       {
+         const SparseMatrix *out_R = fes_out.GetRestrictionMatrix();
          const ParFiniteElementSpace *pfes_in =
             dynamic_cast<const ParFiniteElementSpace *>(&fes_in);
          const ParFiniteElementSpace *pfes_out =
@@ -2638,6 +2638,7 @@ const Operator &GridTransfer::MakeTrueOperator(
       }
       else if (oper_type == Operator::ANY_TYPE)
       {
+         const Operator *out_R = fes_out.GetRestrictionOperator();
          t_oper.Reset(new TripleProductOperator(
                          out_R, &oper, fes_in.GetProlongationMatrix(),
                          false, false, false));

fem/fespace.hpp

@@ -330,6 +330,18 @@ public:
    virtual const Operator *GetProlongationMatrix() const
    { return GetConformingProlongation(); }
 
+   /// Return an operator that performs the transpose of GetRestrictionOperator
+   /** The returned operator is owned by the FiniteElementSpace. In serial this
+       is the same as GetProlongationMatrix() */
+   virtual const Operator *GetRestrictionTransposeOperator() const
+   { return GetConformingProlongation(); }
+
+   /// An abstract operator that performs the same action as GetRestrictionMatrix
+   /** In some cases this is an optimized matrix-free implementation. The
+       returned operator is owned by the FiniteElementSpace. */
+   virtual const Operator *GetRestrictionOperator() const
+   { return GetConformingRestriction(); }
+
    /// The returned SparseMatrix is owned by the FiniteElementSpace.
    virtual const SparseMatrix *GetRestrictionMatrix() const
    { return GetConformingRestriction(); }

fem/pbilinearform.cpp

@@ -565,6 +565,38 @@ HypreParMatrix* ParDiscreteLinearOperator::ParallelAssemble() const
    return RAP;
 }
 
+void ParDiscreteLinearOperator::ParallelAssemble(OperatorHandle &A)
+{
+   // construct the rectangular block-diagonal matrix dA
+   OperatorHandle dA(A.Type());
+   dA.MakeRectangularBlockDiag(domain_fes->GetComm(),
+                               range_fes->GlobalVSize(),
+                               domain_fes->GlobalVSize(),
+                               range_fes->GetDofOffsets(),
+                               domain_fes->GetDofOffsets(),
+                               mat);
+
+   OperatorHandle R_test_transpose(A.Type()), P_trial(A.Type());
+
+   // TODO - construct the Dof_TrueDof_Matrix directly in the required format.
+   R_test_transpose.ConvertFrom(range_fes->Dof_TrueDof_Matrix());
+   P_trial.ConvertFrom(domain_fes->Dof_TrueDof_Matrix());
+
+   A.MakeRAP(R_test_transpose, dA, P_trial);
+}
+
+void ParDiscreteLinearOperator::FormRectangularSystemMatrix(OperatorHandle &A)
+{
+   if (ext)
+   {
+      Array<int> empty;
+      ext->FormRectangularSystemOperator(empty, empty, A);
+      return;
+   }
+
+   mfem_error("not implemented!");
+}
+
 void ParDiscreteLinearOperator::GetParBlocks(Array2D<HypreParMatrix *> &blocks)
 const
 {

fem/pbilinearform.hpp

@@ -160,6 +160,9 @@ public:
    /// Get the parallel finite element space prolongation matrix
    virtual const Operator *GetProlongation() const
    { return pfes->GetProlongationMatrix(); }
+   /// Get the transpose of GetRestriction, useful for matrix-free RAP
+   virtual const Operator *GetRestrictionTranspose() const
+   { return pfes->GetRestrictionTransposeOperator(); }
    /// Get the parallel finite element space restriction matrix
    virtual const Operator *GetRestriction() const
    { return pfes->GetRestrictionMatrix(); }
@@ -301,10 +304,18 @@ public:
    /// Returns the matrix "assembled" on the true dofs
    HypreParMatrix *ParallelAssemble() const;
 
+   /** @brief Returns the matrix assembled on the true dofs, i.e.
+       @a A = R_test A_local P_trial, in the format (type id) specified by
+       @a A. */
+   void ParallelAssemble(OperatorHandle &A);
+
    /** Extract the parallel blocks corresponding to the vector dimensions of the
        domain and range parallel finite element spaces */
    void GetParBlocks(Array2D<HypreParMatrix *> &blocks) const;
 
+   /** @brief Return in @a A a parallel (on truedofs) version of this operator. */
+   virtual void FormRectangularSystemMatrix(OperatorHandle &A);
+
    virtual ~ParDiscreteLinearOperator() { }
 };
 

fem/pfespace.cpp

@@ -101,6 +101,8 @@ void ParFiniteElementSpace::ParInit(ParMesh *pm)
 
    P = NULL;
    Pconf = NULL;
+   Rconf = NULL;
+   R_transpose = NULL;
    R = NULL;
 
    num_face_nbr_dofs = -1;
@@ -927,6 +929,45 @@ const Operator *ParFiniteElementSpace::GetProlongationMatrix() const
    }
 }
 
+const Operator *ParFiniteElementSpace::GetRestrictionOperator() const
+{
+   if (Conforming())
+   {
+      if (Rconf) { return Rconf; }
+
+      if (NRanks == 1)
+      {
+         R_transpose = new IdentityOperator(GetTrueVSize());
+      }
+      else
+      {
+         if (!Device::Allows(Backend::DEVICE_MASK))
+         {
+            R_transpose = new ConformingProlongationOperator(*this, true);
+         }
+         else
+         {
+            R_transpose =
+               new DeviceConformingProlongationOperator(*this, true);
+         }
+      }
+      Rconf = new TransposeOperator(R_transpose);
+      return Rconf;
+   }
+   else
+   {
+      Dof_TrueDof_Matrix();
+      R_transpose = new TransposeOperator(R);
+      return R;
+   }
+}
+
+const Operator *ParFiniteElementSpace::GetRestrictionTransposeOperator() const
+{
+   GetRestrictionOperator();
+   return R_transpose;
+}
+
 void ParFiniteElementSpace::ExchangeFaceNbrData()
 {
    if (num_face_nbr_dofs >= 0) { return; }
@@ -2840,6 +2881,8 @@ void ParFiniteElementSpace::Destroy()
 
    delete P; P = NULL;
    delete Pconf; Pconf = NULL;
+   delete Rconf; Rconf = NULL;
+   delete R_transpose; R_transpose = NULL;
    delete R; R = NULL;
 
    delete gcomm; gcomm = NULL;
@@ -2965,12 +3008,12 @@ void ParFiniteElementSpace::Update(bool want_transform)
    }
 }
 
-
 ConformingProlongationOperator::ConformingProlongationOperator(
-   const ParFiniteElementSpace &pfes)
+   const ParFiniteElementSpace &pfes, bool local_)
    : Operator(pfes.GetVSize(), pfes.GetTrueVSize()),
      external_ldofs(),
-     gc(pfes.GroupComm())
+     gc(pfes.GroupComm()),
+     local(local_)
 {
    MFEM_VERIFY(pfes.Conforming(), "");
    const Table &group_ldof = gc.GroupLDofTable();
@@ -3019,7 +3062,14 @@ void ConformingProlongationOperator::Mult(const Vector &x, Vector &y) const
    const int m = external_ldofs.Size();
 
    const int in_layout = 2; // 2 - input is ltdofs array
-   gc.BcastBegin(const_cast<double*>(xdata), in_layout);
+   if (local)
+   {
+      y = 0.0;
+   }
+   else
+   {
+      gc.BcastBegin(const_cast<double*>(xdata), in_layout);
+   }
 
    int j = 0;
    for (int i = 0; i < m; i++)
@@ -3031,7 +3081,10 @@ void ConformingProlongationOperator::Mult(const Vector &x, Vector &y) const
    std::copy(xdata+j-m, xdata+Width(), ydata+j);
 
    const int out_layout = 0; // 0 - output is ldofs array
-   gc.BcastEnd(ydata, out_layout);
+   if (!local)
+   {
+      gc.BcastEnd(ydata, out_layout);
+   }
 }
 
 void ConformingProlongationOperator::MultTranspose(
@@ -3044,7 +3097,10 @@ void ConformingProlongationOperator::MultTranspose(
    double *ydata = y.HostWrite();
    const int m = external_ldofs.Size();
 
-   gc.ReduceBegin(xdata);
+   if (!local)
+   {
+      gc.ReduceBegin(xdata);
+   }
 
    int j = 0;
    for (int i = 0; i < m; i++)
@@ -3056,13 +3112,18 @@ void ConformingProlongationOperator::MultTranspose(
    std::copy(xdata+j, xdata+Height(), ydata+j-m);
 
    const int out_layout = 2; // 2 - output is an array on all ltdofs
-   gc.ReduceEnd<double>(ydata, out_layout, GroupCommunicator::Sum);
+   if (!local)
+   {
+      gc.ReduceEnd<double>(ydata, out_layout, GroupCommunicator::Sum);
+   }
 }
 
 DeviceConformingProlongationOperator::DeviceConformingProlongationOperator(
-   const ParFiniteElementSpace &pfes) :
+   const ParFiniteElementSpace &pfes,
+   bool local_) :
    ConformingProlongationOperator(pfes),
-   mpi_gpu_aware(Device::GetGPUAwareMPI())
+   mpi_gpu_aware(Device::GetGPUAwareMPI()),
+   local(local_)
 {
    MFEM_ASSERT(pfes.Conforming(), "internal error");
    const SparseMatrix *R = pfes.GetRestrictionMatrix();
@@ -3179,32 +3240,42 @@ void DeviceConformingProlongationOperator::Mult(const Vector &x,
                                                 Vector &y) const
 {
    const GroupTopology &gtopo = gc.GetGroupTopology();
-   BcastBeginCopy(x); // copy to 'shr_buf'
    int req_counter = 0;
-   for (int nbr = 1; nbr < gtopo.GetNumNeighbors(); nbr++)
+   if (local)
    {
-      const int send_offset = shr_buf_offsets[nbr];
-      const int send_size = shr_buf_offsets[nbr+1] - send_offset;
-      if (send_size > 0)
+      y = 0.0;
+   }
+   else
+   {
+      BcastBeginCopy(x); // copy to 'shr_buf'
+      for (int nbr = 1; nbr < gtopo.GetNumNeighbors(); nbr++)
       {
-         auto send_buf = mpi_gpu_aware ? shr_buf.Read() : shr_buf.HostRead();
-         MPI_Isend(send_buf + send_offset, send_size, MPI_DOUBLE,
-                   gtopo.GetNeighborRank(nbr), 41822,
-                   gtopo.GetComm(), &requests[req_counter++]);
-      }
-      const int recv_offset = ext_buf_offsets[nbr];
-      const int recv_size = ext_buf_offsets[nbr+1] - recv_offset;
-      if (recv_size > 0)
-      {
-         auto recv_buf = mpi_gpu_aware ? ext_buf.Write() : ext_buf.HostWrite();
-         MPI_Irecv(recv_buf + recv_offset, recv_size, MPI_DOUBLE,
-                   gtopo.GetNeighborRank(nbr), 41822,
-                   gtopo.GetComm(), &requests[req_counter++]);
+         const int send_offset = shr_buf_offsets[nbr];
+         const int send_size = shr_buf_offsets[nbr+1] - send_offset;
+         if (send_size > 0)
+         {
+            auto send_buf = mpi_gpu_aware ? shr_buf.Read() : shr_buf.HostRead();
+            MPI_Isend(send_buf + send_offset, send_size, MPI_DOUBLE,
+                      gtopo.GetNeighborRank(nbr), 41822,
+                      gtopo.GetComm(), &requests[req_counter++]);
+         }
+         const int recv_offset = ext_buf_offsets[nbr];
+         const int recv_size = ext_buf_offsets[nbr+1] - recv_offset;
+         if (recv_size > 0)
+         {
+            auto recv_buf = mpi_gpu_aware ? ext_buf.Write() : ext_buf.HostWrite();
+            MPI_Irecv(recv_buf + recv_offset, recv_size, MPI_DOUBLE,
+                      gtopo.GetNeighborRank(nbr), 41822,
+                      gtopo.GetComm(), &requests[req_counter++]);
+         }
       }
    }
    BcastLocalCopy(x, y);
-   MPI_Waitall(req_counter, requests, MPI_STATUSES_IGNORE);
-   BcastEndCopy(y); // copy from 'ext_buf'
+   if (!local)
+   {
+      MPI_Waitall(req_counter, requests, MPI_STATUSES_IGNORE);
+      BcastEndCopy(y); // copy from 'ext_buf'
+   }
 }
 
 DeviceConformingProlongationOperator::~DeviceConformingProlongationOperator()
@@ -3267,32 +3338,38 @@ void DeviceConformingProlongationOperator::MultTranspose(const Vector &x,
                                                          Vector &y) const
 {
    const GroupTopology &gtopo = gc.GetGroupTopology();
-   ReduceBeginCopy(x); // copy to 'ext_buf'
    int req_counter = 0;
-   for (int nbr = 1; nbr < gtopo.GetNumNeighbors(); nbr++)
+   if (!local)
    {
-      const int send_offset = ext_buf_offsets[nbr];
-      const int send_size = ext_buf_offsets[nbr+1] - send_offset;
-      if (send_size > 0)
+      ReduceBeginCopy(x); // copy to 'ext_buf'
+      for (int nbr = 1; nbr < gtopo.GetNumNeighbors(); nbr++)
       {
-         auto send_buf = mpi_gpu_aware ? ext_buf.Read() : ext_buf.HostRead();
-         MPI_Isend(send_buf + send_offset, send_size, MPI_DOUBLE,
-                   gtopo.GetNeighborRank(nbr), 41823,
-                   gtopo.GetComm(), &requests[req_counter++]);
-      }
-      const int recv_offset = shr_buf_offsets[nbr];
-      const int recv_size = shr_buf_offsets[nbr+1] - recv_offset;
-      if (recv_size > 0)
-      {
-         auto recv_buf = mpi_gpu_aware ? shr_buf.Write() : shr_buf.HostWrite();
-         MPI_Irecv(recv_buf + recv_offset, recv_size, MPI_DOUBLE,
-                   gtopo.GetNeighborRank(nbr), 41823,
-                   gtopo.GetComm(), &requests[req_counter++]);
+         const int send_offset = ext_buf_offsets[nbr];
+         const int send_size = ext_buf_offsets[nbr+1] - send_offset;
+         if (send_size > 0)
+         {
+            auto send_buf = mpi_gpu_aware ? ext_buf.Read() : ext_buf.HostRead();
+            MPI_Isend(send_buf + send_offset, send_size, MPI_DOUBLE,
+                      gtopo.GetNeighborRank(nbr), 41823,
+                      gtopo.GetComm(), &requests[req_counter++]);
+         }
+         const int recv_offset = shr_buf_offsets[nbr];
+         const int recv_size = shr_buf_offsets[nbr+1] - recv_offset;
+         if (recv_size > 0)
+         {
+            auto recv_buf = mpi_gpu_aware ? shr_buf.Write() : shr_buf.HostWrite();
+            MPI_Irecv(recv_buf + recv_offset, recv_size, MPI_DOUBLE,
+                      gtopo.GetNeighborRank(nbr), 41823,
+                      gtopo.GetComm(), &requests[req_counter++]);
+         }
       }
    }
    ReduceLocalCopy(x, y);
-   MPI_Waitall(req_counter, requests, MPI_STATUSES_IGNORE);
-   ReduceEndAssemble(y); // assemble from 'shr_buf'
+   if (!local)
+   {
+      MPI_Waitall(req_counter, requests, MPI_STATUSES_IGNORE);
+      ReduceEndAssemble(y); // assemble from 'shr_buf'
+   }
 }
 
 } // namespace mfem

fem/pfespace.hpp

@@ -75,6 +75,12 @@ private:
 
    /// The (block-diagonal) matrix R (restriction of dof to true dof). Owned.
    mutable SparseMatrix *R;
+   /// Optimized action-only restriction operator for conforming meshes. Owned.
+   mutable Operator *Rconf;
+   /** Transpose of R or Rconf. For conforming mesh, this is a matrix-free
+       (Device)ConformingProlongationOperator, for a non-conforming mesh
+       this is a TransposeOperator wrapping R. */
+   mutable Operator *R_transpose;
 
    ParNURBSExtension *pNURBSext() const
    { return dynamic_cast<ParNURBSExtension *>(NURBSext); }
@@ -341,6 +347,16 @@ public:
    HYPRE_Int GetMyTDofOffset() const;
 
    virtual const Operator *GetProlongationMatrix() const;
+   /** @brief Return logical transpose of restriction matrix, but in
+       non-assembled optimized matrix-free form.
+
+       The implementation is like GetProlongationMatrix, but it sets local
+       DOFs to the true DOF values if owned locally, otherwise zero. */
+   virtual const Operator *GetRestrictionTransposeOperator() const;
+   /** Get an Operator that performs the action of GetRestrictionMatrix(),
+       but potentially with a non-assembled optimized matrix-free
+       implementation. */
+   virtual const Operator *GetRestrictionOperator() const;
    /// Get the R matrix which restricts a local dof vector to true dof vector.
    virtual const SparseMatrix *GetRestrictionMatrix() const
    { Dof_TrueDof_Matrix(); return R; }
@@ -395,9 +411,11 @@ class ConformingProlongationOperator : public Operator
 protected:
    Array<int> external_ldofs;
    const GroupCommunicator &gc;
+   bool local;
 
 public:
-   ConformingProlongationOperator(const ParFiniteElementSpace &pfes);
+   ConformingProlongationOperator(const ParFiniteElementSpace &pfes,
+                                  bool local_=false);
 
    virtual void Mult(const Vector &x, Vector &y) const;
 
@@ -416,6 +434,8 @@ protected:
    Array<int> ltdof_ldof, unq_ltdof;
    Array<int> unq_shr_i, unq_shr_j;
    MPI_Request *requests;
+   bool local;
+
    // Kernel: copy ltdofs from 'src' to 'shr_buf' - prepare for send.
    //         shr_buf[i] = src[shr_ltdof[i]]
    void BcastBeginCopy(const Vector &src) const;
@@ -441,7 +461,8 @@ protected:
    void ReduceEndAssemble(Vector &dst) const;
 
 public:
-   DeviceConformingProlongationOperator(const ParFiniteElementSpace &pfes);
+   DeviceConformingProlongationOperator(const ParFiniteElementSpace &pfes,
+                                        bool local_=false);
 
    virtual ~DeviceConformingProlongationOperator();
 

fem/restriction.cpp

@@ -195,6 +195,31 @@ void ElementRestriction::MultTransposeUnsigned(const Vector& x, Vector& y) const
    });
 }
 
+void ElementRestriction::MultLeftInverse(const Vector& x, Vector& y) const
+{
+   // Assumes all elements have the same number of dofs
+   const int nd = dof;
+   const int vd = vdim;
+   const bool t = byvdim;
+   auto d_offsets = offsets.Read();
+   auto d_indices = indices.Read();
+   auto d_x = Reshape(x.Read(), nd, vd, ne);
+   auto d_y = Reshape(y.Write(), t?vd:ndofs, t?ndofs:vd);
+   MFEM_FORALL(i, ndofs,
+   {
+      const int nextOffset = d_offsets[i + 1];
+      for (int c = 0; c < vd; ++c)
+      {
+         double dofValue = 0;
+         const int j = nextOffset - 1;
+         const int idx_j = (d_indices[j] >= 0) ? d_indices[j] : -1 - d_indices[j];
+         dofValue = (d_indices[j] >= 0) ? d_x(idx_j % nd, c, idx_j / nd) :
+         -d_x(idx_j % nd, c, idx_j / nd);
+         d_y(t?c:i,t?i:c) = dofValue;
+      }
+   });
+}
+
 void ElementRestriction::BooleanMask(Vector& y) const
 {
    // Assumes all elements have the same number of dofs

fem/restriction.hpp

@@ -57,6 +57,10 @@ public:
    /// Compute MultTranspose without applying signs based on DOF orientations.
    void MultTransposeUnsigned(const Vector &x, Vector &y) const;
 
+   /// Compute MultTranspose by setting (rather than adding) element
+   /// contributions; this is a left inverse of the Mult() operation
+   void MultLeftInverse(const Vector &x, Vector &y) const;
+
    /// @brief Fills the E-vector y with `boolean` values 0.0 and 1.0 such that each
    /// each entry of the L-vector is uniquely represented in `y`.
    /** This means, the sum of the E-vector `y` is equal to the sum of the

linalg/CMakeLists.txt

@@ -10,6 +10,7 @@
 # CONTRIBUTING.md for details.
 
 list(APPEND SRCS
+  auxiliary.cpp
   blockmatrix.cpp
   blockoperator.cpp
   blockvector.cpp
@@ -27,6 +28,7 @@ list(APPEND SRCS
   )
 
 list(APPEND HDRS
+  auxiliary.hpp
   blockmatrix.hpp
   blockoperator.hpp
   blockvector.hpp

linalg/auxiliary.hpp

@@ -0,0 +1,285 @@
+// 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_AUXILIARY
+#define MFEM_AUXILIARY
+
+#include "../config/config.hpp"
+
+#ifdef MFEM_USE_MPI
+
+#include "../general/tic_toc.hpp"
+#include "solvers.hpp"
+
+namespace mfem
+{
+
+// forward declarations
+class Coefficient;
+class MatrixCoefficient;
+class ParMesh;
+class ParBilinearForm;
+class ParDiscreteLinearOperator;
+
+/** @brief Auxiliary space solvers for MatrixFreeAMS preconditioner
+
+    Given an operator A and a transfer G, this will create a solver
+    that approximates (G^T A G)^{-1}. Used for two different
+    auxiliary spaces in the AMS cycle.
+
+    The produced solver is based on a low-order refined discretization
+    for the high-order H1 problem. */
+class MatrixFreeAuxiliarySpace : public Solver
+{
+public:
+   /** @brief Pi space constructor
+
+       In the AMS framework this auxiliary space has two coefficients.
+
+       @param mesh_lor       Low-order refined auxiliary mesh
+       @param alpha_coeff    coefficient on curl-curl term (1 if null)
+       @param beta_coeff     coefficient on mass term (1 if null)
+       @param beta_mcoeff    matrix coefficient on mass term
+       @param ess_bdr        attributes for essential boundaries
+       @param curlcurl_oper  High-order operator for the system
+       @param pi             Intentity interpolation operator
+       @param useAmgX_       Use AmgX instead of hypre for auxiliary solves
+       @param cg_iterations  number of CG iterations used to invert
+                             auxiliary system, choosing 0 means to use a
+                             single V-cycle
+   */
+   MatrixFreeAuxiliarySpace(
+      ParMesh& mesh_lor, Coefficient* alpha_coeff, Coefficient* beta_coeff,
+      MatrixCoefficient* beta_mcoeff,
+      Array<int>& ess_bdr, Operator& curlcurl_oper, Operator& pi,
+#ifdef MFEM_USE_AMGX
+      bool useAmgX_,
+#endif
+      int cg_iterations = 0);
+
+   // Complex Pi space constructor
+   MatrixFreeAuxiliarySpace(
+      ParMesh& mesh_lor, Coefficient* alpha_coeff, Coefficient* beta_coeff,
+      Coefficient* beta_imag, Coefficient* abs_beta_imag,
+      MatrixCoefficient* beta_mcoeff,
+      Array<int>& ess_bdr, Operator& curlcurl_oper, Operator *oper_complex,
+      Operator& pi,
+#ifdef MFEM_USE_AMGX
+      bool useAmgX_,
+#endif
+      int cg_iterations = 0);
+
+   /** @brief G space constructor
+
+       This has one coefficient in the AMS framework.
+
+       @param mesh_lor       Low-order refined auxiliary mesh
+       @param beta_coeff     coefficient on mass term (1 if null)
+       @param beta_mcoeff    matrix coefficient on mass term
+       @param ess_bdr        attributes for essential boundaries
+       @param curlcurl_oper  High-order operator for the system
+       @param g              Gradient interpolation operator
+       @param useAmgX_       Use AmgX instead of hypre for auxiliary solves
+       @param cg_iterations  number of CG iterations used to invert
+                             auxiliary system, choosing 0 means to
+                             use a single V-cycle
+   */
+   MatrixFreeAuxiliarySpace(
+      ParMesh& mesh_lor, Coefficient* beta_coeff,
+      MatrixCoefficient* beta_mcoeff, Array<int>& ess_bdr,
+      Operator& curlcurl_oper, Operator& g,
+#ifdef MFEM_USE_AMGX
+      bool useAmgX_,
+#endif
+      int cg_iterations = 1);
+
+   // Complex G space constructor
+   MatrixFreeAuxiliarySpace(
+      ParMesh& mesh_lor, Coefficient* beta_coeff, Coefficient* beta_imag,
+      Coefficient* abs_beta_imag,
+      MatrixCoefficient* beta_mcoeff, Array<int>& ess_bdr,
+      Operator& curlcurl_oper, Operator *oper_complex, Operator& g,
+#ifdef MFEM_USE_AMGX
+      bool useAmgX_,
+#endif
+      int cg_iterations = 1);
+
+   ~MatrixFreeAuxiliarySpace();
+
+   void Mult(const Vector& x, Vector& y) const;
+
+   void SetOperator(const Operator& op) {}
+
+private:
+   /** @brief Helper routine for constructors.
+
+       @param system_dimension is passed to HypreBoomerAMG::SetSystemsOptions
+   */
+   void SetupAMG(int system_dimension);
+   void SetupVCycle();
+
+   /// inner_cg_iterations > 99 applies an exact solve here
+   void SetupCG(Operator& curlcurl_oper, Operator& conn,
+                int inner_cg_iterations);
+
+   void SetupGMRES(Operator& curlcurl_oper, Operator& conn);
+
+   void SetupPMHSS();
+
+   MPI_Comm comm;
+   Array<int> ess_tdof_list;
+   HypreParMatrix * aspacematrix;
+   HypreParMatrix * aspacematrix_complex;
+   HypreParMatrix * aspacematrix_imag;
+   Solver * aspacepc;
+   Operator* matfree;
+   CGSolver* cg;
+   GMRESSolver* gmres;
+   GMRESSolver* gmres_PMHSS;
+   Operator* aspacewrapper;
+#ifdef MFEM_USE_AMGX
+   const bool useAmgX;
+#endif
+   mutable int inner_aux_iterations;
+
+   const bool imagBdry;
+
+   Complex_PMHSS *PMHSS = NULL;
+
+   Array<int> offsets;
+   Array<int> offsets_nd;
+   BlockDiagonalPreconditioner *BlockDP;
+
+   BlockOperator *conn_block;
+};
+
+
+/** @brief Perform AMS cycle with generic Operator objects.
+
+    Most users should use MatrixFreeAMS, which wraps this. */
+class GeneralAMS : public Solver
+{
+public:
+   /** @brief Constructor.
+
+       Most of these arguments just need a Mult() operation,
+       but pi and g also require MultTranspose() */
+   GeneralAMS(const Operator& curlcurl_op_,
+              Operator *oper_complex,
+              const Operator& pi_,
+              const Operator& gradient_,
+              const Operator& pispacesolver_,
+              const Operator& gspacesolver_,
+              const Operator& smoother_,
+              const Array<int>& ess_tdof_list_);
+   virtual ~GeneralAMS();
+
+   /// in principle this should set A_ = op;
+   void SetOperator(const Operator &op) {}
+
+   virtual void Mult(const Vector& x, Vector& y) const;
+
+private:
+   const Operator& curlcurl_op;
+   Operator *oper_complex;
+   const Operator& pi;
+   const Operator& gradient;
+   const Operator& pispacesolver;
+   const Operator& gspacesolver;
+   const Operator& smoother;
+   const Array<int> ess_tdof_list;
+
+   void FormResidual(const Vector& rhs, const Vector& x,
+                     Vector& residual) const;
+};
+
+
+/** @brief An auxiliary Maxwell solver for a high-order curl-curl
+    system without high-order assembly.
+
+    The auxiliary space solves are done using a low-order refined approach,
+    but all the interpolation operators, residuals, etc. are done in a
+    matrix-free manner.
+
+    See Barker and Kolev, Matrix-free preconditioning for high-order H(curl)
+    discretizations (https://doi.org/10.1002/nla.2348) */
+class MatrixFreeAMS : public Solver
+{
+public:
+   /** @brief Construct matrix-free AMS preconditioner
+
+       @param aform        BilinearForm for curl-curl problem, generally will
+                           have a CurlCurlIntegrator and possibly a
+                           VectorFEMassIntegrator.
+       @param oper         Operator to precondition.
+       @param nd_fespace   Underlying Nedelec finite element space.
+       @param alpha_coeff  coefficient on curl-curl term in Maxwell problem
+                           (can be null, in which case constant 1 is assumed)
+       @param beta_coeff   (scalar) coefficient on mass term in Maxwell problem
+       @param beta_mcoeff  (matrix) coefficient on mass term
+       @param ess_bdr      boundary *attributes* that are marked essential. In
+                           contrast to other MFEM cases, these are *attributes*
+                           not dofs, because we need to apply these boundary
+                           conditions to different bilinear forms.
+       @param useAmgX      use AmgX (instead of hypre) for LOR problems
+       @param inner_pi_its number of CG iterations on auxiliary pi space,
+                           may need more for difficult coefficients
+       @param inner_g_its  number of CG iterations on auxiliary g space,
+                           may need more for difficult coefficients
+       @param nd_smoother  optional user-provided smoother for Nedelec space,
+                           this object takes ownership and will delete.
+   */
+   MatrixFreeAMS(ParBilinearForm& aform, Operator& oper, Operator *oper_complex,
+                 ParFiniteElementSpace& nd_fespace, Coefficient* alpha_coeff,
+                 Coefficient* beta_coeff, Coefficient* beta_imag,
+                 Coefficient* abs_beta_imag, MatrixCoefficient* beta_mcoeff,
+                 Array<int>& ess_bdr,
+#ifdef MFEM_USE_AMGX
+                 bool useAmgX = false,
+#endif
+                 int inner_pi_its = 0, int inner_g_its = 1,
+                 Solver* nd_smoother = NULL);
+
+   ~MatrixFreeAMS();
+
+   void SetOperator(const Operator &op) {}
+
+   void Mult(const Vector& x, Vector& y) const { general_ams->Mult(x, y); }
+
+private:
+   GeneralAMS * general_ams;
+
+   Solver * smoother;
+   ParDiscreteLinearOperator * pa_grad;
+   OperatorPtr Gradient;
+   ParDiscreteLinearOperator * pa_interp;
+   OperatorPtr Pi;
+
+   Solver * Gspacesolver;
+   Solver * Pispacesolver;
+
+   ParFiniteElementSpace * h1_fespace;
+   ParFiniteElementSpace * h1_fespace_d;
+
+   Array<int> offsets_nd;
+   Array<int> offsets_vector;
+   Array<int> offsets_scalar;
+
+   BlockOperator *Pi_block;
+   BlockOperator *Gradient_block;
+   BlockOperator *smoother_block;
+};
+
+} // namespace mfem
+
+#endif // MFEM_USE_MPI
+
+#endif

linalg/linalg.hpp

@@ -29,6 +29,7 @@
 #include "solvers.hpp"
 #include "handle.hpp"
 #include "invariants.hpp"
+#include "auxiliary.hpp"
 
 #ifdef MFEM_USE_AMGX
 #include "amgxsolver.hpp"

linalg/operator.hpp

@@ -37,7 +37,8 @@ protected:
       const Array<int> &test_tdof_list,
       RectangularConstrainedOperator* &Aout);
 
-   /// Returns RAP Operator of this, taking in input/output Prolongation matrices
+   /** @brief Returns RAP Operator of this, using input/output Prolongation matrices
+       @a Pi corresponds to "P", @a Po corresponds to "Rt" */
    Operator *SetupRAP(const Operator *Pi, const Operator *Po);
 
 public:
@@ -112,6 +113,11 @@ public:
    {
       return GetProlongation(); // Assume square unless specialized
    }
+   /** @brief Transpose of GetOutputRestriction, directly available in this
+       form to facilitate matrix-free RAP-type operators.
+
+       `NULL` means identity. */
+   virtual const Operator *GetOutputRestrictionTranspose() const { return NULL; }
    /** @brief Restriction operator from output vectors for the operator to linear
        algebra (linear system) vectors. `NULL` means identity. */
    virtual const Operator *GetOutputRestriction() const
@@ -690,6 +696,42 @@ public:
    { A_.Mult(x, y); y *= a_; }
 };
 
+/// General sum operator: x -> A(x)+B(x)
+class SumOperator : public Operator
+{
+   const Operator *A, *B;
+   bool ownA, ownB;
+   mutable Vector z, w;
+   double cA, cB;
+
+public:
+   SumOperator(const Operator *A_, const Operator *B_,
+               bool ownA_, bool ownB_, double cA_, double cB_)
+      : Operator(A_->Height(), B_->Width()),
+        A(A_), B(B_), ownA(ownA_), ownB(ownB_), z(A_->Height()), w(A_->Width()),
+        cA(cA_), cB(cB_)
+   {
+      MFEM_VERIFY(A->Width() == B->Width() && A->Height() == B->Height(),
+                  "incompatible Operators: A->Width() = " << A->Width()
+                  << ", B->Height() = " << B->Height());
+
+      z.UseDevice(true);
+      w.UseDevice(true);
+   }
+
+   ~SumOperator()
+   {
+      if (ownA) { delete A; }
+      if (ownB) { delete B; }
+   }
+
+   virtual void Mult(const Vector &x, Vector &y) const
+   { B->Mult(x, z); A->Mult(x, y); y *= cA; z *= cB; y += z;}
+
+   virtual void MultTranspose(const Vector &x, Vector &y) const
+   { B->MultTranspose(x, w); A->MultTranspose(x, y); y *= cA; w *= cB; y += w;}
+
+};
 
 /** @brief The transpose of a given operator. Switches the roles of the methods
     Mult() and MultTranspose(). */

linalg/solvers.hpp

@@ -852,6 +852,169 @@ public:
 
 #endif // MFEM_USE_SUITESPARSE
 
+class Complex_PMHSS : public Solver
+{
+public:
+   Complex_PMHSS(Operator *Re, Operator *Im, Solver *prec_Re, Solver *prec_Im,
+                 double a_)
+      : Solver(2*Re->Height()), a(a_), A(Re, Im, false, false),
+        A_Re(Re, NULL, false, false),
+        A_Im(Im, NULL, false, false), u(2*Re->Height()), rhs(2*Re->Height()),
+        n(Re->Height())
+   {
+      MFEM_VERIFY(Re->Height() == Im->Height() && Re->Height() == Re->Width() &&
+                  Im->Height() == Im->Width(), "");
+      MFEM_VERIFY(this->Height() == A.Height(), "");
+
+      // Create CG solver for real operator aV + A_Re in complex space.
+
+      V = useIdentityV ? (Operator*) new IdentityOperator(this->Height()) :
+          (Operator*) &A_Re;
+
+      // In the case V = A_Re, it is faster to use a scaled operator than a SumOperator
+      Operator *sumOpRe = useIdentityV ? (Operator*) new SumOperator(V, &A_Re, false,
+                                                                     false, a, 1.0)
+                          : (Operator*) new ScaledOperator(&A_Re, a + 1.0);
+
+      SumOperator *sumOpIm = new SumOperator(V, &A_Im, false, false, a, 1.0);
+
+      CGSolver *cg = new CGSolver(MPI_COMM_WORLD);
+      cg->SetRelTol(1e-6);
+      cg->SetMaxIter(1000);
+      cg->SetPrintLevel(0);
+      cg->SetOperator(*sumOpRe);
+      cg->SetPreconditioner(*prec_Re);
+      cg->iterative_mode = false;
+
+      SRe = cg;
+
+      CGSolver *cgi = new CGSolver(MPI_COMM_WORLD);
+      cgi->SetRelTol(1e-6);
+      cgi->SetMaxIter(1000);
+      cgi->SetPrintLevel(0);
+      cgi->SetOperator(*sumOpIm);
+      if (prec_Im && useIdentityV) { cgi->SetPreconditioner(*prec_Im); }
+      if (!useIdentityV) { cgi->SetPreconditioner(*prec_Re); }
+      cgi->iterative_mode = false;
+
+      /*
+      // For negative definite imaginary part, but then PMHSS does not work?
+      MINRESSolver *cgi = new MINRESSolver(MPI_COMM_WORLD);
+      cgi->SetRelTol(1e-12);
+      cgi->SetMaxIter(1000);
+      cgi->SetPrintLevel(0);
+      cgi->SetOperator(*sumOpIm);
+      if (prec_Im) cgi->SetPreconditioner(*prec_Im);
+      */
+
+      SIm = cgi;
+   }
+
+   void SetOperator(const Operator &op)
+   {
+      MFEM_VERIFY(false, "Don't call SetOperator");
+   }
+
+   void ComputeResidual(const Vector &b, const Vector &sol, Vector &res) const
+   {
+      A.Mult(sol, res);
+      res -= b;
+   }
+
+   void Mult(const Vector &x, Vector &y) const
+   {
+      if (!(x.Size() == Height() && y.Size() == Height()))
+      {
+         std::cout << "bug";
+      }
+
+      MFEM_VERIFY(x.Size() == Height() && y.Size() == Height(), "");
+
+      const double initNorm = x.Norml2();
+      mfem::out << "MHSS RHS norm " << initNorm << '\n';
+
+      // With V = I, use modified HSS (MHSS) from Bai, Benzi, Chen 2010.
+      y = 0.0;
+
+      for (int it=0; it<maxiter; ++it)
+      {
+         // Solve (aI + Re) u = (aI - i Im) y + x
+
+         if (it == 0)
+         {
+            // Optimize the first iteration, when the initial guess is y=0.
+            SRe->Mult(x, u);
+         }
+         else
+         {
+            A_Im.Mult(y, u);  // u = Im y
+            // Set rhs = -i Im y = -i u
+            for (int j=0; j<n; ++j)
+            {
+               rhs[j] = u[n+j];
+               rhs[n+j] = -u[j];
+            }
+
+            rhs += x;
+
+            V->Mult(y, u);
+            rhs.Add(a, u);
+
+            SRe->Mult(rhs, u);
+         }
+
+         // Solve (aI + Im) y = (aI + i Re) u - i x
+
+         A_Re.Mult(u, y);  // y = Re u
+         // Set rhs = i (Re u - x) = i (y - x)
+         for (int j=0; j<n; ++j)
+         {
+            rhs[j] = -(y[n+j] - x[n+j]);
+            rhs[n+j] = y[j] - x[j];
+         }
+
+         if (useIdentityV)
+         {
+            //V->Mult(u, y);
+            //rhs.Add(a, y);
+            rhs.Add(a, u);
+         }
+         else
+         {
+            // Using V = A_Re
+            rhs.Add(a, y);
+         }
+
+         SIm->Mult(rhs, y);
+
+         ComputeResidual(x, y, rhs);
+         const double resNorm = rhs.Norml2();
+         mfem::out << "MHSS iter " << it << " residual norm " << resNorm << '\n';
+
+         if (resNorm / initNorm < tol)
+         {
+            mfem::out << "MHSS converged\n";
+            break;
+         }
+      }
+   }
+
+private:
+   const double a;
+   const int maxiter = 1;
+   ComplexOperator A, A_Re, A_Im;
+   mutable Vector u, rhs;
+   const int n;
+
+   const double tol = 1.0e-8;
+
+   const bool useIdentityV = false;
+   Operator *V = NULL;
+
+   Solver *SRe = NULL;
+   Solver *SIm = NULL;
+};
+
 }
 
 #endif // MFEM_SOLVERS

tests/unit/CMakeLists.txt

@@ -55,6 +55,8 @@ set(UNIT_TESTS_SRCS
   fem/test_operatorjacobismoother.cpp
   fem/test_pa_coeff.cpp
   fem/test_pa_kernels.cpp
+  fem/test_pa_grad.cpp
+  fem/test_pa_idinterp.cpp
   fem/test_quadf_coef.cpp
   fem/test_quadraturefunc.cpp
   fem/test_blocknonlinearform.cpp

tests/unit/fem/test_pa_grad.cpp

@@ -0,0 +1,219 @@
+// 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 "catch.hpp"
+#include "mfem.hpp"
+#include "unit_tests.hpp"
+
+using namespace mfem;
+
+double compare_pa_assembly(int dim, int num_elements, int order, bool transpose)
+{
+   Mesh * mesh;
+   if (num_elements == 0)
+   {
+      if (dim == 2)
+      {
+         mesh = new Mesh("../../data/star.mesh", order);
+      }
+      else
+      {
+         mesh = new Mesh("../../data/beam-hex.mesh", order);
+      }
+   }
+   else
+   {
+      if (dim == 2)
+      {
+         mesh = new Mesh(num_elements, num_elements, Element::QUADRILATERAL, true);
+      }
+      else
+      {
+         mesh = new Mesh(num_elements, num_elements, num_elements,
+                         Element::HEXAHEDRON, true);
+      }
+   }
+   FiniteElementCollection *h1_fec = new H1_FECollection(order, dim);
+   FiniteElementCollection *nd_fec = new ND_FECollection(order, dim);
+   FiniteElementSpace h1_fespace(mesh, h1_fec);
+   FiniteElementSpace nd_fespace(mesh, nd_fec);
+
+   DiscreteLinearOperator assembled_grad(&h1_fespace, &nd_fespace);
+   assembled_grad.AddDomainInterpolator(new GradientInterpolator);
+   const int skip_zeros = 1;
+   assembled_grad.Assemble(skip_zeros);
+   assembled_grad.Finalize(skip_zeros);
+   const SparseMatrix& assembled_grad_mat = assembled_grad.SpMat();
+
+   DiscreteLinearOperator pa_grad(&h1_fespace, &nd_fespace);
+   pa_grad.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+   pa_grad.AddDomainInterpolator(new GradientInterpolator);
+   pa_grad.Assemble();
+   pa_grad.Finalize();
+
+   int insize, outsize;
+   if (transpose)
+   {
+      insize = nd_fespace.GetVSize();
+      outsize = h1_fespace.GetVSize();
+   }
+   else
+   {
+      insize = h1_fespace.GetVSize();
+      outsize = nd_fespace.GetVSize();
+   }
+   Vector xv(insize);
+   Vector assembled_y(outsize);
+   Vector pa_y(outsize);
+
+   xv.Randomize();
+   if (transpose)
+   {
+      assembled_grad_mat.BuildTranspose();
+      assembled_grad_mat.MultTranspose(xv, assembled_y);
+      pa_grad.MultTranspose(xv, pa_y);
+   }
+   else
+   {
+      assembled_grad_mat.Mult(xv, assembled_y);
+      pa_grad.Mult(xv, pa_y);
+   }
+
+   pa_y -= assembled_y;
+   double error = pa_y.Norml2() / assembled_y.Norml2();
+   INFO("dim " << dim << " ne " << num_elements << " order " << order
+        << (transpose ? " T:" : ":") << " error in PA gradient: " << error);
+
+   delete h1_fec;
+   delete nd_fec;
+   delete mesh;
+
+   return error;
+}
+
+TEST_CASE("PAGradient", "[CUDA]")
+{
+   auto transpose = GENERATE(true, false);
+   auto order = GENERATE(1, 2, 3, 4);
+   auto dim = GENERATE(2, 3);
+   auto num_elements = GENERATE(0, 1, 2, 3, 4);
+
+   double error = compare_pa_assembly(dim, num_elements, order, transpose);
+   REQUIRE(error == MFEM_Approx(0.0, 1.0e-14));
+}
+
+#ifdef MFEM_USE_MPI
+
+double par_compare_pa_assembly(int dim, int num_elements, int order,
+                               bool transpose)
+{
+   int rank;
+   MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+   int size;
+   MPI_Comm_size(MPI_COMM_WORLD, &size);
+
+   Mesh * smesh;
+   if (dim == 2)
+   {
+      smesh = new Mesh(num_elements, num_elements, Element::QUADRILATERAL, true);
+   }
+   else
+   {
+      smesh = new Mesh(num_elements, num_elements, num_elements,
+                       Element::HEXAHEDRON, true);
+   }
+   ParMesh * mesh = new ParMesh(MPI_COMM_WORLD, *smesh);
+   delete smesh;
+   FiniteElementCollection *h1_fec = new H1_FECollection(order, dim);
+   FiniteElementCollection *nd_fec = new ND_FECollection(order, dim);
+   ParFiniteElementSpace h1_fespace(mesh, h1_fec);
+   ParFiniteElementSpace nd_fespace(mesh, nd_fec);
+
+   ParDiscreteLinearOperator assembled_grad(&h1_fespace, &nd_fespace);
+   assembled_grad.AddDomainInterpolator(new GradientInterpolator);
+   const int skip_zeros = 1;
+   assembled_grad.Assemble(skip_zeros);
+   assembled_grad.Finalize(skip_zeros);
+   HypreParMatrix * assembled_grad_mat = assembled_grad.ParallelAssemble();
+
+   ParDiscreteLinearOperator pa_grad(&h1_fespace, &nd_fespace);
+   pa_grad.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+   pa_grad.AddDomainInterpolator(new GradientInterpolator);
+   pa_grad.Assemble();
+   OperatorPtr pa_grad_oper;
+   pa_grad.FormRectangularSystemMatrix(pa_grad_oper);
+
+   int insize, outsize;
+   if (transpose)
+   {
+      insize = assembled_grad_mat->Height();
+      outsize = assembled_grad_mat->Width();
+   }
+   else
+   {
+      insize = assembled_grad_mat->Width();
+      outsize = assembled_grad_mat->Height();
+   }
+   Vector xv(insize);
+   Vector assembled_y(outsize);
+   Vector pa_y(outsize);
+   assembled_y = 0.0;
+   pa_y = 0.0;
+
+   xv.Randomize();
+   if (transpose)
+   {
+      assembled_grad_mat->MultTranspose(xv, assembled_y);
+      pa_grad_oper->MultTranspose(xv, pa_y);
+   }
+   else
+   {
+      assembled_grad_mat->Mult(xv, assembled_y);
+      pa_grad_oper->Mult(xv, pa_y);
+   }
+
+   Vector error_vec(pa_y);
+   error_vec -= assembled_y;
+   // serial norms and serial error; we are enforcing equality on each processor
+   // in the test
+   double error = error_vec.Norml2() / assembled_y.Norml2();
+
+   for (int p = 0; p < size; ++p)
+   {
+      if (rank == p)
+      {
+         INFO("[" << rank << "][par] dim " << dim << " ne " << num_elements
+              << " order " << order << (transpose ? " T:" : ":")
+              << " error in PA gradient: " << error);
+      }
+      MPI_Barrier(MPI_COMM_WORLD);
+   }
+
+   delete h1_fec;
+   delete nd_fec;
+   delete assembled_grad_mat;
+   delete mesh;
+
+   return error;
+}
+
+TEST_CASE("ParallelPAGradient", "[Parallel], [ParallelPAGradient]")
+{
+   auto transpose = GENERATE(true, false);
+   auto order = GENERATE(1, 2, 3, 4);
+   auto dim = GENERATE(2, 3);
+   auto num_elements = GENERATE(4, 5);
+
+   double error = par_compare_pa_assembly(dim, num_elements, order, transpose);
+   REQUIRE(error == MFEM_Approx(0.0, 1.0e-14));
+}
+
+#endif

tests/unit/fem/test_pa_idinterp.cpp

@@ -0,0 +1,123 @@
+// 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 "catch.hpp"
+#include "mfem.hpp"
+#include "unit_tests.hpp"
+
+using namespace mfem;
+
+double compare_pa_id_assembly(int dim, int num_elements, int order,
+                              bool transpose)
+{
+   Mesh * mesh;
+   if (num_elements == 0)
+   {
+      if (dim == 2)
+      {
+         mesh = new Mesh("../../data/star.mesh", order);
+      }
+      else
+      {
+         mesh = new Mesh("../../data/beam-hex.mesh", order);
+
+         // Transform mesh vertices to test without alignment with coordinate axes.
+         for (int i=0; i<mesh->GetNV(); ++i)
+         {
+            double *v = mesh->GetVertex(i);
+            const double yscale = 1.0 + v[1];
+            const double zscale = 1.0 + v[2];
+            v[0] *= zscale;
+            v[1] *= zscale;
+            v[2] *= yscale;
+         }
+      }
+   }
+   else
+   {
+      if (dim == 2)
+      {
+         mesh = new Mesh(num_elements, num_elements, Element::QUADRILATERAL, true);
+      }
+      else
+      {
+         mesh = new Mesh(num_elements, num_elements, num_elements,
+                         Element::HEXAHEDRON, true);
+      }
+   }
+   FiniteElementCollection *h1_fec = new H1_FECollection(order, dim);
+   FiniteElementCollection *nd_fec = new ND_FECollection(order, dim);
+   FiniteElementSpace h1_fespace(mesh, h1_fec, dim);
+   FiniteElementSpace nd_fespace(mesh, nd_fec);
+
+   DiscreteLinearOperator assembled_id(&h1_fespace, &nd_fespace);
+   assembled_id.AddDomainInterpolator(new IdentityInterpolator);
+   const int skip_zeros = 1;
+   assembled_id.Assemble(skip_zeros);
+   assembled_id.Finalize(skip_zeros);
+   const SparseMatrix& assembled_id_mat = assembled_id.SpMat();
+
+   DiscreteLinearOperator pa_id(&h1_fespace, &nd_fespace);
+   pa_id.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+   pa_id.AddDomainInterpolator(new IdentityInterpolator);
+   pa_id.Assemble();
+   pa_id.Finalize();
+
+   int insize, outsize;
+   if (transpose)
+   {
+      insize = nd_fespace.GetVSize();
+      outsize = h1_fespace.GetVSize();
+   }
+   else
+   {
+      insize = h1_fespace.GetVSize();
+      outsize = nd_fespace.GetVSize();
+   }
+   Vector x(insize);
+   Vector assembled_y(outsize);
+   Vector pa_y(outsize);
+
+   x.Randomize();
+   if (transpose)
+   {
+      assembled_id_mat.BuildTranspose();
+      assembled_id_mat.MultTranspose(x, assembled_y);
+      pa_id.MultTranspose(x, pa_y);
+   }
+   else
+   {
+      assembled_id.Mult(x, assembled_y);
+      pa_id.Mult(x, pa_y);
+   }
+
+   pa_y -= assembled_y;
+   double error = pa_y.Norml2() / assembled_y.Norml2();
+   INFO("dim " << dim << " ne " << num_elements << " order " << order
+        << (transpose ? " T:" : ":") << " error in PA identity: " << error);
+
+   delete h1_fec;
+   delete nd_fec;
+   delete mesh;
+
+   return error;
+}
+
+TEST_CASE("PAIdentityInterp", "[CUDA]")
+{
+   auto transpose = GENERATE(true, false);
+   auto order = GENERATE(1, 2, 3, 4);
+   auto dim = GENERATE(2, 3);
+   auto num_elements = GENERATE(0, 1, 2, 3, 4);
+
+   double error = compare_pa_id_assembly(dim, num_elements, order, transpose);
+   REQUIRE(error == MFEM_Approx(0.0, 1.0e-14));
+}