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

examples/ex1p_complex.cpp

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

examples/ex3p_complex.cpp

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

linalg/auxiliary.cpp

@@ -16,11 +16,13 @@
 #include "linalg.hpp"
 #include "../fem/pfespace.hpp"
 #include "../fem/pbilinearform.hpp"
+#include "../fem/complex_fem.hpp"
 
 namespace mfem
 {
 
 GeneralAMS::GeneralAMS(const Operator& curlcurl_op_,
+                       Operator *oper_complex_,
                        const Operator& pi_,
                        const Operator& gradient_,
                        const Operator& pispacesolver_,
@@ -30,6 +32,7 @@ GeneralAMS::GeneralAMS(const Operator& curlcurl_op_,
    :
    Solver(curlcurl_op_.Height()),
    curlcurl_op(curlcurl_op_),
+   oper_complex(oper_complex_),
    pi(pi_),
    gradient(gradient_),
    pispacesolver(pispacesolver_),
@@ -46,7 +49,15 @@ GeneralAMS::~GeneralAMS()
 void GeneralAMS::FormResidual(const Vector& rhs, const Vector& x,
                               Vector& residual) const
 {
-   curlcurl_op.Mult(x, residual);
+   if (oper_complex)
+   {
+      oper_complex->Mult(x, residual);
+   }
+   else
+   {
+      curlcurl_op.Mult(x, residual);
+   }
+
    residual *= -1.0;
    residual += rhs;
 }
@@ -70,7 +81,8 @@ void GeneralAMS::FormResidual(const Vector& rhs, const Vector& x,
 void GeneralAMS::Mult(const Vector& x, Vector& y) const
 {
    MFEM_ASSERT(x.Size() == y.Size(), "Sizes don't match!");
-   MFEM_ASSERT(curlcurl_op.Height() == x.Size(), "Sizes don't match!");
+   MFEM_ASSERT(x.Size() == ((oper_complex != NULL) ? oper_complex->Height() :
+                            curlcurl_op.Height()), "Sizes don't match!");
 
    Vector residual(x.Size());
    residual = 0.0;
@@ -118,8 +130,8 @@ void GeneralAMS::Mult(const Vector& x, Vector& y) const
 
 // Pi-space constructor
 MatrixFreeAuxiliarySpace::MatrixFreeAuxiliarySpace(
-   ParMesh& mesh_lor, Coefficient* alpha_coeff,
-   Coefficient* beta_coeff, MatrixCoefficient* beta_mcoeff, Array<int>& ess_bdr,
+   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_,
@@ -132,7 +144,8 @@ MatrixFreeAuxiliarySpace::MatrixFreeAuxiliarySpace(
 #ifdef MFEM_USE_AMGX
    useAmgX(useAmgX_),
 #endif
-   inner_aux_iterations(0)
+   inner_aux_iterations(0),
+   imagBdry(false)
 {
    H1_FECollection * fec_lor = new H1_FECollection(1, mesh_lor.Dimension());
    ParFiniteElementSpace fespace_lor_d(&mesh_lor, fec_lor, mesh_lor.Dimension(),
@@ -143,6 +156,7 @@ MatrixFreeAuxiliarySpace::MatrixFreeAuxiliarySpace(
    {
       fespace_lor_d.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
    }
+
    ParBilinearForm a_space(&fespace_lor_d);
 
    // this choice of policy is important for the G-space solver, but
@@ -176,6 +190,7 @@ MatrixFreeAuxiliarySpace::MatrixFreeAuxiliarySpace(
    a_space.EliminateEssentialBC(ess_bdr, policy);
    a_space.Finalize();
    aspacematrix = a_space.ParallelAssemble();
+
    aspacematrix->CopyRowStarts();
    aspacematrix->CopyColStarts();
 
@@ -192,6 +207,205 @@ MatrixFreeAuxiliarySpace::MatrixFreeAuxiliarySpace(
    delete fec_lor;
 }
 
+// Complex Pi-space constructor
+MatrixFreeAuxiliarySpace::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) :
+   Solver(2*pi.Width()),
+   comm(mesh_lor.GetComm()),
+   matfree(NULL),
+   cg(NULL),
+#ifdef MFEM_USE_AMGX
+   useAmgX(useAmgX_),
+#endif
+   inner_aux_iterations(0),
+   imagBdry(beta_imag != NULL && abs_beta_imag != NULL)
+{
+   MFEM_VERIFY(imagBdry, "");
+   MFEM_VERIFY(2*curlcurl_oper.Height() == oper_complex->Height(), "");
+
+   H1_FECollection * fec_lor = new H1_FECollection(1, mesh_lor.Dimension());
+   ParFiniteElementSpace fespace_lor_d(&mesh_lor, fec_lor, mesh_lor.Dimension(),
+                                       Ordering::byVDIM);
+
+   offsets.SetSize(3);
+   offsets[0] = 0;
+   offsets[1] = fespace_lor_d.GetTrueVSize();
+   offsets[2] = offsets[1];
+   offsets.PartialSum();
+
+   offsets_nd.SetSize(3);
+   offsets_nd[0] = 0;
+   offsets_nd[1] = curlcurl_oper.Height();
+   offsets_nd[2] = offsets_nd[1];
+   offsets_nd.PartialSum();
+
+   // build LOR AMG v-cycle
+   if (ess_bdr.Size())
+   {
+      fespace_lor_d.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+
+   {
+      // Assemble real system
+      ParBilinearForm a_space(&fespace_lor_d);
+
+      // this choice of policy is important for the G-space solver, but
+      // also can make some difference here
+      const Matrix::DiagonalPolicy policy = Matrix::DIAG_KEEP;
+      a_space.SetDiagonalPolicy(policy);
+      if (alpha_coeff == NULL)
+      {
+         a_space.AddDomainIntegrator(new VectorDiffusionIntegrator);
+      }
+      else
+      {
+         a_space.AddDomainIntegrator(new VectorDiffusionIntegrator(*alpha_coeff));
+      }
+
+      if (beta_mcoeff != NULL)
+      {
+         MFEM_VERIFY(beta_coeff == NULL, "Only one beta coefficient should be defined.");
+         a_space.AddDomainIntegrator(new VectorMassIntegrator(*beta_mcoeff));
+      }
+      else if (beta_coeff != NULL)
+      {
+         a_space.AddDomainIntegrator(new VectorMassIntegrator(*beta_coeff));
+      }
+      else
+      {
+         a_space.AddDomainIntegrator(new VectorMassIntegrator);
+      }
+      a_space.UsePrecomputedSparsity();
+      a_space.Assemble();
+      a_space.EliminateEssentialBC(ess_bdr, policy);
+      a_space.Finalize();
+      aspacematrix = a_space.ParallelAssemble();
+
+      aspacematrix->CopyRowStarts();
+      aspacematrix->CopyColStarts();
+   }
+
+   {
+      // Assemble complex system
+      ParSesquilinearForm a_space(&fespace_lor_d);
+
+      // this choice of policy is important for the G-space solver, but
+      // also can make some difference here
+      // NOTE: without essential BC, the policy is not applicable anyway
+      //const Matrix::DiagonalPolicy policy = Matrix::DIAG_KEEP;
+      //a_space.SetDiagonalPolicy(policy);
+      if (alpha_coeff == NULL)
+      {
+         a_space.AddDomainIntegrator(new VectorDiffusionIntegrator, NULL);
+      }
+      else
+      {
+         a_space.AddDomainIntegrator(new VectorDiffusionIntegrator(*alpha_coeff), NULL);
+      }
+
+      if (beta_mcoeff != NULL)
+      {
+         MFEM_VERIFY(beta_coeff == NULL, "Only one beta coefficient should be defined.");
+         a_space.AddDomainIntegrator(new VectorMassIntegrator(*beta_mcoeff), NULL);
+      }
+      else if (beta_coeff != NULL)
+      {
+         a_space.AddDomainIntegrator(new VectorMassIntegrator(*beta_coeff), NULL);
+      }
+      else
+      {
+         a_space.AddDomainIntegrator(new VectorMassIntegrator, NULL);
+      }
+
+      a_space.AddBoundaryIntegrator(NULL,
+                                    new VectorMassIntegrator(*beta_imag));  // im part
+
+      //a_space.UsePrecomputedSparsity();
+      a_space.Assemble();
+      //a_space.EliminateEssentialBC(ess_bdr, policy);
+      a_space.Finalize();
+
+      //aspacematrix = a_space.ParallelAssemble();
+      aspacematrix_complex = a_space.ParallelAssemble()->GetSystemMatrix();
+      /*
+      {
+         Array<int> empty_ess_tdof_list;
+         OperatorPtr Aptr;
+         a_space.FormSystemMatrix(empty_ess_tdof_list, Aptr);
+
+         aspacematrix_complex = Aptr.As<HypreParMatrix>();
+      }
+      */
+      aspacematrix_complex->CopyRowStarts();
+      aspacematrix_complex->CopyColStarts();
+   }
+
+   {
+      // Assemble imaginary system
+      ParBilinearForm a_space(&fespace_lor_d);
+
+      a_space.AddBoundaryIntegrator(new VectorMassIntegrator(
+                                       *abs_beta_imag));  // im part
+
+      //a_space.UsePrecomputedSparsity();
+      a_space.Assemble();
+      //a_space.EliminateEssentialBC(ess_bdr, policy);
+      a_space.Finalize();
+
+      aspacematrix_imag = a_space.ParallelAssemble();
+      /*
+      {
+         Array<int> empty_ess_tdof_list;
+         OperatorPtr Aptr;
+         a_space.FormSystemMatrix(empty_ess_tdof_list, Aptr);
+
+         aspacematrix_imag = Aptr.As<HypreParMatrix>();
+      }
+      */
+
+      aspacematrix_imag->CopyRowStarts();
+      aspacematrix_imag->CopyColStarts();
+   }
+
+   //SetupAMG(fespace_lor_d.GetMesh()->Dimension());
+   {
+      HypreBoomerAMG *amg = new HypreBoomerAMG(*aspacematrix);
+      const int system_dimension = fespace_lor_d.GetMesh()->Dimension();
+      amg->SetSystemsOptions(system_dimension);
+      amg->SetPrintLevel(0);
+      aspacepc = amg;
+   }
+
+   /*
+
+   if (cg_iterations > 0)
+   {
+   SetupCG(curlcurl_oper, pi, cg_iterations);
+   }
+   else
+   {
+   SetupVCycle();
+   }
+   */
+
+   SetupPMHSS();
+
+   conn_block = new BlockOperator(offsets_nd, offsets);
+   conn_block->SetDiagonalBlock(0, &pi);
+   conn_block->SetDiagonalBlock(1, &pi);
+
+   SetupGMRES(*oper_complex, *conn_block);
+
+   delete fec_lor;
+}
+
 /* G-space constructor
 
    The auxiliary space solves in general, and this one in particular,
@@ -216,7 +430,8 @@ MatrixFreeAuxiliarySpace::MatrixFreeAuxiliarySpace(
 #ifdef MFEM_USE_AMGX
    useAmgX(useAmgX_),
 #endif
-   inner_aux_iterations(0)
+   inner_aux_iterations(0),
+   imagBdry(false)
 {
    H1_FECollection * fec_lor = new H1_FECollection(1, mesh_lor.Dimension());
    ParFiniteElementSpace fespace_lor(&mesh_lor, fec_lor);
@@ -275,6 +490,218 @@ MatrixFreeAuxiliarySpace::MatrixFreeAuxiliarySpace(
    delete fec_lor;
 }
 
+/* Complex G-space constructor
+
+   The auxiliary space solves in general, and this one in particular,
+   seem to be quite sensitive to handling of boundary conditions. Note
+   some careful choices for Matrix::DiagonalPolicy and the ZeroWrap
+   object, as well as the use of a single CG iteration (instead of just
+   an AMG V-cycle). Just a V-cycle may be more efficient in some cases,
+   but we recommend the CG wrapper for robustness here. */
+MatrixFreeAuxiliarySpace::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)
+   :
+   Solver(curlcurl_oper.Height()),
+   comm(mesh_lor.GetComm()),
+   matfree(NULL),
+   cg(NULL),
+#ifdef MFEM_USE_AMGX
+   useAmgX(useAmgX_),
+#endif
+   inner_aux_iterations(0),
+   imagBdry(beta_imag != NULL && abs_beta_imag != NULL)
+{
+   MFEM_VERIFY(imagBdry, "");
+   MFEM_VERIFY(2*curlcurl_oper.Height() == oper_complex->Height(), "");
+
+   H1_FECollection * fec_lor = new H1_FECollection(1, mesh_lor.Dimension());
+   ParFiniteElementSpace fespace_lor(&mesh_lor, fec_lor);
+
+   offsets.SetSize(3);
+   offsets[0] = 0;
+   offsets[1] = fespace_lor.GetTrueVSize();
+   offsets[2] = offsets[1];
+   offsets.PartialSum();
+
+   offsets_nd.SetSize(3);
+   offsets_nd[0] = 0;
+   offsets_nd[1] = curlcurl_oper.Height();
+   offsets_nd[2] = offsets_nd[1];
+   offsets_nd.PartialSum();
+
+   const double regeps = 1.0e-6;
+   ConstantCoefficient epscoef(regeps);
+
+   {
+      // Assemble real system
+
+      // build LOR AMG v-cycle
+      ParBilinearForm a_space(&fespace_lor);
+
+      // we need something like DIAG_ZERO in the solver, but explicitly doing
+      // that makes BoomerAMG setup complain, so instead we constrain the boundary
+      // in the CG solver
+      const Matrix::DiagonalPolicy policy = Matrix::DIAG_ONE;
+
+      a_space.SetDiagonalPolicy(policy);
+
+      if (beta_mcoeff != NULL)
+      {
+         MFEM_VERIFY(beta_coeff == NULL, "Only one beta coefficient should be defined.");
+         a_space.AddDomainIntegrator(new DiffusionIntegrator(*beta_mcoeff));
+      }
+      else if (beta_coeff != NULL)
+      {
+         a_space.AddDomainIntegrator(new DiffusionIntegrator(*beta_coeff));
+      }
+      else
+      {
+         a_space.AddDomainIntegrator(new DiffusionIntegrator);
+      }
+      a_space.AddDomainIntegrator(new MassIntegrator(epscoef));
+
+      a_space.UsePrecomputedSparsity();
+      a_space.Assemble();
+      if (ess_bdr.Size())
+      {
+         fespace_lor.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+      }
+
+      // you have to use (serial) BilinearForm eliminate routines to get
+      // diag policy DIAG_ZERO all the ParallelEliminateTDofs etc. routines
+      // implicitly have a Matrix::DIAG_KEEP policy
+      a_space.EliminateEssentialBC(ess_bdr, policy);
+      a_space.Finalize();
+      aspacematrix = a_space.ParallelAssemble();
+
+      aspacematrix->CopyRowStarts();
+      aspacematrix->CopyColStarts();
+   }
+
+   {
+      // Assemble complex system
+
+      // build LOR AMG v-cycle
+      ParSesquilinearForm a_space(&fespace_lor);
+
+      // we need something like DIAG_ZERO in the solver, but explicitly doing
+      // that makes BoomerAMG setup complain, so instead we constrain the boundary
+      // in the CG solver
+      const Matrix::DiagonalPolicy policy = Matrix::DIAG_ONE;
+
+      //a_space.SetDiagonalPolicy(policy);
+
+      if (beta_mcoeff != NULL)
+      {
+         MFEM_VERIFY(beta_coeff == NULL, "Only one beta coefficient should be defined.");
+         a_space.AddDomainIntegrator(new DiffusionIntegrator(*beta_mcoeff), NULL);
+      }
+      else if (beta_coeff != NULL)
+      {
+         a_space.AddDomainIntegrator(new DiffusionIntegrator(*beta_coeff), NULL);
+      }
+      else
+      {
+         a_space.AddDomainIntegrator(new DiffusionIntegrator, NULL);
+      }
+      a_space.AddDomainIntegrator(new MassIntegrator(epscoef), NULL);
+
+      a_space.AddBoundaryIntegrator(NULL, new MassIntegrator(*beta_imag));  // im part
+
+      //a_space.UsePrecomputedSparsity();
+      a_space.Assemble();
+      /*
+      if (ess_bdr.Size())
+      {
+         fespace_lor.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+      }
+      */
+
+      // you have to use (serial) BilinearForm eliminate routines to get
+      // diag policy DIAG_ZERO all the ParallelEliminateTDofs etc. routines
+      // implicitly have a Matrix::DIAG_KEEP policy
+      //a_space.EliminateEssentialBC(ess_bdr, policy);
+      a_space.Finalize();
+      //ComplexHypreParMatrix *complex_aspacematrix = a_space.ParallelAssemble();
+      aspacematrix_complex = a_space.ParallelAssemble()->GetSystemMatrix();
+
+      /*
+      {
+         Array<int> empty_ess_tdof_list;
+         OperatorPtr Aptr;
+         a_space.FormSystemMatrix(empty_ess_tdof_list, Aptr);
+
+         aspacematrix_complex = Aptr.As<HypreParMatrix>();
+      }
+      */
+
+      aspacematrix_complex->CopyRowStarts();
+      aspacematrix_complex->CopyColStarts();
+      // TODO: aspacematrix_complex is not used?
+   }
+
+   {
+      // Assemble imaginary system
+      ParBilinearForm a_space(&fespace_lor);
+
+      a_space.AddBoundaryIntegrator(new MassIntegrator(*abs_beta_imag));  // im part
+
+      //a_space.UsePrecomputedSparsity();
+      a_space.Assemble();
+      //a_space.EliminateEssentialBC(ess_bdr, policy);
+      a_space.Finalize();
+
+      aspacematrix_imag = a_space.ParallelAssemble();
+      /*
+      {
+         Array<int> empty_ess_tdof_list;
+         OperatorPtr Aptr;
+         a_space.FormSystemMatrix(empty_ess_tdof_list, Aptr);
+
+         aspacematrix_imag = Aptr.As<HypreParMatrix>();
+      }
+      */
+
+      aspacematrix_imag->CopyRowStarts();
+      aspacematrix_imag->CopyColStarts();
+   }
+
+   //SetupAMG(0);
+   {
+      HypreBoomerAMG *amg = new HypreBoomerAMG(*aspacematrix);
+      amg->SetPrintLevel(0);
+      aspacepc = amg;
+   }
+
+   /*
+   if (cg_iterations > 0)
+   {
+      SetupCG(curlcurl_oper, g, cg_iterations);
+   }
+   else
+   {
+      SetupVCycle();
+   }
+   */
+
+   SetupPMHSS();
+
+   conn_block = new BlockOperator(offsets_nd, offsets);
+   conn_block->SetDiagonalBlock(0, &g);
+   conn_block->SetDiagonalBlock(1, &g);
+
+   SetupGMRES(*oper_complex, *conn_block);
+
+   delete fec_lor;
+}
+
 void MatrixFreeAuxiliarySpace::SetupCG(
    Operator& curlcurl_oper, Operator& conn,
    int inner_cg_iterations)
@@ -303,11 +730,49 @@ void MatrixFreeAuxiliarySpace::SetupCG(
    aspacewrapper = cg;
 }
 
+void MatrixFreeAuxiliarySpace::SetupGMRES(Operator& curlcurl_oper,
+                                          Operator& conn)
+{
+   MFEM_ASSERT(conn.Height() == curlcurl_oper.Width(),
+               "Operators don't match!");
+   matfree = new RAPOperator(conn, curlcurl_oper, conn);
+   MFEM_ASSERT(matfree->Height() == PMHSS->Height(),
+               "Operators don't match!");
+
+   {
+      gmres_PMHSS = new GMRESSolver(comm);
+      gmres_PMHSS->SetPrintLevel(1);
+      gmres_PMHSS->SetKDim(100);
+      gmres_PMHSS->SetMaxIter(100);
+      gmres_PMHSS->SetRelTol(1e-8);
+      gmres_PMHSS->SetAbsTol(0.0);
+
+      gmres_PMHSS->SetOperator(*aspacematrix_complex);
+      gmres_PMHSS->SetPreconditioner(*PMHSS);
+   }
+
+   gmres = new GMRESSolver(comm);
+   gmres->SetPrintLevel(1);
+   gmres->SetKDim(100);
+   gmres->SetMaxIter(100);
+   gmres->SetRelTol(1e-8);
+   gmres->SetAbsTol(0.0);
+
+   gmres->SetOperator(*matfree);
+   //gmres->SetPreconditioner(*aspacepc);
+   gmres->SetPreconditioner(
+      *PMHSS);  // TODO: which is better, PMHSS or gmres_PMHSS?
+   //gmres->SetPreconditioner(*gmres_PMHSS);
+
+   aspacewrapper = gmres;
+}
+
 void MatrixFreeAuxiliarySpace::SetupVCycle()
 {
    aspacewrapper = aspacepc;
 }
 
+// NOTE: if ess_tdof_list is empty, this just does the AMG mult.
 class ZeroWrap : public Solver
 {
 public:
@@ -360,6 +825,16 @@ private:
    Array<int>& ess_tdof_list;
 };
 
+void MatrixFreeAuxiliarySpace::SetupPMHSS()
+{
+   BlockDP = new BlockDiagonalPreconditioner(offsets);
+   BlockDP->SetDiagonalBlock(0, aspacepc);
+   BlockDP->SetDiagonalBlock(1, aspacepc);
+
+   PMHSS = new Complex_PMHSS(aspacematrix, aspacematrix_imag, BlockDP, NULL,
+                             1.0);
+}
+
 void MatrixFreeAuxiliarySpace::SetupAMG(int system_dimension)
 {
    if (system_dimension == 0)
@@ -397,15 +872,25 @@ void MatrixFreeAuxiliarySpace::SetupAMG(int system_dimension)
 
 void MatrixFreeAuxiliarySpace::Mult(const Vector& x, Vector& y) const
 {
-   int rank;
-   MPI_Comm_rank(comm, &rank);
-
-   y = 0.0;
-   aspacewrapper->Mult(x, y);
-   if (cg && rank == 0)
+   if (imagBdry)
    {
-      int q = cg->GetNumIterations();
-      inner_aux_iterations += q;
+      y = 0.0; // TODO: remove?
+      gmres->Mult(x, y);
+      //PMHSS->Mult(x, y); // TODO: solver?
+      //aspacewrapper->Mult(x, y);  // same as gmres?
+   }
+   else
+   {
+      int rank;
+      MPI_Comm_rank(comm, &rank);
+
+      y = 0.0;
+      aspacewrapper->Mult(x, y);
+      if (cg && rank == 0)
+      {
+         int q = cg->GetNumIterations();
+         inner_aux_iterations += q;
+      }
    }
 }
 
@@ -423,19 +908,26 @@ MatrixFreeAuxiliarySpace::~MatrixFreeAuxiliarySpace()
    inner iteration counts may need to be increased. Boundary conditions can
    matter as well (see DIAG_ZERO policy). */
 MatrixFreeAMS::MatrixFreeAMS(
-   ParBilinearForm& aform, Operator& oper, ParFiniteElementSpace& nd_fespace,
+   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,
 #endif
    int inner_pi_iterations, int inner_g_iterations, Solver * nd_smoother) :
-   Solver(oper.Height())
+   Solver((oper_complex != NULL) ? oper_complex->Height() : oper.Height())
 {
    int order = nd_fespace.GetFE(0)->GetOrder();
    ParMesh *mesh = nd_fespace.GetParMesh();
    int dim = mesh->Dimension();
 
+   const bool imagBdry = (beta_imag != NULL && abs_beta_imag != NULL);
+
+   MFEM_VERIFY((imagBdry && oper_complex != NULL) || (!imagBdry &&
+                                                      oper_complex == NULL), "");
+
    // smoother
    Array<int> ess_tdof_list;
    nd_fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
@@ -470,27 +962,85 @@ MatrixFreeAMS::MatrixFreeAMS(
 
    // build LOR space
    ParMesh mesh_lor(mesh, order, BasisType::GaussLobatto);
-
-   // build G space solver
-   Gspacesolver = new MatrixFreeAuxiliarySpace(mesh_lor, beta_coeff,
-                                               beta_mcoeff, ess_bdr, oper,
-                                               *Gradient,
+   if (imagBdry)
+   {
+      // build G space solver
+      Gspacesolver = new MatrixFreeAuxiliarySpace(mesh_lor, beta_coeff, beta_imag,
+                                                  abs_beta_imag,
+                                                  beta_mcoeff, ess_bdr, oper, oper_complex,
+                                                  *Gradient,
 #ifdef MFEM_USE_AMGX
-                                               useAmgX,
+                                                  useAmgX,
 #endif
-                                               inner_g_iterations);
+                                                  inner_g_iterations);
 
-   // build Pi space solver
-   Pispacesolver = new MatrixFreeAuxiliarySpace(mesh_lor, alpha_coeff,
-                                                beta_coeff, beta_mcoeff,
-                                                ess_bdr, oper, *Pi,
+      // build Pi space solver
+      Pispacesolver = new MatrixFreeAuxiliarySpace(mesh_lor, alpha_coeff,
+                                                   beta_coeff, beta_imag, abs_beta_imag, beta_mcoeff,
+                                                   ess_bdr, oper, oper_complex, *Pi,
 #ifdef MFEM_USE_AMGX
-                                                useAmgX,
+                                                   useAmgX,
 #endif
-                                                inner_pi_iterations);
+                                                   inner_pi_iterations);
+
+      offsets_nd.SetSize(3);
+      offsets_nd[0] = 0;
+      offsets_nd[1] = nd_fespace.GetTrueVSize();
+      offsets_nd[2] = offsets_nd[1];
+      offsets_nd.PartialSum();
+
+      offsets_vector.SetSize(3);
+      offsets_vector[0] = 0;
+      offsets_vector[1] = h1_fespace_d->GetTrueVSize();
+      offsets_vector[2] = offsets_vector[1];
+      offsets_vector.PartialSum();
+
+      offsets_scalar.SetSize(3);
+      offsets_scalar[0] = 0;
+      offsets_scalar[1] = h1_fespace->GetTrueVSize();
+      offsets_scalar[2] = offsets_scalar[1];
+      offsets_scalar.PartialSum();
+
+      Pi_block  = new BlockOperator(offsets_nd, offsets_vector);
+      Pi_block->SetDiagonalBlock(0, Pi.Ptr());
+      Pi_block->SetDiagonalBlock(1, Pi.Ptr());
+
+      Gradient_block  = new BlockOperator(offsets_nd, offsets_scalar);
+      Gradient_block->SetDiagonalBlock(0, Gradient.Ptr());
+      Gradient_block->SetDiagonalBlock(1, Gradient.Ptr());
+
+      smoother_block  = new BlockOperator(offsets_nd);
+      smoother_block->SetDiagonalBlock(0, smoother);
+      smoother_block->SetDiagonalBlock(1, smoother);
+
+      general_ams = new GeneralAMS(oper, oper_complex, *Pi_block, *Gradient_block,
+                                   *Pispacesolver,
+                                   *Gspacesolver, *smoother_block, ess_tdof_list);
+   }
+   else
+   {
+      // build G space solver
+      Gspacesolver = new MatrixFreeAuxiliarySpace(mesh_lor, beta_coeff,
+                                                  beta_mcoeff, ess_bdr, oper,
+                                                  *Gradient,
+#ifdef MFEM_USE_AMGX
+                                                  useAmgX,
+#endif
+                                                  inner_g_iterations);
+
+      // build Pi space solver
+      Pispacesolver = new MatrixFreeAuxiliarySpace(mesh_lor, alpha_coeff,
+                                                   beta_coeff, beta_mcoeff,
+                                                   ess_bdr, oper, *Pi,
+#ifdef MFEM_USE_AMGX
+                                                   useAmgX,
+#endif
+                                                   inner_pi_iterations);
+
+      general_ams = new GeneralAMS(oper, NULL, *Pi, *Gradient, *Pispacesolver,
+                                   *Gspacesolver, *smoother, ess_tdof_list);
+   }
 
-   general_ams = new GeneralAMS(oper, *Pi, *Gradient, *Pispacesolver,
-                                *Gspacesolver, *smoother, ess_tdof_list);
 
    delete h1_fec;
 }

linalg/auxiliary.hpp

@@ -57,14 +57,26 @@ public:
                              single V-cycle
    */
    MatrixFreeAuxiliarySpace(
-      ParMesh& mesh_lor, Coefficient* alpha_coeff,
-      Coefficient* beta_coeff, MatrixCoefficient* beta_mcoeff,
+      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.
@@ -89,6 +101,17 @@ public:
 #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;
@@ -107,17 +130,35 @@ private:
    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;
 };
 
 
@@ -132,6 +173,7 @@ public:
        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_,
@@ -147,6 +189,7 @@ public:
 
 private:
    const Operator& curlcurl_op;
+   Operator *oper_complex;
    const Operator& pi;
    const Operator& gradient;
    const Operator& pispacesolver;
@@ -194,9 +237,10 @@ public:
        @param nd_smoother  optional user-provided smoother for Nedelec space,
                            this object takes ownership and will delete.
    */
-   MatrixFreeAMS(ParBilinearForm& aform, Operator& oper,
+   MatrixFreeAMS(ParBilinearForm& aform, Operator& oper, Operator *oper_complex,
                  ParFiniteElementSpace& nd_fespace, Coefficient* alpha_coeff,
-                 Coefficient* beta_coeff, MatrixCoefficient* beta_mcoeff,
+                 Coefficient* beta_coeff, Coefficient* beta_imag,
+                 Coefficient* abs_beta_imag, MatrixCoefficient* beta_mcoeff,
                  Array<int>& ess_bdr,
 #ifdef MFEM_USE_AMGX
                  bool useAmgX = false,
@@ -224,6 +268,14 @@ private:
 
    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

linalg/operator.hpp

@@ -696,6 +696,42 @@ public:
    { A_.Mult(x, y); y *= a_; }
 };
 
+/// General sum operator: x -> A(x)+B(x)
+class SumOperator : public Operator
+{
+   const Operator *A, *B;
+   bool ownA, ownB;
+   mutable Vector z, w;
+   double cA, cB;
+
+public:
+   SumOperator(const Operator *A_, const Operator *B_,
+               bool ownA_, bool ownB_, double cA_, double cB_)
+      : Operator(A_->Height(), B_->Width()),
+        A(A_), B(B_), ownA(ownA_), ownB(ownB_), z(A_->Height()), w(A_->Width()),
+        cA(cA_), cB(cB_)
+   {
+      MFEM_VERIFY(A->Width() == B->Width() && A->Height() == B->Height(),
+                  "incompatible Operators: A->Width() = " << A->Width()
+                  << ", B->Height() = " << B->Height());
+
+      z.UseDevice(true);
+      w.UseDevice(true);
+   }
+
+   ~SumOperator()
+   {
+      if (ownA) { delete A; }
+      if (ownB) { delete B; }
+   }
+
+   virtual void Mult(const Vector &x, Vector &y) const
+   { B->Mult(x, z); A->Mult(x, y); y *= cA; z *= cB; y += z;}
+
+   virtual void MultTranspose(const Vector &x, Vector &y) const
+   { B->MultTranspose(x, w); A->MultTranspose(x, y); y *= cA; w *= cB; y += w;}
+
+};
 
 /** @brief The transpose of a given operator. Switches the roles of the methods
     Mult() and MultTranspose(). */

linalg/solvers.hpp

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