Updated explanation from PR

CHANGELOG

@@ -11,12 +11,6 @@
 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/ex16p.cpp

@@ -65,7 +65,7 @@ protected:
    HypreSmoother M_prec; // Preconditioner for the mass matrix M
 
    CGSolver T_solver;    // Implicit solver for T = M + dt K
-   HypreSmoother T_prec; // Preconditioner for the implicit solver
+   HypreBoomerAMG T_prec; // Preconditioner for the implicit solver
 
    double alpha, kappa;
 
@@ -353,7 +353,6 @@ int main(int argc, char *argv[])
          }
 #endif
       }
-      oper.SetParameters(u);
    }
 
 #ifdef MFEM_USE_ADIOS2
@@ -414,6 +413,7 @@ ConductionOperator::ConductionOperator(ParFiniteElementSpace &f, double al,
    T_solver.SetPrintLevel(0);
    T_solver.SetPreconditioner(T_prec);
 
+   T_prec.SetPrintLevel(0);
    SetParameters(u);
 }
 
@@ -430,19 +430,53 @@ void ConductionOperator::Mult(const Vector &u, Vector &du_dt) const
 void ConductionOperator::ImplicitSolve(const double dt,
                                        const Vector &u, Vector &du_dt)
 {
-   // Solve the equation:
-   //    du_dt = M^{-1}*[-K(u + dt*du_dt)]
-   // for du_dt
-   if (!T)
-   {
-      T = Add(1.0, Mmat, dt, Kmat);
-      current_dt = dt;
+   // Here we use Picard iterations to solve a nonlinear equation
+   // for the Runge-Kutta stage vector k,
+   //
+   //    M*k = N(u+dt*k)                     (1)
+   //
+   // for nonlinear operator N. We assume N can be written as
+   //
+   //    N(u+dt*k) := L[u+dt*k](u+dt*k) + f(t)
+   //
+   // where L is a matrix-valued operator evaluated at u+dt*k and f(t)
+   // a (potentially zero) time-dependent forcing vector. (1) can be
+   // rewritten as a fixed-point equation
+   //
+   //    x = (M - dt*L[x])^{-1} (Mu + f)     (2)
+   //
+   // where x := u + dt*k, which can be solved using a Picard iteration,
+   // where a function G(x) = x is solved via iteraitons x_{k+1} = G(x_k).
+
+   double tol = 1e-6;
+   int maxiter = 100;
+
+   // Right-hand side for nonlinear iteration
+   Mmat.Mult(u, z);  // Add forcing function if one exists
+   du_dt = u;        // Set u as initial guess for x (2)
+   Vector temp(u);   // Vector to measure error
+   temp = u;
+   double error = 1;
+   int iter = 0;
+   while (error > tol) {
+      iter ++;
+      this->SetParameters(du_dt);   // Update nonlinear operator L[x]
+      T = Add(1.0, Mmat, dt, Kmat); // Form matrix (M - dt*L[x])
       T_solver.SetOperator(*T);
+      T_solver.Mult(z, du_dt);      // Apply (M - dt*L[x])^{-1}
+      temp -= du_dt;                // Measure error
+      error = std::sqrt(InnerProduct(MPI_COMM_WORLD, temp, temp));
+      temp = du_dt;
+      if (iter >= maxiter) {
+         mfem_warning("Nonlinear iteration did not converge!");
+         break;
+      }
    }
-   MFEM_VERIFY(dt == current_dt, ""); // SDIRK methods use the same dt
-   Kmat.Mult(u, z);
-   z.Neg();
-   T_solver.Mult(z, du_dt);
+
+   // Above we solved for x = u + dt*k, where k is the desired update
+   // Map du_dt -> k.
+   du_dt -= u;
+   du_dt /= dt;
 }
 
 void ConductionOperator::SetParameters(const Vector &u)
@@ -483,4 +517,4 @@ double InitialTemperature(const Vector &x)
    {
       return 1.0;
    }
-}
+}

examples/ex1p.cpp

@@ -231,10 +231,17 @@ int main(int argc, char *argv[])
    else
    {
       prec = new HypreBoomerAMG;
+      prec->SetOperator(*A);
    }
-   CGSolver cg(MPI_COMM_WORLD);
+   // CGSolver cg(MPI_COMM_WORLD);
+   AndersonAcceleration cg(MPI_COMM_WORLD);
+   cg.SetKDim(10);
+   cg.SetRestart(true);      // WORKS
+   cg.SetAAStart(0);          // WORKS
+   cg.SetWeight(1.0);         // Not robust but seems to work
+
    cg.SetRelTol(1e-12);
-   cg.SetMaxIter(2000);
+   cg.SetMaxIter(50);
    cg.SetPrintLevel(1);
    if (prec) { cg.SetPreconditioner(*prec); }
    cg.SetOperator(*A);

examples/ex1p_complex.cpp

@@ -1,435 +0,0 @@
-//                       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,10 +69,7 @@ int main(int argc, char *argv[])
    bool static_cond = false;
    bool pa = false;
    const char *device_config = "cpu";
-   bool visualization = true;
-#ifdef MFEM_USE_AMGX
-   bool useAmgX = false;
-#endif
+   bool visualization = 1;
 
    OptionsParser args(argc, argv);
    args.AddOption(&mesh_file, "-m", "--mesh",
@@ -90,11 +87,6 @@ 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())
@@ -167,10 +159,9 @@ 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())
    {
-      ess_bdr.SetSize(pmesh->bdr_attributes.Max());
+      Array<int> ess_bdr(pmesh->bdr_attributes.Max());
       ess_bdr = 1;
       fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
    }
@@ -214,20 +205,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 an AMS preconditioner.
-   if (pa)
+   // 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
    {
-#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
+      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(ams);
+      cg.SetPreconditioner(Jacobi);
       cg.Mult(B, X);
    }
    else

examples/ex3p_complex.cpp

@@ -1,517 +0,0 @@
-//                       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,41 +1770,9 @@ 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,13 +376,6 @@ 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(); }
@@ -984,18 +977,9 @@ 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,6 +1021,7 @@ void PAMixedBilinearFormExtension::Update()
       localTrial.UseDevice(true);
       localTrial.SetSize(elem_restrict_trial->Height(),
                          Device::GetMemoryType());
+
    }
    if (elem_restrict_test)
    {
@@ -1220,134 +1221,4 @@ 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,7 +21,6 @@ namespace mfem
 
 class BilinearForm;
 class MixedBilinearForm;
-class DiscreteLinearOperator;
 
 /// Class extending the BilinearForm class to support different AssemblyLevels.
 /**  FA - Full Assembly
@@ -213,7 +212,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,
@@ -259,35 +258,6 @@ 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(fes)\n"
+   mfem_error ("BilinearFormIntegrator::AssemblePA(...)\n"
                "   is not implemented for this class.");
 }
 
 void BilinearFormIntegrator::AssemblePA(const FiniteElementSpace&,
                                         const FiniteElementSpace&)
 {
-   mfem_error ("BilinearFormIntegrator::AssemblePA(fes, fes)\n"
+   mfem_error ("BilinearFormIntegrator::AssemblePA(...)\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::AddMultTransposePA(...)\n"
+   mfem_error ("BilinearFormIntegrator::MultAssembledTranspose(...)\n"
                "   is not implemented for this class.");
 }
 

fem/bilininteg.hpp

@@ -1844,10 +1844,8 @@ 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 (H1) space
-    as u.
-
-    See also MixedVectorGradientIntegrator when v is in H(curl). */
+    scalar coefficient, and v is a vector with components v_i in the same space
+    as u. */
 class GradientIntegrator : public BilinearFormIntegrator
 {
 protected:
@@ -2991,36 +2989,11 @@ 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;
 };
 
 
@@ -3035,24 +3008,6 @@ 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, "PAGradientApplyTranspose2D not implemented.");
+   MFEM_ASSERT(false, "GradientPAApplyTranspose 3D 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,7 +2638,6 @@ 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,18 +330,6 @@ 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,38 +565,6 @@ 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,9 +160,6 @@ 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(); }
@@ -304,18 +301,10 @@ 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,8 +101,6 @@ void ParFiniteElementSpace::ParInit(ParMesh *pm)
 
    P = NULL;
    Pconf = NULL;
-   Rconf = NULL;
-   R_transpose = NULL;
    R = NULL;
 
    num_face_nbr_dofs = -1;
@@ -929,45 +927,6 @@ 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; }
@@ -2881,8 +2840,6 @@ 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;
@@ -3008,12 +2965,12 @@ void ParFiniteElementSpace::Update(bool want_transform)
    }
 }
 
+
 ConformingProlongationOperator::ConformingProlongationOperator(
-   const ParFiniteElementSpace &pfes, bool local_)
+   const ParFiniteElementSpace &pfes)
    : Operator(pfes.GetVSize(), pfes.GetTrueVSize()),
      external_ldofs(),
-     gc(pfes.GroupComm()),
-     local(local_)
+     gc(pfes.GroupComm())
 {
    MFEM_VERIFY(pfes.Conforming(), "");
    const Table &group_ldof = gc.GroupLDofTable();
@@ -3062,14 +3019,7 @@ 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
-   if (local)
-   {
-      y = 0.0;
-   }
-   else
-   {
-      gc.BcastBegin(const_cast<double*>(xdata), in_layout);
-   }
+   gc.BcastBegin(const_cast<double*>(xdata), in_layout);
 
    int j = 0;
    for (int i = 0; i < m; i++)
@@ -3081,10 +3031,7 @@ 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
-   if (!local)
-   {
-      gc.BcastEnd(ydata, out_layout);
-   }
+   gc.BcastEnd(ydata, out_layout);
 }
 
 void ConformingProlongationOperator::MultTranspose(
@@ -3097,10 +3044,7 @@ void ConformingProlongationOperator::MultTranspose(
    double *ydata = y.HostWrite();
    const int m = external_ldofs.Size();
 
-   if (!local)
-   {
-      gc.ReduceBegin(xdata);
-   }
+   gc.ReduceBegin(xdata);
 
    int j = 0;
    for (int i = 0; i < m; i++)
@@ -3112,18 +3056,13 @@ 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
-   if (!local)
-   {
-      gc.ReduceEnd<double>(ydata, out_layout, GroupCommunicator::Sum);
-   }
+   gc.ReduceEnd<double>(ydata, out_layout, GroupCommunicator::Sum);
 }
 
 DeviceConformingProlongationOperator::DeviceConformingProlongationOperator(
-   const ParFiniteElementSpace &pfes,
-   bool local_) :
+   const ParFiniteElementSpace &pfes) :
    ConformingProlongationOperator(pfes),
-   mpi_gpu_aware(Device::GetGPUAwareMPI()),
-   local(local_)
+   mpi_gpu_aware(Device::GetGPUAwareMPI())
 {
    MFEM_ASSERT(pfes.Conforming(), "internal error");
    const SparseMatrix *R = pfes.GetRestrictionMatrix();
@@ -3240,42 +3179,32 @@ void DeviceConformingProlongationOperator::Mult(const Vector &x,
                                                 Vector &y) const
 {
    const GroupTopology &gtopo = gc.GetGroupTopology();
+   BcastBeginCopy(x); // copy to 'shr_buf'
    int req_counter = 0;
-   if (local)
+   for (int nbr = 1; nbr < gtopo.GetNumNeighbors(); nbr++)
    {
-      y = 0.0;
-   }
-   else
-   {
-      BcastBeginCopy(x); // copy to 'shr_buf'
-      for (int nbr = 1; nbr < gtopo.GetNumNeighbors(); nbr++)
+      const int send_offset = shr_buf_offsets[nbr];
+      const int send_size = shr_buf_offsets[nbr+1] - send_offset;
+      if (send_size > 0)
       {
-         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++]);
-         }
+         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);
-   if (!local)
-   {
-      MPI_Waitall(req_counter, requests, MPI_STATUSES_IGNORE);
-      BcastEndCopy(y); // copy from 'ext_buf'
-   }
+   MPI_Waitall(req_counter, requests, MPI_STATUSES_IGNORE);
+   BcastEndCopy(y); // copy from 'ext_buf'
 }
 
 DeviceConformingProlongationOperator::~DeviceConformingProlongationOperator()
@@ -3338,38 +3267,32 @@ void DeviceConformingProlongationOperator::MultTranspose(const Vector &x,
                                                          Vector &y) const
 {
    const GroupTopology &gtopo = gc.GetGroupTopology();
+   ReduceBeginCopy(x); // copy to 'ext_buf'
    int req_counter = 0;
-   if (!local)
+   for (int nbr = 1; nbr < gtopo.GetNumNeighbors(); nbr++)
    {
-      ReduceBeginCopy(x); // copy to 'ext_buf'
-      for (int nbr = 1; nbr < gtopo.GetNumNeighbors(); nbr++)
+      const int send_offset = ext_buf_offsets[nbr];
+      const int send_size = ext_buf_offsets[nbr+1] - send_offset;
+      if (send_size > 0)
       {
-         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++]);
-         }
+         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);
-   if (!local)
-   {
-      MPI_Waitall(req_counter, requests, MPI_STATUSES_IGNORE);
-      ReduceEndAssemble(y); // assemble from 'shr_buf'
-   }
+   MPI_Waitall(req_counter, requests, MPI_STATUSES_IGNORE);
+   ReduceEndAssemble(y); // assemble from 'shr_buf'
 }
 
 } // namespace mfem

fem/pfespace.hpp

@@ -75,12 +75,6 @@ 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); }
@@ -347,16 +341,6 @@ 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; }
@@ -411,11 +395,9 @@ class ConformingProlongationOperator : public Operator
 protected:
    Array<int> external_ldofs;
    const GroupCommunicator &gc;
-   bool local;
 
 public:
-   ConformingProlongationOperator(const ParFiniteElementSpace &pfes,
-                                  bool local_=false);
+   ConformingProlongationOperator(const ParFiniteElementSpace &pfes);
 
    virtual void Mult(const Vector &x, Vector &y) const;
 
@@ -434,8 +416,6 @@ 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;
@@ -461,8 +441,7 @@ protected:
    void ReduceEndAssemble(Vector &dst) const;
 
 public:
-   DeviceConformingProlongationOperator(const ParFiniteElementSpace &pfes,
-                                        bool local_=false);
+   DeviceConformingProlongationOperator(const ParFiniteElementSpace &pfes);
 
    virtual ~DeviceConformingProlongationOperator();
 

fem/restriction.cpp

@@ -195,31 +195,6 @@ 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,10 +57,6 @@ 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,7 +10,6 @@
 # CONTRIBUTING.md for details.
 
 list(APPEND SRCS
-  auxiliary.cpp
   blockmatrix.cpp
   blockoperator.cpp
   blockvector.cpp
@@ -28,7 +27,6 @@ list(APPEND SRCS
   )
 
 list(APPEND HDRS
-  auxiliary.hpp
   blockmatrix.hpp
   blockoperator.hpp
   blockvector.hpp

linalg/auxiliary.hpp

@@ -1,285 +0,0 @@
-// 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,7 +29,6 @@
 #include "solvers.hpp"
 #include "handle.hpp"
 #include "invariants.hpp"
-#include "auxiliary.hpp"
 
 #ifdef MFEM_USE_AMGX
 #include "amgxsolver.hpp"

linalg/operator.hpp

@@ -37,8 +37,7 @@ protected:
       const Array<int> &test_tdof_list,
       RectangularConstrainedOperator* &Aout);
 
-   /** @brief Returns RAP Operator of this, using input/output Prolongation matrices
-       @a Pi corresponds to "P", @a Po corresponds to "Rt" */
+   /// Returns RAP Operator of this, taking in input/output Prolongation matrices
    Operator *SetupRAP(const Operator *Pi, const Operator *Po);
 
 public:
@@ -113,11 +112,6 @@ 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
@@ -657,7 +651,6 @@ public:
    virtual void SetOperator(const Operator &op) = 0;
 };
 
-
 /// Identity Operator I: x -> x.
 class IdentityOperator : public Operator
 {
@@ -696,42 +689,6 @@ 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.cpp

@@ -1733,6 +1733,323 @@ void NewtonSolver::AdaptiveLinRtolPostSolve(const Vector &x,
    }
 }
 
+void AndersonAcceleration::QRdelete(std::deque<Vector *> &Q, DenseMatrix &R) const
+{
+   Vector temp(Q[0]->Size());
+   for (int i=0; i<(maxVecs-1); i++) {
+      double d = sqrt( R(i,i+1)*R(i,i+1) + R(i+1,i+1)*R(i+1,i+1) );
+      double c = R(i,i+1) / d;
+      double s = R(i+1,i+1) / d;
+      R(i,i+1) = d;
+      R(i+1,i+1) = 0;
+
+      if (i < (maxVecs-2)) {
+         for (int j=(i+2); j<maxVecs; j++) {
+            d = c*R(i,j) + s*R(i+1,j);
+            R(i+1,j) = -s*R(i,j) + c*R(i+1,j);
+            R(i,j) = d;
+         }
+      }
+      // temp = c*Q[i] + s*Q[i+1];
+      add(c, *(Q[i]), s, *(Q[i+1]), temp);
+      // Q[i+1] = -s*Q[i] + c*Q[i+1];
+      *(Q[i+1]) *= c;
+      Q[i+1] -> Add(-s, *(Q[i]));
+      *(Q[i]) = temp;
+   }
+   
+   // Shift Q <- Q[:,0:(m-2)], i.e., delete last column of Q
+   delete Q.back();
+   Q.pop_back();
+
+   // Shift columns of R to the left by one, R = R[0:(m−2), 1:(m-1)]
+   for (int j=1; j<maxVecs; j++) {
+      for (int i=0; i<maxVecs; i++) {
+         R(i,j-1) = R(i,j);
+      }
+   }
+}
+
+void AndersonAcceleration::SetOperator(const Operator &op)
+{
+   oper = &op;
+   height = op.Height();
+   width = op.Width();
+   MFEM_ASSERT(height == width, "square Operator is required.");
+}
+
+void AndersonAcceleration::FixedPointMult(const Vector &b,
+   const Vector &x, Vector &y, Vector &r) const
+{
+   // Assume Operator represents a fixed-point operator *with
+   // right-hand side*, so we iterate x_{k+1} = M^{-1}G(x_k)
+   if (isFixedPointOp) {
+      oper->Mult(x, y);
+      if (prec) {
+         prec->Mult(y,r);
+         y = r;
+         r -= x;
+      }
+      else {
+         r = y;
+         r -= x;
+      }
+   }
+   // Otherwise add x to y and not to residual r = y - x
+   else {
+      oper->Mult(x, r);
+      r *= -1;
+      r += b;
+      if (prec) {
+         prec->Mult(r,y);
+         r = y;
+         y += x;
+      }
+      else {
+         y = r;
+         y += x;         
+      }
+   }
+}
+
+void AndersonAcceleration::Mult(const Vector &b, Vector &x) const
+{
+   MFEM_ASSERT(oper != NULL, "the Operator is not set (use SetOperator).");
+  
+   int n = width;
+   int numVecs = 0;
+   double resid, norm_df;
+   double min_diag = 1e-13;
+   final_norm = -1;
+
+   // Check vector is initialized, set to zero for
+   // iterative_mode = false
+   if (x.Size() != n) {
+      x.SetSize(n);
+      x = 0.0;
+   }
+   else if (!iterative_mode) {
+      x = 0.0;
+   }
+
+   // Storage containers for acceleration
+   std::deque<Vector *> G;
+   std::deque<Vector *> Q;
+   DenseMatrix R(maxVecs);
+   R = 0.0;
+   Vector g_old(n);
+   Vector g_current(n);
+   Vector f_old(n);
+   Vector f_current(n);
+   Vector gamma(maxVecs);
+   Vector rhs(maxVecs);
+   Vector correction;
+   if (omega > 0 && std::abs(omega -  1) > 1e-14) {
+      correction.SetSize(n);
+   }
+   Vector *dg;
+   Vector *df;
+
+   // Loop over AA iterations
+   int k;
+   for (k=0; k<max_iter; k++) {
+
+      // Compute g_current = G(x), f_current = G(x) - x
+      this->FixedPointMult(b, x, g_current, f_current); 
+
+      // Check norm of current approximation to fixed point G(u) = u
+      resid = Norm(f_current);
+      MFEM_ASSERT(IsFinite(resid), "||G(u) - u|| = " << resid);
+      if (print_level == 1)
+      {
+         mfem::out << "  Iteration : " << setw(3) << k
+                 << "  ||G(u) - u|| = " << resid << endl;
+      }
+      
+      // Set stopping tolerance on first iteration.
+      if (final_norm < 0) {
+         final_norm = std::max(rel_tol*resid, abs_tol);
+      }
+
+      // Check for convergence
+      if (resid <= final_norm)
+      {
+         final_norm = resid;
+         final_iter = k;
+         converged = 1;
+         goto finish;
+      }
+
+      // Start Anderson Acceleration after AAstart FP iterations
+      if (k > AAstart) {
+         // df = f_current - f_old; 
+         df = new Vector(n);
+         add(1.0, f_current, -1.0, f_old, *df);
+         // dg = g_current - g_old;
+         dg = new Vector(n);
+         add(1.0, g_current, -1.0, g_old, *dg);
+      
+         if (numVecs < maxVecs) {
+            G.push_back(dg);
+         }
+         else {
+            delete G[0];
+            G.pop_front();
+            G.push_back(dg);
+         }
+         numVecs++;
+         dg = NULL;
+      }
+      
+      f_old = f_current;
+      g_old = g_current;
+      
+      // First iteration or initial fixed-point iterations
+      if (numVecs == 0) {
+         x = g_current;
+         continue;
+      }
+
+      // All later iterations: orthogonalize and find best approximation
+      if (numVecs == 1) {
+         norm_df = Norm(*df);
+         MFEM_ASSERT(IsFinite(norm_df), "norm_df = " << norm_df);
+         (*df) /= norm_df;
+         Q.push_back(df);
+         R(0,0) = norm_df;
+         df = NULL;
+      }
+      else {
+         // Remove first column in basis F and R, reorthogonalize
+         if (numVecs > maxVecs) {
+            this->QRdelete(Q, R);
+            numVecs--;
+         }
+         // Compute last column of R
+         for (int i=0; i<(numVecs-1); i++) {
+            R(i,numVecs-1) = Dot(*(Q[i]), *df);
+            // df -= R(i,numVecs-1) * Q[i]
+            df -> Add(-R(i,numVecs-1), *(Q[i]));
+         }
+         norm_df = Norm(*df);
+         MFEM_ASSERT(IsFinite(norm_df), "norm_df = " << norm_df);
+         (*df) /= norm_df;
+         Q.push_back(df);
+         R(numVecs-1, numVecs-1) = norm_df;
+         df = NULL;
+      }
+
+      // Back solve for new weights, R\gamma = Q^T * f_current
+      rhs = 0.0;
+      gamma = 0.0;
+      for (int i=0; i<numVecs; i++) {
+         rhs(i) = Dot(*(Q[i]), f_current);   // Form right hand side
+      }
+      for (int i=(numVecs-1); i>=0; i--) {
+         double temp = rhs(i);
+         for (int j=(i+1); j<numVecs; j++) {
+            temp -= R(i,j)*gamma(j);
+         }
+         if (std::abs(R(i,i)) < min_diag) {
+            gamma(i) = 0.0;
+            std::cout << "Diagonal of R -- " << R(i,i) << " ~ 0.\n";
+         }
+         else {
+            gamma(i) = temp / R(i,i);            
+         }
+      }
+
+      /// DEBUG --> test backsolve
+      Vector test(numVecs);
+      for (int i=0; i<numVecs; i++) {
+         test(i) = 0;
+         for (int j=i; j<numVecs; j++) {
+            test(i) += R(i,j) * gamma(j);
+         }
+         if (std::abs(rhs(i) - test(i)) > 1e-10) {
+            std::cout << "Bad solve! Err = " << rhs(i) - test(i) << "\n";
+         }
+      }
+
+      // Compute updated solution x = g_current − G*\gamma 
+      x = g_current;
+      for (int i=0; i<numVecs; i++) {
+         // x -= gamma(i)*G[i]
+         x.Add(-gamma(i), *(G[i]));
+      }
+
+      // Apply damped iteration for \omega \in (0,1),
+      //   x -= (1−omega) * (f_current − Q*R*gamma);
+      if (omega > 0 && std::abs(omega -  1) > 1e-14) {
+         // Redefine rhs = R*gamma
+         for(int i=0; i<numVecs; i++) {
+            rhs(i) = 0;
+            for (int j=i; j<numVecs; j++) {
+               rhs(i) += R(i,j)*gamma(j);
+            }
+         }
+         correction = f_current;
+         for (int i=0; i<numVecs; i++) {
+            // correction -= rhs(i)*Q[i]
+            correction.Add(-rhs(i), *(Q[i]));
+         }
+         // x -= (1 - omega) * correction;
+         x.Add( -(1 - omega), correction);
+      }
+
+      // Restart AA minimization by eliminating all vectors but the most recent
+      if (restart && (numVecs == maxVecs)) {
+         for (int i=0; i<(maxVecs-1); i++) {
+            delete G[0];
+            G.pop_front();          
+            delete Q[0];
+            Q.pop_front();
+         }
+         R = 0.0;
+         R(0,0) = norm_df;
+         numVecs = 1;
+         if (print_level == 1)
+         {
+            mfem::out << "Restarting..." << '\n';
+         }
+      }
+   }
+
+   // Compute final residual, save counts for solve
+   this->FixedPointMult(b, x, g_current, f_current); 
+   resid = Norm(f_current);
+   MFEM_ASSERT(IsFinite(resid), "||G(u) - u|| = " << resid);
+   final_norm = resid;
+   final_iter = max_iter;
+   if (resid <= final_norm) converged = 1;
+   else converged = 0;
+
+finish:
+   if (print_level == 3)
+   {
+      mfem::out << "  Iteration : " << setw(3) << k
+              << "  ||G(u) - u|| = " << resid << endl;
+   }
+   else if (print_level == 2)
+   {
+      mfem::out << "Anderson Acceleration: Number of iterations: " << final_iter << '\n';
+   }
+   if (print_level >= 0 && !converged)
+   {
+      mfem::out << "Anderson Acceleration: No convergence!\n";
+   }
+
+   // Cleanup pointers
+   for (int i=0; i<numVecs; i++) {
+      delete G[0];
+      G.pop_front();          
+      delete Q[0];
+      Q.pop_front();
+   }
+   delete dg;
+   delete df;
+}
+
 void LBFGSSolver::Mult(const Vector &b, Vector &x) const
 {
    MFEM_VERIFY(oper != NULL, "the Operator is not set (use SetOperator).");

linalg/solvers.hpp

@@ -14,6 +14,7 @@
 
 #include "../config/config.hpp"
 #include "densemat.hpp"
+#include <deque>
 
 #ifdef MFEM_USE_MPI
 #include <mpi.h>
@@ -489,6 +490,64 @@ public:
                            const double gamma = 1.0);
 };
 
+
+/// Nonlinear Anderson Acceleration
+class AndersonAcceleration : public IterativeSolver
+{
+protected:
+   int maxVecs; // see SetKDim()
+   int AAstart;
+   bool restart;
+   bool isFixedPointOp;
+   double omega;
+
+   /// Apply fixed-point Mult and compute residual.
+   //  For !isFixedPointOp:
+   //    y <-- G(x) = x + A(x) - b and r <-- G(x) - x = A(x) - b
+   //  For isFixedPointOp:
+   //    y <-- G(x) = A(x) - b and r <-- G(x) - x = A(x) - b - x
+   void FixedPointMult(const Vector &b, const Vector &x, Vector &y, Vector &r) const;
+
+   // Helper function for Anderson Acceleration
+   void QRdelete(std::deque<Vector *> &Q, DenseMatrix &R) const;
+
+public:
+   AndersonAcceleration() : maxVecs(25), AAstart(0), omega(1),
+      restart(false), isFixedPointOp(false) { }
+
+#ifdef MFEM_USE_MPI
+   AndersonAcceleration(MPI_Comm _comm) : IterativeSolver(_comm),
+      maxVecs(25), AAstart(0), omega(1), restart(false),
+      isFixedPointOp(false) { }
+#endif
+
+   /// Boolean describing whether the action of the operator is such
+   //  that we want to solve G(x) = x (true) or G(x) = 0 (false) for
+   //  zero right-hand side vector b passed into Mult(). Default is
+   //  false in construction of class.
+   void IsFixedPointOperator(bool isFixedPointOp_)
+   { isFixedPointOp = isFixedPointOp_; }
+
+   /// Maximum number of vectors to store in Krylov-like space
+   void SetKDim(int dim) { maxVecs = dim; }
+
+   /// Number of fixed-point iterations to do before starting AA
+   void SetAAStart(int start_) { AAstart = start_; }
+
+   /// Boolean to restart, that is, erase entire space after maxVecs
+   //  are stored (AAstart=true) or use a sliding space (AAstart=false)
+   //  where one vector is deleted to make room for a new one.
+   void SetRestart(bool restart_) { restart = restart_; }
+
+   /// Set relaxation weight
+   void SetWeight(double omega_) { omega = omega_; }
+
+   virtual void SetOperator(const Operator &op);
+
+   virtual void Mult(const Vector &b, Vector &x) const;
+};
+
+
 /** L-BFGS method for solving F(x)=b for a given operator F, by minimizing
     the norm of F(x) - b. Requires only the action of the operator F. */
 class LBFGSSolver : public NewtonSolver
@@ -852,169 +911,6 @@ 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,8 +55,6 @@ 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

@@ -1,219 +0,0 @@
-// 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

@@ -1,123 +0,0 @@
-// 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));
-}