fixed bug in vertex-based-patched defined on the coarser grids

examples/ComplexMG/NDtest.cpp

@@ -0,0 +1,300 @@
+//                       MFEM Example 3 - Parallel Version
+//
+// Compile with: make ex3p
+//
+// Sample runs:  mpirun -np 4 ex3p -m ../data/star.mesh
+//               mpirun -np 4 ex3p -m ../data/square-disc.mesh -o 2
+//               mpirun -np 4 ex3p -m ../data/beam-tet.mesh
+//               mpirun -np 4 ex3p -m ../data/beam-hex.mesh
+//               mpirun -np 4 ex3p -m ../data/beam-hex.mesh -o 2 -pa
+//               mpirun -np 4 ex3p -m ../data/escher.mesh
+//               mpirun -np 4 ex3p -m ../data/escher.mesh -o 2
+//               mpirun -np 4 ex3p -m ../data/fichera.mesh
+//               mpirun -np 4 ex3p -m ../data/fichera-q2.vtk
+//               mpirun -np 4 ex3p -m ../data/fichera-q3.mesh
+//               mpirun -np 4 ex3p -m ../data/square-disc-nurbs.mesh
+//               mpirun -np 4 ex3p -m ../data/beam-hex-nurbs.mesh
+//               mpirun -np 4 ex3p -m ../data/amr-quad.mesh -o 2
+//               mpirun -np 4 ex3p -m ../data/amr-hex.mesh
+//               mpirun -np 4 ex3p -m ../data/ref-prism.mesh -o 1
+//               mpirun -np 4 ex3p -m ../data/octahedron.mesh -o 1
+//               mpirun -np 4 ex3p -m ../data/star-surf.mesh -o 2
+//               mpirun -np 4 ex3p -m ../data/mobius-strip.mesh -o 2 -f 0.1
+//               mpirun -np 4 ex3p -m ../data/klein-bottle.mesh -o 2 -f 0.1
+//
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include "mg.hpp"
+
+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 f_exact(const Vector &, Vector &);
+double freq = 1.0, kappa;
+int dim;
+
+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 visualization = true;
+   int href = 1;
+   int pref = 1;
+#ifdef MFEM_USE_AMGX
+   bool useAmgX = false;
+#endif
+
+   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(&href, "-gr", "--geometric-ref", "Number of Geometric refinements");
+   args.AddOption(&pref, "-or", "--order-ref", "Number of Order refinements");
+   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);
+      }
+      // HYPRE_Finalize();
+      MPI_Finalize();
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+   kappa = freq * M_PI;
+
+   Mesh *mesh = new Mesh(mesh_file, 1, 1);
+   dim = mesh->Dimension();
+   int sdim = mesh->SpaceDimension();
+
+   {
+      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.
+   ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
+   delete mesh;
+
+   // 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);
+
+   std::vector<ParFiniteElementSpace * > fespaces(href+1);
+   std::vector<FiniteElementCollection * > fecs(href+1);
+   std::vector<ParMesh * > ParMeshes(href+1);
+   std::vector<HypreParMatrix*>  P(href);
+   // for (int i = 0; i < href; i++)
+   // {
+   //    ParMeshes[i] = new ParMesh(*pmesh);
+   //    fespaces[i]  = new ParFiniteElementSpace(*fespace, *ParMeshes[i]);
+   //    pmesh->UniformRefinement();
+   //    // Update fespace
+   //    fespace->Update();
+   //    OperatorHandle Tr(Operator::Hypre_ParCSR);
+   //    fespace->GetTrueTransferOperator(*fespaces[i], Tr);
+   //    Tr.SetOperatorOwner(false);
+   //    Tr.Get(P[i]);
+   // }
+
+   // fecs[0]  = new ND_FECollection(order, dim);
+   fecs[0]  = new ND_FECollection(order, dim);
+   fespaces[0]  = new ParFiniteElementSpace(pmesh, fecs[0]);
+   VectorFunctionCoefficient Eex(sdim, E_exact);
+   ParGridFunction gf0(fespaces[0]);
+   gf0.ProjectCoefficient(Eex);
+
+   // cout << "order = " << order << endl;
+   Operator *Tr;
+   for (int i = 0; i < href; i++)
+   {
+      int new_order = std::pow(2,i+1); 
+      // cout << "new order = " << new_order << endl;
+      fecs[i+1] = new ND_FECollection(std::pow(2,i+1),dim);
+      fespaces[i+1]  = new ParFiniteElementSpace(pmesh,fecs[i+1]);
+      // cout << "test 0" << endl;
+      
+      Tr = new TrueTransferOperator(*fespaces[i], *fespaces[i+1]);
+      // cout << "test 1" << endl;
+   }
+
+   ParGridFunction gf1(fespaces[1]);
+
+   Tr->Mult(gf0,gf1);
+
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream sol_sock(vishost, visport);
+      sol_sock << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock.precision(8);
+      sol_sock << "solution\n" << *pmesh << gf0 << flush;
+      socketstream sol_sock1(vishost, visport);
+      sol_sock1 << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock1.precision(8);
+      sol_sock1 << "solution\n" << *pmesh << gf1 << flush;
+   }
+
+
+   return 0;
+
+   // fespaces[pref] = new ParFiniteElementSpace(*fespace);
+   FiniteElementCollection *fec = fecs[href];
+   ParFiniteElementSpace *fespace = fespaces[href];
+
+
+   HYPRE_BigInt size = fespace->GlobalTrueVSize();
+   if (myid == 0)
+   {
+      cout << "Number of finite element unknowns: " << size << endl;
+   }
+
+   Array<int> ess_tdof_list;
+   Array<int> ess_bdr;
+   if (pmesh->bdr_attributes.Size())
+   {
+      ess_bdr.SetSize(pmesh->bdr_attributes.Max());
+      ess_bdr = 1;
+      fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+
+   VectorFunctionCoefficient f(sdim, f_exact);
+   ParLinearForm *b = new ParLinearForm(fespace);
+   b->AddDomainIntegrator(new VectorFEDomainLFIntegrator(f));
+   b->Assemble();
+
+   ParGridFunction x(fespace);
+   VectorFunctionCoefficient E(sdim, E_exact);
+   x.ProjectCoefficient(E);
+
+   Coefficient *muinv = new ConstantCoefficient(1.0);
+   Coefficient *sigma = new ConstantCoefficient(1.0);
+   ParBilinearForm *a = new ParBilinearForm(fespace);
+   a->AddDomainIntegrator(new CurlCurlIntegrator(*muinv));
+   a->AddDomainIntegrator(new VectorFEMassIntegrator(*sigma));
+
+   a->Assemble();
+
+   OperatorPtr A;
+   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 (myid == 0)
+   {
+      cout << "Size of linear system: "
+            << A.As<HypreParMatrix>()->GetGlobalNumRows() << endl;
+   }
+
+   ParFiniteElementSpace *prec_fespace = fespace;
+   HypreAMS M(*A.As<HypreParMatrix>(), prec_fespace);
+
+
+   {
+      MGSolver M(A.As<HypreParMatrix>(),P,fespaces);
+      M.SetTheta(0.25);
+      CGSolver cg(MPI_COMM_WORLD);
+      cg.SetRelTol(1e-12);
+      cg.SetMaxIter(2000);
+      cg.SetPrintLevel(1);
+      cg.SetPreconditioner(M);
+      cg.SetOperator(*A);
+      cg.Mult(B, X);
+   }
+
+   a->RecoverFEMSolution(X, *b, x);
+
+   // 15. Compute and print the L^2 norm of the error.
+   {
+      double err = x.ComputeL2Error(E);
+      if (myid == 0)
+      {
+         cout << "\n|| E_h - E ||_{L^2} = " << err << '\n' << endl;
+      }
+   }
+
+   // 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);
+      sol_sock << "solution\n" << *pmesh << x << flush;
+   }
+
+   // 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 f_exact(const Vector &x, Vector &f)
+{
+   if (dim == 3)
+   {
+      f(0) = (1. + kappa * kappa) * sin(kappa * x(1));
+      f(1) = (1. + kappa * kappa) * sin(kappa * x(2));
+      f(2) = (1. + kappa * kappa) * 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; }
+   }
+}

examples/ComplexMG/complex_mgtest.cpp

@@ -0,0 +1,507 @@
+//
+// Description:  This example code demonstrates the use of MFEM to define a
+//               simple finite element discretization of the Helmholtz problem
+//               -Delta p - omega^2 p = f with impedance boundary conditiones.
+//
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include "mg.hpp"
+
+using namespace std;
+using namespace mfem;
+
+// Exact solution and r.h.s., see below for implementation.
+void get_helmholtz_solution_Re(const Vector &x, double & p, double dp[], double & d2p);
+void get_helmholtz_solution_Im(const Vector &x, double & p, double dp[], double & d2p);
+double p_exact_Re(const Vector &x);
+double p_exact_Im(const Vector &x);
+double f_exact_Re(const Vector &x);
+double f_exact_Im(const Vector &x);
+double g_exact_Re(const Vector &x);
+double g_exact_Im(const Vector &x);
+void grad_exact_Re(const Vector &x, Vector &grad_Re);
+void grad_exact_Im(const Vector &x, Vector &grad_Im);
+
+
+int dim;
+double omega;
+int sol = 1;
+double length = 1.0;
+
+int main(int argc, char *argv[])
+{
+
+   // 1. Initialise MPI
+   int num_procs, myid;
+   MPI_Init(&argc, &argv);                    // Initialise MPI
+   MPI_Comm_size(MPI_COMM_WORLD, &num_procs); //total number of processors available
+   MPI_Comm_rank(MPI_COMM_WORLD, &myid);      // Determine process identifier
+
+   //-----------------------------------------------------------------------------
+
+   // 2. Parse command-line options.
+   // geometry file
+   const char *mesh_file = "../../data/one-hex.mesh";
+   // finite element order of approximation
+   int order = 1;
+   bool visualization = 1;
+   // number of wavelengths
+   double k = 0.5;
+   // number of mg levels
+   int href = 1;
+   // number of initial ref
+   int initref = 1;
+   // dimension
+   int nd = 2;
+
+
+   // optional command line inputs
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   args.AddOption(&nd, "-nd", "--dim","Problem space dimension");               
+   args.AddOption(&sol, "-sol", "--exact",
+                  "Exact solution flag - 0:polynomial, 1: plane wave, -1: unknown exact");               
+   args.AddOption(&k, "-k", "--wavelengths",
+                  "Number of wavelengths.");
+   args.AddOption(&length, "-length", "--length",
+                  "length of the domainin in each direction.");                                       
+   args.AddOption(&href, "-href", "--href",
+                  "Number of Geometric Refinements.");
+   args.AddOption(&initref, "-initref", "--initref",
+                  "Number of initial refinements.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   // check if the inputs are correct
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      MPI_Finalize();
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+   // Angular frequency
+   omega = 2.0 * M_PI * k;
+
+   //-----------------------------------------------------------------------------
+
+   // if (scatter) pml = true; // for now only scattering problems with pml
+
+   // 3. Read the mesh from the given mesh file.
+   Mesh mesh;
+   
+   if (nd == 2) 
+   {
+      mesh = Mesh::MakeCartesian2D(1, 1, Element::QUADRILATERAL, true, length, length, false);
+
+   }
+   else
+   {
+      mesh = Mesh::MakeCartesian3D(1, 1, 1, Element::HEXAHEDRON, true, length, length, length);
+   }
+
+   dim = mesh.Dimension();
+   int sdim = mesh.SpaceDimension();
+
+   // 3. Executing uniform h-refinement
+   for (int i = 0; i < initref; i++ )
+   {
+      mesh.UniformRefinement();
+   }
+
+   // 5. Define a parallel mesh and delete the serial mesh.
+   ParMesh pmesh(MPI_COMM_WORLD, mesh);
+   mesh.Clear(); // the serial mesh is no longer needed
+   pmesh.UniformRefinement();
+
+   // ----------------------------------------------------------------------------
+
+   // 6. Define a finite element space on the mesh.
+   H1_FECollection fec(order, pmesh.Dimension());
+   ParFiniteElementSpace fespace(&pmesh, &fec);
+   std::vector<ParFiniteElementSpace * > fespaces(href+1);
+   std::vector<ParMesh * > ParMeshes(href+1);
+   std::vector<HypreParMatrix*>  P(href);
+   for (int i = 0; i < href; i++)
+   {
+      ParMeshes[i] = new ParMesh(pmesh);
+      fespaces[i]  = new ParFiniteElementSpace(fespace, *ParMeshes[i]);
+      pmesh.UniformRefinement();
+      // Update fespace
+      fespace.Update();
+      OperatorHandle Tr(Operator::Hypre_ParCSR);
+      fespace.GetTrueTransferOperator(*fespaces[i], Tr);
+      Tr.SetOperatorOwner(false);
+      Tr.Get(P[i]);
+   }
+   fespaces[href] = new ParFiniteElementSpace(fespace);
+
+   // 6. Set up the linear form (Real and Imaginary part)
+   FunctionCoefficient f_Re(f_exact_Re);
+   FunctionCoefficient g_Re(g_exact_Re);
+   VectorFunctionCoefficient grad_Re(sdim, grad_exact_Re);
+   FunctionCoefficient f_Im(f_exact_Im);
+   FunctionCoefficient g_Im(g_exact_Im);
+   VectorFunctionCoefficient grad_Im(sdim, grad_exact_Im);
+
+   // ParLinearForm *b_Re(new ParLinearForm);
+   ParComplexLinearForm b(&fespace, ComplexOperator::HERMITIAN);
+   b.AddDomainIntegrator(new DomainLFIntegrator(f_Re),new DomainLFIntegrator(f_Im));
+   if(sol >=0) // if exact solution exists. Otherwise use homogeneous impedence (gradp . n + i omega p = 0)
+   {
+      b.AddBoundaryIntegrator(new BoundaryNormalLFIntegrator(grad_Re),
+                              new BoundaryNormalLFIntegrator(grad_Im));
+      b.AddBoundaryIntegrator(new BoundaryLFIntegrator(g_Re),
+                              new BoundaryLFIntegrator(g_Im));
+   }
+   b.real().Vector::operator=(0.0);
+   b.imag().Vector::operator=(0.0);                        
+   b.Assemble();
+
+   // 7. Set up the bilinear form (Real and Imaginary part)
+   ConstantCoefficient one(1.0);
+   ConstantCoefficient zero(0.0);
+   ConstantCoefficient neg_omega(-pow(omega, 2));
+
+   ParSesquilinearForm a(&fespace,ComplexOperator::HERMITIAN);
+   ConstantCoefficient impedance(omega);
+
+   Array<int> bdr_attr(pmesh.bdr_attributes.Max());
+   bdr_attr = 1;
+   RestrictedCoefficient imp_rest(impedance,bdr_attr);
+
+   a.AddDomainIntegrator(new DiffusionIntegrator(one),NULL);
+   // Just putting imaginary integrator with zero just to keep sparsity pattern of real 
+   // and imaginary part the same 
+   a.AddDomainIntegrator(new MassIntegrator(neg_omega),new MassIntegrator(zero));
+   a.AddBoundaryIntegrator(new BoundaryMassIntegrator(zero),new BoundaryMassIntegrator(imp_rest));
+   a.Assemble(0);
+
+   Array<int> ess_tdof_list;
+   Array<int> ess_bdr(pmesh.bdr_attributes.Max());
+   ess_bdr = 0;
+   fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+
+   // Solution grid function
+   ParComplexGridFunction p_gf(&fespace);
+   ParComplexGridFunction p_gf_ex(&fespace);
+   FunctionCoefficient p_Re(p_exact_Re);
+   FunctionCoefficient p_Im(p_exact_Im);
+   p_gf = 0.0;
+   p_gf_ex.ProjectCoefficient(p_Re,p_Im);
+   if (sol >= 0) 
+   {
+      p_gf.ProjectBdrCoefficient(p_Re,p_Im,ess_bdr);
+   }
+
+   OperatorHandle Ah;
+   Vector X, B;
+   a.FormLinearSystem(ess_tdof_list, p_gf, b, Ah, X, B);
+
+   ComplexHypreParMatrix * AZ = Ah.As<ComplexHypreParMatrix>();
+
+
+   ComplexMGSolver * M = new ComplexMGSolver(AZ,P,fespaces);
+   M->SetTheta(0.1);
+
+
+   // ComplexSchwarzSmoother * M = new ComplexSchwarzSmoother(
+      // fespace.GetParMesh(),0,&fespace,AZ);
+
+   GMRESSolver gmres(MPI_COMM_WORLD);
+   gmres.SetPrintLevel(1);
+   gmres.SetMaxIter(2000);
+   gmres.SetKDim(200);
+   gmres.SetRelTol(1e-10);
+   gmres.SetAbsTol(0.0);
+   gmres.SetOperator(*AZ);
+   gmres.SetPreconditioner(*M);
+   gmres.Mult(B, X);
+
+
+   // {
+   //    ComplexMUMPSSolver cmumps;
+   //    cmumps.SetPrintLevel(0);
+   //    cmumps.SetOperator(*AZ);
+   //    cmumps.Mult(B,X);
+   // }
+
+
+   a.RecoverFEMSolution(X,B,p_gf);
+
+
+   if (sol >= 0)
+   {
+      const int h1_norm_type = 1;
+      double L2error;
+      double H1error;
+
+      double L2err_Re  = p_gf.real().ComputeL2Error(p_Re);
+      double L2err_Im  = p_gf.imag().ComputeL2Error(p_Im);
+
+      double loc_H1err_Re = p_gf.real().ComputeH1Error(&p_Re, &grad_Re, &one, 1.0, h1_norm_type);
+      double loc_H1err_Im = p_gf.imag().ComputeH1Error(&p_Im, &grad_Im, &one, 1.0, h1_norm_type);
+      double H1err_Re = GlobalLpNorm(2.0, loc_H1err_Re, MPI_COMM_WORLD);
+      double H1err_Im = GlobalLpNorm(2.0, loc_H1err_Im, MPI_COMM_WORLD);
+         
+      L2error = sqrt(L2err_Re*L2err_Re + L2err_Im*L2err_Im);
+      H1error = sqrt(H1err_Re*H1err_Re + H1err_Im*H1err_Im);
+      if (myid == 0)
+      {
+         cout << " || p_h - p ||_{H^1} = " << H1error <<  endl;
+         cout << " || p_h - p ||_{L^2} = " << L2error <<  endl;
+      }
+   }
+
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      string keys;
+      if(dim ==2 )
+      {
+         keys = "keys mrRljc\n";
+      }
+      else
+      {
+         keys = "keys mc\n";
+      }
+      socketstream sol_sock_re(vishost, visport);
+      sol_sock_re << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock_re.precision(8);
+      sol_sock_re << "solution\n" << pmesh << p_gf.real() << "window_title 'Numerical Pressure (real part)' " 
+               << keys << flush;
+      socketstream sol_sock_im(vishost, visport);
+      sol_sock_im << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock_im.precision(8);
+      sol_sock_im << "solution\n" << pmesh << p_gf.imag() << "window_title 'Numerical Pressure (imag part)' " 
+               << keys << flush;         
+   }
+   MPI_Finalize();
+   return 0;
+}
+
+//define exact solutions
+void get_helmholtz_solution_Re(const Vector &x, double & p, double dp[], double & d2p)
+{
+
+   if (sol == 0) // polynomial
+   {  
+      if (dim == 3)
+      {
+         p = x[0]*(1.0 - x[0]) * x[1]*(1.0 - x[1]) * x[2]*(1.0 - x[2]);
+         dp[0] = (1.0 - 2.0 *x[0]) * x[1]*(1.0 - x[1]) * x[2]*(1.0 - x[2]);
+         dp[1] = (1.0 - 2.0 *x[1]) * x[0]*(1.0 - x[0]) * x[2]*(1.0 - x[2]);
+         dp[2] = (1.0 - 2.0 *x[2]) * x[0]*(1.0 - x[0]) * x[1]*(1.0 - x[1]);
+         d2p = -2.0*(-1.0 + x[0]) * x[0] * (-1.0 + x[1]) * x[1] 
+               -2.0*(-1.0 + x[0]) * x[0] * (-1.0 + x[2]) * x[2] 
+               -2.0*(-1.0 + x[1]) * x[1] * (-1.0 + x[2]) * x[2];   
+      }
+      else
+      {
+         p = x[1] * (1.0 - x[1])* x[0] * (1.0 - x[0]);
+         dp[0] = (1.0 - 2.0 *x[0]) * x[1]*(1.0 - x[1]);
+         dp[1] = (1.0 - 2.0 *x[1]) * x[0]*(1.0 - x[0]);
+         d2p = - 2.0 * x[1] * (1.0 - x[1]) 
+               - 2.0 * x[0] * (1.0 - x[0]);   
+      }
+   }
+
+   else if(sol == 1) // Plane wave
+   {
+      double alpha;
+      if (dim == 2)
+      {
+         alpha = omega/sqrt(2);
+         p = cos(alpha * ( x(0) + x(1) ) );
+         dp[0] = -alpha * sin(alpha * ( x(0) + x(1) ) );
+         dp[1] = dp[0];
+         d2p = -2.0 * alpha * alpha * p;
+      }
+      else 
+      {
+         alpha = omega/sqrt(3);
+         p = cos(alpha * ( x(0) + x(1) + x(2) ) );
+         dp[0] = -alpha * sin(alpha * ( x(0) + x(1) + x(2) ) );
+         dp[1] = dp[0];
+         dp[2] = dp[0];
+         d2p = -3.0 * alpha * alpha * p;
+      }
+   }
+   else if (sol == 2)
+   {
+      if (dim == 2 )
+      {
+         // shift to avoid singularity
+         double shift = 0.1;
+         double x0 = x(0) + shift;
+         double x1 = x(1) + shift;
+         //
+         double r = sqrt(x0 * x0 + x1 * x1);
+
+         p = cos(omega * r);
+         double r_x = x0 / r;
+         double r_y = x1 / r;
+         double r_xx = (1.0 / r) * (1.0 - r_x * r_x);
+         double r_yy = (1.0 / r) * (1.0 - r_y * r_y);
+
+         dp[0] = - omega * sin(omega * r) * r_x; 
+         dp[1] = - omega * sin(omega * r) * r_y; 
+
+         d2p = -omega*omega * cos(omega * r)*r_x * r_x - omega * sin(omega*r) * r_xx
+               -omega*omega * cos(omega * r)*r_y * r_y - omega * sin(omega*r) * r_yy;
+      }
+      else
+      {
+         // shift to avoid singularity
+         double shift = 0.1;
+         double x0 = x(0) + shift;
+         double x1 = x(1) + shift;
+         double x2 = x(2) + shift;
+         //
+         double r = sqrt(x0 * x0 + x1 * x1 + x2 * x2);
+
+         p = cos(omega * r);
+
+         double r_x = x0 / r;
+         double r_y = x1 / r;
+         double r_z = x2 / r;
+         double r_xx = (1.0 / r) * (1.0 - r_x * r_x);
+         double r_yy = (1.0 / r) * (1.0 - r_y * r_y);
+         double r_zz = (1.0 / r) * (1.0 - r_z * r_z);
+
+         dp[0] = - omega * sin(omega * r) * r_x; 
+         dp[1] = - omega * sin(omega * r) * r_y; 
+         dp[2] = - omega * sin(omega * r) * r_z; 
+
+         d2p = -omega*omega * cos(omega * r)*r_x * r_x - omega * sin(omega*r) * r_xx
+               -omega*omega * cos(omega * r)*r_y * r_y - omega * sin(omega*r) * r_yy
+               -omega*omega * cos(omega * r)*r_z * r_z - omega * sin(omega*r) * r_zz;
+      }
+   }
+}
+
+void get_helmholtz_solution_Im(const Vector &x, double & p, double dp[], double & d2p)
+{
+   if (sol == 0) // polynomial
+   {  
+      if (dim == 3)
+      {
+         p = x[0]*(1.0 - x[0]) * x[1]*(1.0 - x[1]) * x[2]*(1.0 - x[2]);
+         dp[0] = (1.0 - 2.0 *x[0]) * x[1]*(1.0 - x[1]) * x[2]*(1.0 - x[2]);
+         dp[1] = (1.0 - 2.0 *x[1]) * x[0]*(1.0 - x[0]) * x[2]*(1.0 - x[2]);
+         dp[2] = (1.0 - 2.0 *x[2]) * x[0]*(1.0 - x[0]) * x[1]*(1.0 - x[1]);
+         d2p = -2.0*(-1.0 + x[0]) * x[0] * (-1.0 + x[1]) * x[1] 
+               -2.0*(-1.0 + x[0]) * x[0] * (-1.0 + x[2]) * x[2] 
+               -2.0*(-1.0 + x[1]) * x[1] * (-1.0 + x[2]) * x[2];   
+      }
+      else
+      {
+         p = x[1] * (1.0 - x[1])* x[0] * (1.0 - x[0]);
+         dp[0] = (1.0 - 2.0 *x[0]) * x[1]*(1.0 - x[1]);
+         dp[1] = (1.0 - 2.0 *x[1]) * x[0]*(1.0 - x[0]);
+         d2p = - 2.0 * x[1] * (1.0 - x[1]) 
+               - 2.0 * x[0] * (1.0 - x[0]);   
+      }
+   }
+   else if (sol == 1)// plane wave
+   {
+      double alpha; 
+      if (dim == 2)
+      {
+         alpha = omega/sqrt(2);
+         p = -sin(alpha * ( x(0) + x(1) ) );
+         dp[0] = -alpha * cos(alpha * ( x(0) + x(1) ) );
+         dp[1] = dp[0];
+         d2p = -2.0 * alpha * alpha * p;
+      }
+      else 
+      {
+         alpha = omega/sqrt(3);
+         p = -sin(alpha * ( x(0) + x(1) + x(2) ) );
+         dp[0] = -alpha * cos(alpha * ( x(0) + x(1) + x(2) ) );
+         dp[1] = dp[0];
+         dp[2] = dp[0];
+         d2p = -3.0 * alpha * alpha * p;
+      }   
+   }
+}
+
+
+double p_exact_Re(const Vector &x)
+{
+   double p, d2p;
+   double dp[3];
+   get_helmholtz_solution_Re(x, p, dp, d2p);
+   return p;
+}
+double p_exact_Im(const Vector &x)
+{
+   double p, d2p;
+   double dp[3];
+   get_helmholtz_solution_Im(x, p, dp, d2p);
+   return p;
+}
+
+//calculate RHS from exact solution f = - \Delta u
+double f_exact_Re(const Vector &x)
+{
+   double p_re, d2p_re;
+   double dp_re[3];
+   double f_re;
+   f_re = 0.0;
+   get_helmholtz_solution_Re(x, p_re, dp_re, d2p_re);
+   f_re = -d2p_re - omega * omega * p_re;
+   return f_re;
+}
+double f_exact_Im(const Vector &x)
+{
+   double p_im, d2p_im;
+   double dp_im[3];
+   double f_im;
+   f_im = 0.0;
+   get_helmholtz_solution_Im(x, p_im, dp_im, d2p_im);
+   f_im = -d2p_im - omega * omega * p_im;
+   return f_im;   
+}
+
+void grad_exact_Re(const Vector &x, Vector &dp)
+{
+   double p, d2p;
+   get_helmholtz_solution_Re(x, p, dp, d2p);
+}
+void grad_exact_Im(const Vector &x, Vector &dp)
+{
+   double p, d2p;
+   get_helmholtz_solution_Im(x, p, dp, d2p);
+}
+
+//define impedence coefficient: i omega p
+double g_exact_Re(const Vector &x)
+{
+   double p, d2p;
+   double dp[3];
+   get_helmholtz_solution_Im(x, p, dp, d2p);
+
+   return -omega * p;
+}
+double g_exact_Im(const Vector &x)
+{
+   double p, d2p;
+   double dp[3];
+   get_helmholtz_solution_Re(x, p, dp, d2p);
+
+   return omega * p;
+}

examples/ComplexMG/makefile

@@ -0,0 +1,58 @@
+# Copyright (c) 2010-2021, 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.
+
+# Use the MFEM build directory
+MFEM_DIR ?= ../..
+MFEM_BUILD_DIR ?= ../..
+SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/examples/ComplexMG/,)
+CONFIG_MK = $(MFEM_BUILD_DIR)/config/config.mk
+
+MFEM_LIB_FILE = mfem_is_not_built
+-include $(CONFIG_MK)
+
+SEQ_EXAMPLES = 
+PAR_EXAMPLES = mgtest complex_mgtest
+ifeq ($(MFEM_USE_MPI),NO)
+   EXAMPLES = $(SEQ_EXAMPLES)
+else
+   EXAMPLES = $(PAR_EXAMPLES) $(SEQ_EXAMPLES)
+endif
+
+.SUFFIXES:
+.SUFFIXES: .o .cpp .mk
+.PHONY: all clean
+.PRECIOUS: %.o
+
+COMMON_O= schwarz.o mg.o
+
+# Remove built-in rules
+%: %.cpp
+%.o: %.cpp
+
+all: $(EXAMPLES)
+
+# Rules for building the EXAMPLES
+
+%: $(SRC)%.cpp $(COMMON_O) $(MFEM_LIB_FILE) $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) $< -o $@ $(COMMON_O) $(MFEM_LIBS)
+
+# Rules for compiling miniapp dependencies
+$(COMMON_O) $($(EXAMPLES)): \
+%.o: $(SRC)%.cpp $(SRC)%.hpp $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) -c $(<) -o $(@)
+
+# Generate an error message if the MFEM library is not built and exit
+$(MFEM_LIB_FILE):
+	$(error The MFEM library is not built)
+
+clean:
+	rm -f *.o *~ $(SEQ_EXAMPLES) $(PAR_EXAMPLES)
+	rm -rf *.dSYM *.TVD.*breakpoints

examples/ComplexMG/mg.cpp

@@ -0,0 +1,225 @@
+#include "mfem.hpp"
+#include "mg.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+
+MGSolver::MGSolver(HypreParMatrix * Af_, std::vector<HypreParMatrix *> P_,std::vector<ParFiniteElementSpace * > fespaces)
+   : Solver(Af_->Height(), Af_->Width()), Af(Af_), P(P_) {
+
+   StopWatch chrono;
+   NumGrids = P.size();
+   S.resize(NumGrids);
+   A.resize(NumGrids + 1);
+
+   A[NumGrids] = Af;
+   for (int i = NumGrids ; i > 0; i--)
+   {
+      A[i - 1] = RAP(A[i], P[i - 1]);
+   }
+
+   mumps = new MUMPSSolver;
+   mumps->SetPrintLevel(0);
+   mumps->SetMatrixSymType(MUMPSSolver::MatType::UNSYMMETRIC);
+   mumps->SetOperator(*A[0]);
+   invAc = mumps;
+
+
+   for (int i = NumGrids - 1; i >= 0 ; i--)
+   {
+      // S[i] = new SchwarzSmoother(fespaces[i]->GetParMesh(),1,fespaces[i+1],A[i+1]);
+      S[i] = new SchwarzSmoother(fespaces[i+1]->GetParMesh(),0,fespaces[i+1],A[i+1]);
+   }
+}
+
+void MGSolver::Mult(const Vector &r, Vector &z) const
+{
+   // Residual vectors
+   std::vector<Vector> rv(NumGrids + 1);
+   // correction vectors
+   std::vector<Vector> zv(NumGrids + 1);
+   // allocation
+   for (int i = 0; i <= NumGrids ; i++)
+   {
+      int n = A[i]->Width();
+      rv[i].SetSize(n);
+      zv[i].SetSize(n);
+   }
+   // Initial residual
+   rv[NumGrids] = r;
+   // smooth and update residuals down to the coarsest level
+   for (int i = NumGrids; i > 0 ; i--)
+   {
+      // Pre smooth
+      S[i - 1]->Mult(rv[i], zv[i]); zv[i] *= theta;
+      // compute residual
+      Vector w(A[i]->Height());
+      A[i]->Mult(zv[i], w);
+      rv[i] -= w;
+      // Restrict
+      P[i - 1]->MultTranspose(rv[i], rv[i - 1]);
+   }
+
+   // Coarse grid Solve
+   invAc->Mult(rv[0], zv[0]);
+   //
+   for (int i = 1; i <= NumGrids ; i++)
+   {
+      // Prolong correction
+      Vector u(P[i - 1]->Height());
+      P[i - 1]->Mult(zv[i - 1], u);
+      // Update correction
+      zv[i] += u;
+      // Update residual
+      Vector v(A[i]->Height());
+      A[i]->Mult(u, v); rv[i] -= v;
+      // Post smooth
+      S[i - 1]->Mult(rv[i], v); v *= theta;
+      // Update correction
+      zv[i] += v;
+   }
+   z = zv[NumGrids];
+}
+
+MGSolver::~MGSolver() 
+{
+   int n = S.size();
+   for (int i = n - 1; i >= 0 ; i--)
+   {
+      delete S[i];
+      delete A[i];
+   }
+   S.clear();
+   A.clear();
+   delete invAc;
+}
+
+ComplexMGSolver::ComplexMGSolver(ComplexHypreParMatrix * Af_, 
+std::vector<HypreParMatrix *> P_,std::vector<ParFiniteElementSpace * > fespaces)
+   : Solver(Af_->Height(), Af_->Width()), Af(Af_), P(P_) {
+
+   NumGrids = P.size();
+   S.resize(NumGrids);
+   A.resize(NumGrids + 1);
+
+   A[NumGrids] = Af;
+   for (int i = NumGrids ; i > 0; i--)
+   {
+      // A[i - 1] = RAP(A[i], P[i - 1]);
+      A[i - 1] = new ComplexHypreParMatrix(RAP(&A[i]->real(), P[i - 1]), 
+                                           RAP(&A[i]->imag(), P[i - 1]), true, true);
+   }
+
+   mumps = new ComplexMUMPSSolver;
+   mumps->SetPrintLevel(0);
+   mumps->SetOperator(*A[0]);
+   invAc = mumps;
+
+
+   for (int i = NumGrids - 1; i >= 0 ; i--)
+   {
+      S[i] = new ComplexSchwarzSmoother(fespaces[i]->GetParMesh(),1,fespaces[i+1],A[i+1]);
+      // S[i] = new ComplexSchwarzSmoother(fespaces[i+1]->GetParMesh(),0,fespaces[i+1],A[i+1]);
+   }
+
+}
+
+void ComplexMGSolver::Mult(const Vector &r, Vector &z) const
+{
+   // Residual vectors
+   std::vector<Vector> rv(NumGrids + 1);
+   // correction vectors
+   std::vector<Vector> zv(NumGrids + 1);
+   // allocation
+   for (int i = 0; i <= NumGrids ; i++)
+   {
+      int n = A[i]->Width();
+      rv[i].SetSize(n);
+      zv[i].SetSize(n);
+   }
+   // Initial residual
+   rv[NumGrids] = r;
+   // smooth and update residuals down to the coarsest level
+   Vector rv0, rv1;
+   for (int i = NumGrids; i > 0 ; i--)
+   {
+      // Pre smooth
+      S[i - 1]->Mult(rv[i], zv[i]); zv[i] *= theta;
+      // compute residual
+      Vector w(A[i]->Height());
+      A[i]->Mult(zv[i], w);
+      rv[i] -= w;
+      // Restrict
+
+      double * data1 = rv[i].GetData();
+      double * data0 = rv[i-1].GetData();
+      int size1 = rv[i].Size();
+      int size0 = rv[i-1].Size();
+
+      // Real part
+      rv1.SetDataAndSize(data1, size1/2);
+      rv0.SetDataAndSize(data0, size0/2);
+      P[i - 1]->MultTranspose(rv1, rv0);
+
+      // Imag part
+      rv1.SetDataAndSize(&data1[size1/2], size1/2);
+      rv0.SetDataAndSize(&data0[size0/2], size0/2);
+      P[i - 1]->MultTranspose(rv1, rv0);
+      // P[i - 1]->MultTranspose(rv[i], rv[i - 1]);
+   }
+
+   // Coarse grid Solve
+   invAc->Mult(rv[0], zv[0]);
+   //
+   Vector u0,u1;
+   for (int i = 1; i <= NumGrids ; i++)
+   {
+      // Prolong correction
+      Vector u(2*P[i - 1]->Height());
+      // real part
+      double * data1 = u.GetData();
+      double * data0 = zv[i - 1].GetData();
+      int size1 = u.Size();
+      int size0 = zv[i-1].Size();
+
+      u1.SetDataAndSize(data1, size1/2);
+      u0.SetDataAndSize(data0, size0/2);
+
+      P[i - 1]->Mult(u0, u1);
+
+      // imag part
+
+      u1.SetDataAndSize(&data1[size1/2], size1/2);
+      u0.SetDataAndSize(&data0[size0/2], size0/2);
+
+      P[i - 1]->Mult(u0, u1);
+
+      // P[i - 1]->Mult(zv[i - 1], u);
+      // Update correction
+      zv[i] += u;
+      // Update residual
+      Vector v(A[i]->Height());
+      A[i]->Mult(u, v); rv[i] -= v;
+      // Post smooth
+      S[i - 1]->Mult(rv[i], v); v *= theta;
+      // Update correction
+      zv[i] += v;
+   }
+   z = zv[NumGrids];
+}
+
+ComplexMGSolver::~ComplexMGSolver() 
+{
+   int n = S.size();
+   for (int i = n - 1; i >= 0 ; i--)
+   {
+      delete S[i];
+      delete A[i];
+   }
+   S.clear();
+   A.clear();
+   delete invAc;
+}

examples/ComplexMG/mg.hpp

@@ -0,0 +1,55 @@
+
+   
+#pragma once
+
+#include "mfem.hpp"
+#include "schwarz.hpp"
+
+class MGSolver : public Solver
+{
+private:
+   /// The linear system matrix
+   HypreParMatrix *Af;
+   std::vector<HypreParMatrix *> A;
+   std::vector<HypreParMatrix *> P;
+   std::vector<SchwarzSmoother *> S;
+   int NumGrids;
+   MUMPSSolver *mumps = nullptr;
+   Solver *invAc = nullptr;
+   double theta = 1.0;
+
+public:
+   MGSolver(HypreParMatrix *Af_, std::vector<HypreParMatrix *> P_, std::vector<ParFiniteElementSpace *> fespaces);
+
+   virtual void SetOperator(const Operator &op) {}
+
+   virtual void SetTheta(const double a) { theta = a; }
+
+   virtual void Mult(const Vector &r, Vector &z) const;
+   virtual ~MGSolver();
+};
+
+class ComplexMGSolver : public Solver
+{
+private:
+   /// The linear system matrix
+   ComplexHypreParMatrix *Af;
+   std::vector<ComplexHypreParMatrix *> A;
+   std::vector<HypreParMatrix *> P;
+   std::vector<ComplexSchwarzSmoother *> S;
+   int NumGrids;
+   ComplexMUMPSSolver *mumps = nullptr;
+   Solver *invAc = nullptr;
+   double theta = 1.0;
+
+public:
+   ComplexMGSolver(ComplexHypreParMatrix *Af_, 
+   std::vector<HypreParMatrix *> P_, std::vector<ParFiniteElementSpace *> fespaces);
+
+   virtual void SetOperator(const Operator &op) {}
+
+   virtual void SetTheta(const double a) { theta = a; }
+
+   virtual void Mult(const Vector &r, Vector &z) const;
+   virtual ~ComplexMGSolver();
+};

examples/ComplexMG/mgtest.cpp

@@ -0,0 +1,138 @@
+//                       MFEM Example 0 - Parallel Version
+//
+// Compile with: make ex0p
+//
+// Sample runs:  mpirun -np 4 ex0p
+//               mpirun -np 4 ex0p -m ../data/fichera.mesh
+//               mpirun -np 4 ex0p -m ../data/square-disc.mesh -o 2
+//
+// Description: This example code demonstrates the most basic parallel usage of
+//              MFEM to define a simple finite element discretization of the
+//              Laplace problem -Delta u = 1 with zero Dirichlet boundary
+//              conditions. General 2D/3D serial mesh files and finite element
+//              polynomial degrees can be specified by command line options.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+#include "mg.hpp"
+
+using namespace std;
+using namespace mfem;
+
+int main(int argc, char *argv[])
+{
+   // 1. Initialize MPI
+   MPI_Session mpi(argc, argv);
+
+   // 2. Parse command line options
+   const char *mesh_file = "../data/star.mesh";
+   int order = 1;
+   int href = 1;
+   int pref = 0;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order", "Finite element polynomial degree");
+   args.AddOption(&href, "-gr", "--geometric-ref", "Number of Geometric refinements");
+   args.AddOption(&pref, "-or", "--order-ref", "Number of Order refinements");
+   args.ParseCheck();
+
+   // 3. Read the serial mesh from the given mesh file.
+   Mesh mesh(mesh_file);
+
+   // 4. Define a parallel mesh by a partitioning of the serial mesh. Refine
+   //    this mesh once in parallel to increase the resolution.
+   ParMesh pmesh(MPI_COMM_WORLD, mesh);
+   mesh.Clear(); // the serial mesh is no longer needed
+   pmesh.UniformRefinement();
+
+   // 5. Define a finite element space on the mesh. Here we use H1 continuous
+   //    high-order Lagrange finite elements of the given order.
+   H1_FECollection fec(order, pmesh.Dimension());
+   ParFiniteElementSpace fespace(&pmesh, &fec);
+
+
+   std::vector<ParFiniteElementSpace * > fespaces(href+1);
+   std::vector<ParMesh * > ParMeshes(href+1);
+   std::vector<HypreParMatrix*>  P(href);
+   for (int i = 0; i < href; i++)
+   {
+      ParMeshes[i] = new ParMesh(pmesh);
+      fespaces[i]  = new ParFiniteElementSpace(fespace, *ParMeshes[i]);
+      pmesh.UniformRefinement();
+      // Update fespace
+      fespace.Update();
+      OperatorHandle Tr(Operator::Hypre_ParCSR);
+      fespace.GetTrueTransferOperator(*fespaces[i], Tr);
+      Tr.SetOperatorOwner(false);
+      Tr.Get(P[i]);
+   }
+   fespaces[href] = new ParFiniteElementSpace(fespace);
+
+
+   HYPRE_BigInt total_num_dofs = fespace.GlobalTrueVSize();
+   if (mpi.Root()) { cout << "Number of unknowns: " << total_num_dofs << endl; }
+
+   // 6. Extract the list of all the boundary DOFs. These will be marked as
+   //    Dirichlet in order to enforce zero boundary conditions.
+   Array<int> boundary_dofs;
+   fespace.GetBoundaryTrueDofs(boundary_dofs);
+
+   // 7. Define the solution x as a finite element grid function in fespace. Set
+   //    the initial guess to zero, which also sets the boundary conditions.
+   ParGridFunction x(&fespace);
+   x = 0.0;
+
+   // 8. Set up the linear form b(.) corresponding to the right-hand side.
+   ConstantCoefficient one(1.0);
+   ParLinearForm b(&fespace);
+   b.AddDomainIntegrator(new DomainLFIntegrator(one));
+   b.Assemble();
+
+   // 9. Set up the bilinear form a(.,.) corresponding to the -Delta operator.
+   ParBilinearForm a(&fespace);
+   a.AddDomainIntegrator(new DiffusionIntegrator);
+   a.Assemble();
+
+   // 10. Form the linear system A X = B. This includes eliminating boundary
+   //     conditions, applying AMR constraints, parallel assembly, etc.
+   HypreParMatrix A;
+   Vector B, X;
+   a.FormLinearSystem(boundary_dofs, x, b, A, X, B);
+
+   // 11. Solve the system using PCG with hypre's BoomerAMG preconditioner.
+   // HypreBoomerAMG M(A);
+   // MGSolver * M = new MGSolver(&A,P,fespaces);
+   // M->SetTheta(0.2);
+
+
+   int j = 3;
+   ParMesh * cpmesh = fespaces[href-j]->GetParMesh();
+   SchwarzSmoother * M = new SchwarzSmoother(cpmesh, j, &fespace, &A);
+   M->SetDumpingParam(0.3);
+   M->SetNumSmoothSteps(1);
+
+   X = 0.;
+   CGSolver cg(MPI_COMM_WORLD);
+   cg.SetRelTol(1e-8);
+   cg.SetMaxIter(2000);
+   cg.SetPrintLevel(1);
+   cg.SetPreconditioner(*M);
+   cg.SetOperator(A);
+   cg.Mult(B, X);
+
+   // 12. Recover the solution x as a grid function and save to file. The output
+   //     can be viewed using GLVis as follows: "glvis -np <np> -m mesh -g sol"
+   a.RecoverFEMSolution(X, b, x);
+
+   char vishost[] = "localhost";
+   int  visport   = 19916;
+   socketstream sol_sock(vishost, visport);
+   sol_sock << "parallel " << mpi.WorldSize() << " " << mpi.WorldRank() << "\n";
+   sol_sock.precision(8);
+   sol_sock << "solution\n" << pmesh << x << flush;
+
+
+   return 0;
+}

examples/ComplexMG/schwarz.hpp

@@ -0,0 +1,159 @@
+#pragma once
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+bool is_a_patch(int iv, Array<int> patch_ids);
+
+bool owned(int tdof, int * offs);
+
+SparseMatrix * GetLocalRestriction(const Array<int> & tdof_i, const int * row_start, 
+                                   int num_rows, int num_cols);                              
+
+void GetLocal2GlobalMap(const Array<int> & tdof_i, const int * row_start, 
+                        int num_rows, int num_cols, Array<int> & l2gmap);
+
+void GetOffdColumnValues(const Array<int> & tdof_i, const Array<int> & tdof_j, SparseMatrix & offd, const int * cmap, 
+                         const int * row_start, SparseMatrix * PatchMat);
+
+void GetArrayIntersection(const Array<int> & A, const Array<int> & B, Array<int>  & C); 
+
+
+
+int GetNumColumns(const int tdof_i, const Array<int> & tdof_j, SparseMatrix & diag,
+SparseMatrix & offd, const int * cmap, const int * row_start);
+
+void GetColumnValues(int tdof_i,const Array<int> & tdof_j, SparseMatrix & diag ,
+SparseMatrix & offd, const int *cmap, const int * row_start, Array<int> &cols, Array<double> &vals);
+
+
+class VertexPatchInfo 
+{
+private:
+   ParMesh * pmesh = nullptr;
+   int ref_levels=0;
+public:
+   int mynrpatch;
+   int nrpatch;
+   std::vector<Array<int>> vert_contr;
+   std::vector<Array<int>> edge_contr;
+   std::vector<Array<int>> face_contr;
+   std::vector<Array<int>> elem_contr;
+   Array<int> host_rank;
+   Array<int> patch_natural_order_idx;
+   Array<int> patch_global_dofs_ids;
+
+   VertexPatchInfo(ParMesh * pmesh_, int ref_levels_);
+   ~VertexPatchInfo() {}
+
+};
+
+class PatchDofInfo 
+{
+public:
+   MPI_Comm comm = MPI_COMM_WORLD;
+   int nrpatch;
+   Array<int> host_rank;
+   std::vector<Array<int>> patch_tdofs;
+   std::vector<Array<int>> patch_local_tdofs;
+   PatchDofInfo(ParMesh * cpmesh_, int ref_levels_,ParFiniteElementSpace *fespace);
+   ~PatchDofInfo(){};
+};
+
+class PatchAssembly
+{
+public:
+   MPI_Comm comm;
+   int nrpatch;
+   std::vector<int>tdof_offsets;
+   std::vector<Array<int>> patch_other_tdofs;
+   std::vector<Array<int>> patch_owned_other_tdofs;
+   std::vector<Array<int>> l2gmaps; // patch to global maps for the dofs owned by the processor
+   Array<SparseMatrix* > PatchMat;
+   PatchDofInfo *patch_tdof_info=nullptr;   
+   Array<int> host_rank; 
+   HypreParMatrix * A = nullptr;
+   int get_rank(int tdof);
+   // constructor
+   PatchAssembly(ParMesh * cpmesh_, int ref_levels_,ParFiniteElementSpace *fespace_, HypreParMatrix * A_);
+   ~PatchAssembly();
+private:
+   void compute_trueoffsets();
+   ParFiniteElementSpace *fespace=nullptr;
+};
+
+class PatchRestriction  {
+private:
+   MPI_Comm comm;
+   int num_procs, myid;
+   Array<int> host_rank;
+   PatchAssembly * P;
+   int nrpatch;
+   Array<int> send_count;
+   Array<int> send_displ;
+   Array<int> recv_count;
+   Array<int> recv_displ;
+   int sbuff_size;
+   int rbuff_size;
+public:
+   PatchRestriction(PatchAssembly * P_);
+   void Mult(const Vector & r , Array<BlockVector *> & res);
+   void MultTranspose(const Array<BlockVector *> & sol, Vector & z);
+   virtual ~PatchRestriction() {}
+};
+
+class SchwarzSmoother : public Solver {
+private:
+   MPI_Comm comm;
+   int nrpatch;
+   int maxit = 1;
+   double theta = 1.0;
+   Array<int> host_rank;
+#ifdef MFEM_USE_SUITESPARSE
+   Array<UMFPackSolver * > PatchInv;
+#else   
+   Array<GMRESSolver * > PatchInv;
+#endif   /// The linear system matrix
+   HypreParMatrix * A;
+   PatchAssembly * P;
+   PatchRestriction * R= nullptr;
+public:
+   SchwarzSmoother(ParMesh * cpmesh_, int ref_levels_,ParFiniteElementSpace *fespace_, HypreParMatrix * A_);
+   void SetNumSmoothSteps(const int iter) {maxit = iter;}
+   void SetDumpingParam(const double dump_param) {theta = dump_param;}
+   virtual void SetOperator(const Operator &op) {}
+   virtual void Mult(const Vector &r, Vector &z) const;
+   virtual ~SchwarzSmoother();
+};
+
+
+class ComplexSchwarzSmoother : public Solver {
+private:
+   MPI_Comm comm;
+   int nrpatch;
+   int maxit = 1;
+   double theta = 1.0;
+   Array<int> host_rank;
+#ifdef MFEM_USE_SUITESPARSE
+   Array<UMFPackSolver * > PatchInv;
+#else   
+   Array<GMRESSolver * > PatchInv;
+#endif   /// The linear system matrix
+   ComplexHypreParMatrix * A;
+   PatchAssembly * P_r;
+   PatchAssembly * P_i;
+   PatchRestriction * R_r= nullptr;
+   PatchRestriction * R_i= nullptr;
+public:
+   ComplexSchwarzSmoother(ParMesh * cpmesh_, int ref_levels_,ParFiniteElementSpace *fespace_, ComplexHypreParMatrix * A_);
+   void SetNumSmoothSteps(const int iter) {maxit = iter;}
+   void SetDumpingParam(const double dump_param) {theta = dump_param;}
+   virtual void SetOperator(const Operator &op) {}
+   virtual void Mult(const Vector &r, Vector &z) const;
+   virtual ~ComplexSchwarzSmoother();
+};
+
+

examples/ex25p.cpp

@@ -453,7 +453,7 @@ int main(int argc, char *argv[])
    //     assembly, eliminating boundary conditions, applying conforming
    //     constraints for non-conforming AMR, etc.
    if (pa) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
-   a.Assemble();
+   a.Assemble(0);
 
    OperatorPtr Ah;
    Vector B, X;
@@ -478,13 +478,30 @@ int main(int argc, char *argv[])
 #ifdef MFEM_USE_MUMPS
    if (!pa && mumps_solver)
    {
+
+      Vector Y(X);
       HypreParMatrix *A = Ah.As<ComplexHypreParMatrix>()->GetSystemMatrix();
       MUMPSSolver mumps;
       mumps.SetPrintLevel(0);
       mumps.SetMatrixSymType(MUMPSSolver::MatType::UNSYMMETRIC);
       mumps.SetOperator(*A);
-      mumps.Mult(B,X);
+
+      mumps.Mult(B,Y);
       delete A;
+
+      // complex mumps
+      ComplexHypreParMatrix * Ac = Ah.As<ComplexHypreParMatrix>();
+      ComplexMUMPSSolver cmumps;
+      cmumps.SetPrintLevel(0);
+      cmumps.SetOperator(*Ac);
+      cmumps.Mult(B,X);
+
+      Y-=X;
+      double diff = Y.Norml2();
+      if (myid == 0)
+      {
+         cout << "||X - Y|| = " << diff << endl;
+      }
    }
 #endif
    // 16a. Set up the parallel Bilinear form a(.,.) for the preconditioner

linalg/complex_operator.cpp

@@ -875,6 +875,432 @@ ComplexHypreParMatrix::getColStartStop(const HypreParMatrix * A_r,
    delete [] stat;
 }
 
+#ifdef MFEM_USE_MUMPS
+
+void ComplexMUMPSSolver::SetOperator(const Operator &op)
+{
+   auto APtr = dynamic_cast<const ComplexHypreParMatrix *>(&op);
+
+   MFEM_VERIFY(APtr, "Not compatible matrix type");
+   height = op.Height();
+   width = op.Width();
+
+   conv = APtr->GetConvention();
+   comm = APtr->real().GetComm();
+   MPI_Comm_size(comm, &numProcs);
+   MPI_Comm_rank(comm, &myid);
+
+   auto parcsr_op_r = (hypre_ParCSRMatrix *) const_cast<HypreParMatrix &>
+                      (APtr->real());
+   auto parcsr_op_i = (hypre_ParCSRMatrix *) const_cast<HypreParMatrix &>
+                      (APtr->imag());
+
+   hypre_CSRMatrix *csr_op_r = hypre_MergeDiagAndOffd(parcsr_op_r);
+   hypre_CSRMatrix *csr_op_i = hypre_MergeDiagAndOffd(parcsr_op_i);
+#if MFEM_HYPRE_VERSION >= 21600
+   hypre_CSRMatrixBigJtoJ(csr_op_r);
+   hypre_CSRMatrixBigJtoJ(csr_op_i);
+#endif
+   MFEM_VERIFY(csr_op_r->num_nonzeros == csr_op_i->num_nonzeros,
+               "Incompatible sparsity partters");
+
+   int *Iptr = csr_op_r->i;
+   int *Jptr = csr_op_r->j;
+   int n_loc = csr_op_r->num_rows;
+   row_start = parcsr_op_i->first_row_index;
+   MUMPS_INT8 nnz = csr_op_r->num_nonzeros;
+
+   int * I = new int[nnz];
+   int * J = new int[nnz];
+
+   // Fill in I and J arrays for
+   // COO format in 1-based indexing
+   int k = 0;
+   double * data_r = csr_op_r->data;
+   double * data_i = csr_op_i->data;
+   mumps_double_complex *zdata = new mumps_double_complex[nnz];
+   for (int i = 0; i < n_loc; i++)
+   {
+      for (int j = Iptr[i]; j < Iptr[i + 1]; j++)
+      {
+         I[k] = row_start + i + 1;
+         J[k] = Jptr[k] + 1;
+         zdata[k].r = data_r[k];
+         zdata[k].i = data_i[k];
+         k++;
+      }
+   }
+
+   // new MUMPS object
+   if (id)
+   {
+      id->job = -2;
+      zmumps_c(id);
+      delete id;
+   }
+   id = new ZMUMPS_STRUC_C;
+   // C to Fortran communicator
+   id->comm_fortran = (MUMPS_INT) MPI_Comm_c2f(comm);
+
+   // Host is involved in computation
+   id->par = 1;
+
+   id->sym = 0;
+
+   // MUMPS init
+   id->job = -1;
+   zmumps_c(id);
+
+   // Set MUMPS default parameters
+   SetParameters();
+
+   id->n = parcsr_op_r->global_num_rows;
+
+   id->nnz_loc = nnz;
+
+   id->irn_loc = I;
+
+   id->jcn_loc = J;
+
+   id->a_loc = zdata;
+
+   // MUMPS Analysis
+   id->job = 1;
+   zmumps_c(id);
+
+   // MUMPS Factorization
+   id->job = 2;
+   zmumps_c(id);
+
+   hypre_CSRMatrixDestroy(csr_op_r);
+   hypre_CSRMatrixDestroy(csr_op_i);
+   delete [] I;
+   delete [] J;
+   delete [] zdata;
+
+#if MFEM_MUMPS_VERSION >= 530
+   delete [] irhs_loc;
+   irhs_loc = new int[n_loc];
+   for (int i = 0; i < n_loc; i++)
+   {
+      irhs_loc[i] = row_start + i + 1;
+   }
+   row_starts.SetSize(numProcs);
+   MPI_Allgather(&row_start, 1, MPI_INT, row_starts, 1, MPI_INT, comm);
+#else
+   if (myid == 0)
+   {
+      delete [] rhs_glob;
+      delete [] recv_counts;
+      global_num_rows = parcsr_op_r->global_num_rows;
+      rhs_glob = new mumps_double_complex[global_num_rows];
+      recv_counts = new int[numProcs];
+   }
+   MPI_Gather(&n_loc, 1, MPI_INT, recv_counts, 1, MPI_INT, 0, comm);
+   if (myid == 0)
+   {
+      delete [] displs;
+      displs = new int[numProcs];
+      displs[0] = 0;
+      int s = 0;
+      for (int k = 0; k < numProcs-1; k++)
+      {
+         s += recv_counts[k];
+         displs[k+1] = s;
+      }
+   }
+#endif
+}
+
+void ComplexMUMPSSolver::Mult(const Vector &x, Vector &y) const
+{
+   int n = x.Size()/2;
+   double * datax = x.GetData();
+   Vector ximag;
+   if (conv == ComplexOperator::Convention::BLOCK_SYMMETRIC)
+   {
+      ximag.SetDataAndSize(&datax[n],n);
+      ximag *=-1.0;
+   }
+#if MFEM_MUMPS_VERSION >= 530
+   id->nloc_rhs = n;
+   id->lrhs_loc = n;
+   mumps_double_complex *zx = new mumps_double_complex[n];
+   for (int i = 0; i<n; i++)
+   {
+      zx[i].r = x[i];
+      zx[i].i = x[n+i];
+   }
+   id->rhs_loc = zx;
+   id->irhs_loc = irhs_loc;
+
+   id->lsol_loc = id->MUMPSC_INFO(23);
+   id->isol_loc = new int[id->MUMPSC_INFO(23)];
+   id->sol_loc = new mumps_double_complex[id->MUMPSC_INFO(23)];
+
+   // MUMPS solve
+   id->job = 3;
+   zmumps_c(id);
+
+   double *zy = new double[2*id->MUMPSC_INFO(23)];
+   for (int i = 0; i<id->MUMPSC_INFO(23); i++)
+   {
+      zy[i] = id->sol_loc[i].r;
+      zy[id->MUMPSC_INFO(23)+i] = id->sol_loc[i].i;
+   }
+
+   RedistributeSol(id->isol_loc, zy, y.GetData());
+
+   delete [] zy;
+   delete [] zx;
+   delete [] id->sol_loc;
+   delete [] id->isol_loc;
+#else
+   // real
+   double * rhs_glob_r = nullptr;
+   double * rhs_glob_i = nullptr;
+   if (myid == 0)
+   {
+      rhs_glob_r = new double[global_num_rows];
+      rhs_glob_i = new double[global_num_rows];
+   }
+   double * xdata = x.GetData();
+   MPI_Gatherv(xdata, n, MPI_DOUBLE,
+               rhs_glob_r, recv_counts,
+               displs, MPI_DOUBLE, 0, comm);
+   MPI_Gatherv(&xdata[n], n, MPI_DOUBLE,
+               rhs_glob_i, recv_counts,
+               displs, MPI_DOUBLE, 0, comm);
+
+   if (myid == 0)
+   {
+      for (int i = 0; i<global_num_rows; i++)
+      {
+         rhs_glob[i].r = rhs_glob_r[i];
+         rhs_glob[i].i = rhs_glob_i[i];
+      }
+      id->rhs = rhs_glob;
+   }
+   // MUMPS solve
+   id->job = 3;
+   zmumps_c(id);
+   if (myid == 0)
+   {
+      for (int i = 0; i<global_num_rows; i++)
+      {
+         rhs_glob_r[i] = rhs_glob[i].r;
+         rhs_glob_i[i] = rhs_glob[i].i;
+      }
+   }
+   double * ydata = y.GetData();
+   MPI_Scatterv(rhs_glob_r, recv_counts, displs,
+                MPI_DOUBLE, ydata, n,
+                MPI_DOUBLE, 0, comm);
+   MPI_Scatterv(rhs_glob_i, recv_counts, displs,
+                MPI_DOUBLE, &ydata[n], n,
+                MPI_DOUBLE, 0, comm);
+
+   if (myid == 0)
+   {
+      delete [] rhs_glob_r;
+      delete [] rhs_glob_i;
+   }
+
+#endif
+
+   if (conv == ComplexOperator::Convention::BLOCK_SYMMETRIC)
+   {
+      ximag *=-1.0;
+   }
+}
+
+void ComplexMUMPSSolver::SetPrintLevel(int print_lvl)
+{
+   print_level = print_lvl;
+}
+
+ComplexMUMPSSolver::~ComplexMUMPSSolver()
+{
+   if (id)
+   {
+#if MFEM_MUMPS_VERSION >= 530
+      delete [] irhs_loc;
+#else
+      delete [] recv_counts;
+      delete [] displs;
+      delete [] rhs_glob;
+#endif
+      id->job = -2;
+      zmumps_c(id);
+      delete id;
+   }
+}
+
+void ComplexMUMPSSolver::SetParameters()
+{
+   // output stream for error messages
+   id->ICNTL(1) = 6;
+   // output stream for diagnosting printing local to each proc
+   id->ICNTL(2) = 6;
+   // output stream for global info
+   id->ICNTL(3) = 6;
+   // Level of error printing
+   id->ICNTL(4) = print_level;
+   //input matrix format (assembled)
+   id->ICNTL(5) = 0;
+   // Use A or A^T
+   id->ICNTL(9) = 1;
+   // Iterative refinement (disabled)
+   id->ICNTL(10) = 0;
+   // Error analysis-statistics (disabled)
+   id->ICNTL(11) = 0;
+   // Use of ScaLAPACK (Parallel factorization on root)
+   id->ICNTL(13) = 0;
+   // Percentage increase of estimated workspace (default = 20%)
+   id->ICNTL(14) = 20;
+   // Number of OpenMP threads (default)
+   id->ICNTL(16) = 0;
+   // Matrix input format (distributed)
+   id->ICNTL(18) = 3;
+   // Schur complement (no Schur complement matrix returned)
+   id->ICNTL(19) = 0;
+
+#if MFEM_MUMPS_VERSION >= 530
+   // Distributed RHS
+   id->ICNTL(20) = 10;
+   // Distributed Sol
+   id->ICNTL(21) = 1;
+#else
+   // Centralized RHS
+   id->ICNTL(20) = 0;
+   // Centralized Sol
+   id->ICNTL(21) = 0;
+#endif
+   // Out of core factorization and solve (disabled)
+   id->ICNTL(22) = 0;
+   // Max size of working memory (default = based on estimates)
+   id->ICNTL(23) = 0;
+}
+
+#if MFEM_MUMPS_VERSION >= 530
+int ComplexMUMPSSolver::GetRowRank(int i, const Array<int> &row_starts_) const
+{
+   if (row_starts_.Size() == 1)
+   {
+      return 0;
+   }
+   auto up = std::upper_bound(row_starts_.begin(), row_starts_.end(), i);
+   return std::distance(row_starts_.begin(), up) - 1;
+}
+
+void ComplexMUMPSSolver::RedistributeSol(const int * row_map,
+                                         const double * x, double * y) const
+{
+   int size = id->MUMPSC_INFO(23);
+   int n = id->nloc_rhs;
+   int * send_count = new int[numProcs]();
+   for (int i = 0; i < size; i++)
+   {
+      int j = row_map[i] - 1;
+      int row_rank = GetRowRank(j, row_starts);
+      if (myid == row_rank) { continue; }
+      send_count[row_rank]++;
+   }
+
+   int * recv_count = new int[numProcs];
+   MPI_Alltoall(send_count, 1, MPI_INT, recv_count, 1, MPI_INT, comm);
+
+   int * send_displ = new int [numProcs]; send_displ[0] = 0;
+   int * recv_displ = new int [numProcs]; recv_displ[0] = 0;
+   int sbuff_size = send_count[numProcs-1];
+   int rbuff_size = recv_count[numProcs-1];
+   for (int k = 0; k < numProcs - 1; k++)
+   {
+      send_displ[k + 1] = send_displ[k] + send_count[k];
+      recv_displ[k + 1] = recv_displ[k] + recv_count[k];
+      sbuff_size += send_count[k];
+      rbuff_size += recv_count[k];
+   }
+
+   int * sendbuf_index = new int[sbuff_size];
+   double * sendbuf_values_r = new double[sbuff_size];
+   double * sendbuf_values_i = new double[sbuff_size];
+   int * soffs = new int[numProcs]();
+
+   for (int i = 0; i < size; i++)
+   {
+      int j = row_map[i] - 1;
+      int row_rank = GetRowRank(j, row_starts);
+      if (myid == row_rank)
+      {
+         int local_index = j - row_start;
+         y[local_index] = x[i];
+         y[local_index+n] = x[i+size];
+      }
+      else
+      {
+         int k = send_displ[row_rank] + soffs[row_rank];
+         sendbuf_index[k] = j;
+         sendbuf_values_r[k] = x[i];
+         sendbuf_values_i[k] = x[i+size];
+         soffs[row_rank]++;
+      }
+   }
+
+   int * recvbuf_index = new int[rbuff_size];
+   double * recvbuf_values_r = new double[rbuff_size];
+   double * recvbuf_values_i = new double[rbuff_size];
+   MPI_Alltoallv(sendbuf_index,
+                 send_count,
+                 send_displ,
+                 MPI_INT,
+                 recvbuf_index,
+                 recv_count,
+                 recv_displ,
+                 MPI_INT,
+                 comm);
+   MPI_Alltoallv(sendbuf_values_r,
+                 send_count,
+                 send_displ,
+                 MPI_DOUBLE,
+                 recvbuf_values_r,
+                 recv_count,
+                 recv_displ,
+                 MPI_DOUBLE,
+                 comm);
+   MPI_Alltoallv(sendbuf_values_i,
+                 send_count,
+                 send_displ,
+                 MPI_DOUBLE,
+                 recvbuf_values_i,
+                 recv_count,
+                 recv_displ,
+                 MPI_DOUBLE,
+                 comm);
+
+   // Unpack recv buffer
+   for (int i = 0; i < rbuff_size; i++)
+   {
+      int local_index = recvbuf_index[i] - row_start;
+      y[local_index] = recvbuf_values_r[i];
+      y[local_index+n] = recvbuf_values_i[i];
+   }
+
+   delete [] recvbuf_values_r;
+   delete [] recvbuf_values_i;
+   delete [] recvbuf_index;
+   delete [] soffs;
+   delete [] sendbuf_values_r;
+   delete [] sendbuf_values_i;
+   delete [] sendbuf_index;
+   delete [] recv_displ;
+   delete [] send_displ;
+   delete [] recv_count;
+   delete [] send_count;
+}
+#endif // MUMPS VERSION
+#endif // MFEM_USE_CMUMPS
+
 #endif // MFEM_USE_MPI
 
 }

linalg/complex_operator.hpp

@@ -22,6 +22,11 @@
 #include <umfpack.h>
 #endif
 
+#ifdef MFEM_USE_MUMPS
+#include "zmumps_c.h"
+#include <vector>
+#endif
+
 namespace mfem
 {
 
@@ -286,6 +291,49 @@ private:
    int nranks_;
 };
 
+#ifdef MFEM_USE_MUMPS
+
+class ComplexMUMPSSolver : public mfem::Solver
+{
+public:
+   ComplexMUMPSSolver() {}
+   void SetOperator(const Operator &op);
+   void Mult(const Vector &x, Vector &y) const;
+   void SetPrintLevel(int print_lvl);
+   ~ComplexMUMPSSolver();
+private:
+   MPI_Comm comm;
+   ComplexOperator::Convention conv;
+   int numProcs;
+   int myid;
+   int print_level = 0;
+   int row_start;
+#define ICNTL(I) icntl[(I) -1]
+#define MUMPSC_INFO(I) info[(I) -1]
+   ZMUMPS_STRUC_C *id=nullptr;
+   void SetParameters();
+
+#if MFEM_MUMPS_VERSION >= 530
+
+   Array<int> row_starts;
+   int * irhs_loc = nullptr;
+
+   int GetRowRank(int i, const Array<int> &row_starts_) const;
+
+   void RedistributeSol(const int * row_map,
+                        const double * x,
+                        double * y) const;
+#else
+   int global_num_rows;
+   int * recv_counts = nullptr;
+   int * displs = nullptr;
+   mumps_double_complex * rhs_glob = nullptr;
+#endif
+
+}; // mfem::ComplexMUMPSSolver class
+
+#endif // MFEM_USE_CMUMPS
+
 #endif // MFEM_USE_MPI
 
 }