complex-jacobi minor cleanup

Co-authored-by: Stefan Henneking <stefan.henneking@gmail.com>

.gitignore

@@ -50,6 +50,7 @@ examples/ex1[04-9]
 examples/ex1[0-9]p
 examples/ex2[0-9]
 examples/ex2[0-9]p
+examples/ex25-gpu
 
 examples/refined.mesh
 examples/displaced.mesh

CHANGELOG

@@ -60,6 +60,10 @@ Improved GPU capabilities
 
 - Added support for BlockOperator on GPU. See the updated Example 5.
 
+- Added partial assembly and GPU support for ComplexOperator,
+  [Par]ComplexGridFunction, [Par]ComplexLinearForm, and [Par]SesquilinearForm.
+  See the updated Example 22.
+
 Discretization improvements
 ---------------------------
 - Added support for matrix-free interpolation and restriction operators between
@@ -149,9 +153,6 @@ New and updated examples and miniapps
 - Added a new meshing miniapp, Minimal Surface, which solves Plateau's problem:
   the Dirichlet problem for the minimal surface equation.
 
-- Added partial assembly support to Example 4/4p and Example 5/5p, with diagonal
-  preconditioning.
-
 - Added full assembly support in Example 9/9p.
 
 - Added a new test problem in Example 24/24p, demonstrating a mixed bilinear
@@ -163,7 +164,8 @@ New and updated examples and miniapps
   mesh based on element attributes. Any newly exposed boundary elements are
   assigned attribute numbers related to the trimmed element attributes.
 
-- Added device support in Example 5/5p.
+- Added partial assembly and device support to Example 4/4p, Example 5/5p,
+  Example 22/22p, and Example 25/25p, with diagonal preconditioning.
 
 Improved testing
 ----------------

examples/ex22.cpp

@@ -13,10 +13,10 @@
 //               ex22 -m ../data/inline-hex.mesh -o 2 -p 2
 //               ex22 -m ../data/star.mesh -r 1 -o 2 -sigma 10.0
 //
-// With partial assembly:
-//               ex22 -m ../data/inline-quad.mesh -o 3 -p 1 -pa
-//               ex22 -m ../data/inline-hex.mesh -o 2 -p 2 -pa
-//               ex22 -m ../data/star.mesh -r 1 -o 2 -sigma 10.0 -pa
+// Device sample runs:
+//               ex22 -m ../data/inline-quad.mesh -o 3 -p 1 -pa -d cuda
+//               ex22 -m ../data/inline-hex.mesh -o 2 -p 2 -pa -d cuda
+//               ex22 -m ../data/star.mesh -r 1 -o 2 -sigma 10.0 -pa -d cuda
 //
 // Description:  This example code demonstrates the use of MFEM to define and
 //               solve simple complex-valued linear systems. It implements three
@@ -82,6 +82,7 @@ int main(int argc, char *argv[])
    bool herm_conv = true;
    bool exact_sol = true;
    bool pa = false;
+   const char *device_config = "cpu";
 
    OptionsParser args(argc, argv);
    args.AddOption(&mesh_file, "-m", "--mesh",
@@ -114,6 +115,8 @@ int main(int argc, char *argv[])
                   "Enable or disable GLVis visualization.");
    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.Parse();
    if (!args.Good())
    {
@@ -143,13 +146,18 @@ int main(int argc, char *argv[])
    ComplexOperator::Convention conv =
       herm_conv ? ComplexOperator::HERMITIAN : ComplexOperator::BLOCK_SYMMETRIC;
 
-   // 2. Read the mesh from the given mesh file. We can handle triangular,
+   // 2. 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);
+   device.Print();
+
+   // 3. Read the mesh from the given mesh file. We can handle triangular,
    //    quadrilateral, tetrahedral, hexahedral, surface and volume meshes
    //    with the same code.
    Mesh *mesh = new Mesh(mesh_file, 1, 1);
    int dim = mesh->Dimension();
 
-   // 3. Refine the mesh to increase resolution. In this example we do
+   // 4. Refine the mesh to increase resolution. In this example we do
    //    'ref_levels' of uniform refinement where the user specifies
    //    the number of levels with the '-r' option.
    for (int l = 0; l < ref_levels; l++)
@@ -157,7 +165,7 @@ int main(int argc, char *argv[])
       mesh->UniformRefinement();
    }
 
-   // 4. Define a finite element space on the mesh. Here we use continuous
+   // 5. Define a finite element space on the mesh. Here we use continuous
    //    Lagrange, Nedelec, or Raviart-Thomas finite elements of the specified
    //    order.
    if (dim == 1 && prob != 0 )
@@ -179,7 +187,7 @@ int main(int argc, char *argv[])
    cout << "Number of finite element unknowns: " << fespace->GetTrueVSize()
         << endl;
 
-   // 5. Determine the list of true (i.e. conforming) essential boundary dofs.
+   // 6. Determine the list of true (i.e. conforming) essential boundary dofs.
    //    In this example, the boundary conditions are defined based on the type
    //    of mesh and the problem type.
    Array<int> ess_tdof_list;
@@ -191,12 +199,12 @@ int main(int argc, char *argv[])
       fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
    }
 
-   // 6. Set up the linear form b(.) which corresponds to the right-hand side of
+   // 7. Set up the linear form b(.) which corresponds to the right-hand side of
    //    the FEM linear system.
    ComplexLinearForm b(fespace, conv);
    b.Vector::operator=(0.0);
 
-   // 7. Define the solution vector u as a complex finite element grid function
+   // 8. Define the solution vector u as a complex finite element grid function
    //    corresponding to fespace. Initialize u with initial guess of 1+0i or
    //    the exact solution if it is known.
    ComplexGridFunction u(fespace);
@@ -218,7 +226,6 @@ int main(int argc, char *argv[])
    VectorConstantCoefficient zeroVecCoef(zeroVec);
    VectorConstantCoefficient oneVecCoef(oneVec);
 
-   u = 0.0;
    switch (prob)
    {
       case 0:
@@ -271,7 +278,7 @@ int main(int argc, char *argv[])
                  << "window_title 'Exact: Imaginary Part'" << flush;
    }
 
-   // 8. Set up the sesquilinear form a(.,.) on the finite element space
+   // 9. Set up the sesquilinear form a(.,.) on the finite element space
    //    corresponding to the damped harmonic oscillator operator of the
    //    appropriate type:
    //
@@ -314,7 +321,7 @@ int main(int argc, char *argv[])
       default: break; // This should be unreachable
    }
 
-   // 8a. Set up the bilinear form for the preconditioner corresponding to the
+   // 9a. Set up the bilinear form for the preconditioner corresponding to the
    //     appropriate operator
    //
    //      0) A scalar H1 field
@@ -349,9 +356,9 @@ int main(int argc, char *argv[])
       default: break; // This should be unreachable
    }
 
-   // 9. Assemble the form and the corresponding linear system, applying any
-   //    necessary transformations such as: assembly, eliminating boundary
-   //    conditions, conforming constraints for non-conforming AMR, etc.
+   // 10. Assemble the form and the corresponding linear system, applying any
+   //     necessary transformations such as: assembly, eliminating boundary
+   //     conditions, conforming constraints for non-conforming AMR, etc.
    a->Assemble();
    pcOp->Assemble();
 
@@ -362,7 +369,7 @@ int main(int argc, char *argv[])
 
    cout << "Size of linear system: " << A->Width() << endl << endl;
 
-   // 10. Define and apply a GMRES solver for AU=B with a block diagonal
+   // 11. Define and apply a GMRES solver for AU=B with a block diagonal
    //     preconditioner based on the appropriate sparse smoother.
    {
       Array<int> blockOffsets;
@@ -419,9 +426,11 @@ int main(int argc, char *argv[])
       gmres.Mult(B, U);
    }
 
-   // 11. Recover the solution as a finite element grid function and compute the
+   // 12. Recover the solution as a finite element grid function and compute the
    //     errors if the exact solution is known.
    a->RecoverFEMSolution(U, b, u);
+   u.real().SyncMemory(u);
+   u.imag().SyncMemory(u);
 
    if (exact_sol)
    {
@@ -451,7 +460,7 @@ int main(int argc, char *argv[])
       cout << endl;
    }
 
-   // 12. Save the refined mesh and the solution. This output can be viewed
+   // 13. Save the refined mesh and the solution. This output can be viewed
    //     later using GLVis: "glvis -m mesh -g sol".
    {
       ofstream mesh_ofs("refined.mesh");
@@ -466,7 +475,7 @@ int main(int argc, char *argv[])
       u.imag().Save(sol_i_ofs);
    }
 
-   // 13. Send the solution by socket to a GLVis server.
+   // 14. Send the solution by socket to a GLVis server.
    if (visualization)
    {
       char vishost[] = "localhost";
@@ -525,7 +534,7 @@ int main(int argc, char *argv[])
       }
    }
 
-   // 14. Free the used memory.
+   // 15. Free the used memory.
    delete a;
    delete u_exact;
    delete pcOp;

examples/ex22p.cpp

@@ -13,10 +13,10 @@
 //               mpirun -np 4 ex22p -m ../data/inline-hex.mesh -o 2 -p 2
 //               mpirun -np 4 ex22p -m ../data/star.mesh -o 2 -sigma 10.0
 //
-// With partial assembly:
-//               mpirun -np 4 ex22p -m ../data/inline-quad.mesh -o 1 -p 1 -pa
-//               mpirun -np 4 ex22p -m ../data/inline-hex.mesh -o 1 -p 2 -pa
-//               mpirun -np 4 ex22p -m ../data/star.mesh -o 2 -sigma 10.0 -pa
+// Device sample runs:
+//               mpirun -np 4 ex22p -m ../data/inline-quad.mesh -o 1 -p 1 -pa -d cuda
+//               mpirun -np 4 ex22p -m ../data/inline-hex.mesh -o 1 -p 2 -pa -d cuda
+//               mpirun -np 4 ex22p -m ../data/star.mesh -o 2 -sigma 10.0 -pa -d cuda
 //
 // Description:  This example code demonstrates the use of MFEM to define and
 //               solve simple complex-valued linear systems. It implements three
@@ -46,7 +46,6 @@
 //               We recommend viewing examples 1, 3 and 4 before viewing this
 //               example.
 
-
 #include "mfem.hpp"
 #include <fstream>
 #include <iostream>
@@ -90,6 +89,7 @@ int main(int argc, char *argv[])
    bool herm_conv = true;
    bool exact_sol = true;
    bool pa = false;
+   const char *device_config = "cpu";
 
    OptionsParser args(argc, argv);
    args.AddOption(&mesh_file, "-m", "--mesh",
@@ -124,6 +124,8 @@ int main(int argc, char *argv[])
                   "Enable or disable GLVis visualization.");
    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.Parse();
    if (!args.Good())
    {
@@ -160,19 +162,24 @@ int main(int argc, char *argv[])
    ComplexOperator::Convention conv =
       herm_conv ? ComplexOperator::HERMITIAN : ComplexOperator::BLOCK_SYMMETRIC;
 
-   // 3. Read the (serial) mesh from the given mesh file on all processors. We
+   // 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);
+   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);
    int dim = mesh->Dimension();
 
-   // 4. Refine the serial mesh on all processors to increase the resolution.
+   // 5. Refine the serial mesh on all processors to increase the resolution.
    for (int l = 0; l < ser_ref_levels; l++)
    {
       mesh->UniformRefinement();
    }
 
-   // 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
+   // 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);
@@ -182,7 +189,7 @@ int main(int argc, char *argv[])
       pmesh->UniformRefinement();
    }
 
-   // 6. Define a parallel finite element space on the parallel mesh. Here we
+   // 7. Define a parallel finite element space on the parallel mesh. Here we
    //    use continuous Lagrange, Nedelec, or Raviart-Thomas finite elements of
    //    the specified order.
    if (dim == 1 && prob != 0 )
@@ -210,7 +217,7 @@ int main(int argc, char *argv[])
       cout << "Number of finite element unknowns: " << size << endl;
    }
 
-   // 7. Determine the list of true (i.e. parallel conforming) essential
+   // 8. Determine the list of true (i.e. parallel conforming) essential
    //    boundary dofs. In this example, the boundary conditions are defined
    //    based on the type of mesh and the problem type.
    Array<int> ess_tdof_list;
@@ -222,14 +229,14 @@ int main(int argc, char *argv[])
       fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
    }
 
-   // 8. Set up the parallel linear form b(.) which corresponds to the
+   // 9. Set up the parallel linear form b(.) which corresponds to the
    //    right-hand side of the FEM linear system.
    ParComplexLinearForm b(fespace, conv);
    b.Vector::operator=(0.0);
 
-   // 9. Define the solution vector u as a parallel complex finite element grid
-   //    function corresponding to fespace. Initialize u with initial guess of
-   //    1+0i or the exact solution if it is known.
+   // 10. Define the solution vector u as a parallel complex finite element grid
+   //     function corresponding to fespace. Initialize u with initial guess of
+   //     1+0i or the exact solution if it is known.
    ParComplexGridFunction u(fespace);
    ParComplexGridFunction * u_exact = NULL;
    if (exact_sol) { u_exact = new ParComplexGridFunction(fespace); }
@@ -249,7 +256,6 @@ int main(int argc, char *argv[])
    VectorConstantCoefficient zeroVecCoef(zeroVec);
    VectorConstantCoefficient oneVecCoef(oneVec);
 
-   u = 0.0;
    switch (prob)
    {
       case 0:
@@ -304,7 +310,7 @@ int main(int argc, char *argv[])
                  << "window_title 'Exact: Imaginary Part'" << flush;
    }
 
-   // 10. Set up the parallel sesquilinear form a(.,.) on the finite element
+   // 11. Set up the parallel sesquilinear form a(.,.) on the finite element
    //     space corresponding to the damped harmonic oscillator operator of the
    //     appropriate type:
    //
@@ -347,7 +353,7 @@ int main(int argc, char *argv[])
       default: break; // This should be unreachable
    }
 
-   // 10a. Set up the parallel bilinear form for the preconditioner
+   // 11a. Set up the parallel bilinear form for the preconditioner
    //      corresponding to the appropriate operator
    //
    //      0) A scalar H1 field
@@ -381,7 +387,7 @@ int main(int argc, char *argv[])
       default: break; // This should be unreachable
    }
 
-   // 11. Assemble the parallel bilinear form and the corresponding linear
+   // 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, etc.
@@ -399,7 +405,7 @@ int main(int argc, char *argv[])
            << 2 * fespace->GlobalTrueVSize() << endl << endl;
    }
 
-   // 12. Define and apply a parallel FGMRES solver for AU=B with a block
+   // 13. Define and apply a parallel FGMRES solver for AU=B with a block
    //     diagonal preconditioner based on the appropriate multigrid
    //     preconditioner from hypre.
    {
@@ -460,9 +466,11 @@ int main(int argc, char *argv[])
       fgmres.SetPrintLevel(1);
       fgmres.Mult(B, U);
    }
-   // 13. Recover the parallel grid function corresponding to U. This is the
+   // 14. Recover the parallel grid function corresponding to U. This is the
    //     local finite element solution on each processor.
    a->RecoverFEMSolution(U, b, u);
+   u.real().SyncMemory(u);
+   u.imag().SyncMemory(u);
 
    if (exact_sol)
    {
@@ -495,7 +503,7 @@ int main(int argc, char *argv[])
       }
    }
 
-   // 14. Save the refined mesh and the solution in parallel. This output can be
+   // 15. 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_r_name, sol_i_name;
@@ -515,7 +523,7 @@ int main(int argc, char *argv[])
       u.imag().Save(sol_i_ofs);
    }
 
-   // 15. Send the solution by socket to a GLVis server.
+   // 16. Send the solution by socket to a GLVis server.
    if (visualization)
    {
       char vishost[] = "localhost";
@@ -580,7 +588,7 @@ int main(int argc, char *argv[])
       }
    }
 
-   // 16. Free the used memory.
+   // 17. Free the used memory.
    delete a;
    delete u_exact;
    delete pcOp;

examples/ex24.cpp

@@ -70,7 +70,7 @@ int main(int argc, char *argv[])
    args.AddOption(&order, "-o", "--order",
                   "Finite element order (polynomial degree).");
    args.AddOption(&prob, "-p", "--problem-type",
-                  "Choose between 0: H(Curl) or 1: H(Div)");
+                  "Choose between 0: grad, 1: curl, 2: div");
    args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
                   "--no-static-condensation", "Enable static condensation.");
    args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",

examples/ex24p.cpp

@@ -76,7 +76,7 @@ int main(int argc, char *argv[])
    args.AddOption(&order, "-o", "--order",
                   "Finite element order (polynomial degree).");
    args.AddOption(&prob, "-p", "--problem-type",
-                  "Choose between 0: H(Curl) or 1: H(Div)");
+                  "Choose between 0: grad, 1: curl, 2: div");
    args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
                   "--no-static-condensation", "Enable static condensation.");
    args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",

examples/ex25.cpp

@@ -10,6 +10,10 @@
 //               ex25 -o 2 -f 8.0 -ref 3 -prob 4 -m ../data/inline-quad.mesh
 //               ex25 -o 2 -f 2.0 -ref 1 -prob 4 -m ../data/inline-hex.mesh
 //
+// Device sample runs:
+//               ex25 -o 2 -f 8.0 -ref 3 -prob 4 -m ../data/inline-quad.mesh -pa -d cuda
+//               ex25 -o 2 -f 2.0 -ref 1 -prob 4 -m ../data/inline-hex.mesh -pa -d cuda
+//
 // Description:  This example code solves a simple electromagnetic wave
 //               propagation problem corresponding to the second order
 //               indefinite Maxwell equation
@@ -82,24 +86,24 @@ public:
 };
 
 // Class for returning the PML coefficients of the bilinear form
-class PMLMatrixCoefficient : public MatrixCoefficient
+class PMLDiagMatrixCoefficient : public VectorCoefficient
 {
 private:
    CartesianPML * pml = nullptr;
-   void (*Function)(const Vector &, CartesianPML * , DenseMatrix &);
+   void (*Function)(const Vector &, CartesianPML * , Vector &);
 public:
-   PMLMatrixCoefficient(int dim, void(*F)(const Vector &, CartesianPML *,
-                                          DenseMatrix &),
-                        CartesianPML * pml_)
-      : MatrixCoefficient(dim), pml(pml_), Function(F)
+   PMLDiagMatrixCoefficient(int dim, void(*F)(const Vector &, CartesianPML *,
+                                              Vector &),
+                            CartesianPML * pml_)
+      : VectorCoefficient(dim), pml(pml_), Function(F)
    {}
-   virtual void Eval(DenseMatrix &K, ElementTransformation &T,
+   virtual void Eval(Vector &K, ElementTransformation &T,
                      const IntegrationPoint &ip)
    {
       double x[3];
       Vector transip(x, 3);
       T.Transform(ip, transip);
-      K.SetSize(height, width);
+      K.SetSize(vdim);
       (*Function)(transip, pml, K);
    }
 };
@@ -116,13 +120,13 @@ void source(const Vector &x, Vector & f);
 
 // Functions for computing the necessary coefficients after PML stretching.
 // J is the Jacobian matrix of the stretching function
-void detJ_JT_J_inv_Re(const Vector &x, CartesianPML * pml, DenseMatrix &M);
-void detJ_JT_J_inv_Im(const Vector &x, CartesianPML * pml, DenseMatrix &M);
-void detJ_JT_J_inv_abs(const Vector &x, CartesianPML * pml, DenseMatrix &M);
+void detJ_JT_J_inv_Re(const Vector &x, CartesianPML * pml, Vector &D);
+void detJ_JT_J_inv_Im(const Vector &x, CartesianPML * pml, Vector &D);
+void detJ_JT_J_inv_abs(const Vector &x, CartesianPML * pml, Vector &D);
 
-void detJ_inv_JT_J_Re(const Vector &x, CartesianPML * pml, DenseMatrix &M);
-void detJ_inv_JT_J_Im(const Vector &x, CartesianPML * pml, DenseMatrix &M);
-void detJ_inv_JT_J_abs(const Vector &x, CartesianPML * pml, DenseMatrix &M);
+void detJ_inv_JT_J_Re(const Vector &x, CartesianPML * pml, Vector &D);
+void detJ_inv_JT_J_Im(const Vector &x, CartesianPML * pml, Vector &D);
+void detJ_inv_JT_J_abs(const Vector &x, CartesianPML * pml, Vector &D);
 
 Array2D<double> comp_domain_bdr;
 Array2D<double> domain_bdr;
@@ -153,6 +157,8 @@ int main(int argc, char *argv[])
    double freq = 5.0;
    bool herm_conv = true;
    bool visualization = 1;
+   bool pa = false;
+   const char *device_config = "cpu";
 
    OptionsParser args(argc, argv);
    args.AddOption(&mesh_file, "-m", "--mesh",
@@ -174,12 +180,21 @@ int main(int argc, char *argv[])
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                   "--no-visualization",
                   "Enable or disable GLVis visualization.");
+   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.Parse();
 
    if (iprob > 4) { iprob = 4; }
    prob = (prob_type)iprob;
 
-   // 2. Setup the mesh
+   // 2. 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);
+   device.Print();
+
+   // 3. Setup the mesh
    if (!mesh_file)
    {
       exact_known = true;
@@ -220,7 +235,7 @@ int main(int argc, char *argv[])
    // Setup PML length
    Array2D<double> length(dim, 2); length = 0.0;
 
-   // 3. Setup the Cartesian PML region.
+   // 4. Setup the Cartesian PML region.
    switch (prob)
    {
       case disc:
@@ -246,19 +261,19 @@ int main(int argc, char *argv[])
    comp_domain_bdr = pml->GetCompDomainBdr();
    domain_bdr = pml->GetDomainBdr();
 
-   // 4. Refine the mesh to increase the resolution.
+   // 5. Refine the mesh to increase the resolution.
    for (int l = 0; l < ref_levels; l++)
    {
       mesh->UniformRefinement();
    }
 
-   // 5. Reorient mesh in case of a tet mesh
+   // 6. Reorient mesh in case of a tet mesh
    mesh->ReorientTetMesh();
 
    // Set element attributes in order to distinguish elements in the PML region
    pml->SetAttributes(mesh);
 
-   // 6. Define a finite element space on the mesh. Here we use the Nedelec
+   // 7. Define a finite element space on the mesh. Here we use the Nedelec
    //    finite elements of the specified order.
    FiniteElementCollection *fec = new ND_FECollection(order, dim);
    FiniteElementSpace *fespace = new FiniteElementSpace(mesh, fec);
@@ -266,7 +281,7 @@ int main(int argc, char *argv[])
 
    cout << "Number of finite element unknowns: " << size << endl;
 
-   // 7. Determine the list of true essential boundary dofs. In this example,
+   // 8. Determine the list of true essential boundary dofs. In this example,
    //    the boundary conditions are defined based on the specific mesh and the
    //    problem type.
    Array<int> ess_tdof_list;
@@ -308,12 +323,12 @@ int main(int argc, char *argv[])
    }
    fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
 
-   // 8. Setup Complex Operator convention
+   // 9. Setup Complex Operator convention
    ComplexOperator::Convention conv =
       herm_conv ? ComplexOperator::HERMITIAN : ComplexOperator::BLOCK_SYMMETRIC;
 
-   // 9. Set up the linear form b(.) which corresponds to the right-hand side of
-   //    the FEM linear system.
+   // 10. Set up the linear form b(.) which corresponds to the right-hand side of
+   //     the FEM linear system.
    VectorFunctionCoefficient f(dim, source);
    ComplexLinearForm b(fespace, conv);
    if (prob == load_src)
@@ -323,7 +338,7 @@ int main(int argc, char *argv[])
    b.Vector::operator=(0.0);
    b.Assemble();
 
-   // 10. Define the solution vector x as a complex finite element grid function
+   // 11. Define the solution vector x as a complex finite element grid function
    //     corresponding to fespace.
    ComplexGridFunction x(fespace);
    x = 0.0;
@@ -331,7 +346,7 @@ int main(int argc, char *argv[])
    VectorFunctionCoefficient E_Im(dim, E_bdr_data_Im);
    x.ProjectBdrCoefficientTangent(E_Re, E_Im, ess_bdr);
 
-   // 11. Set up the sesquilinear form a(.,.)
+   // 12. Set up the sesquilinear form a(.,.)
    //
    //     In Comp
    //     Domain:   1/mu (Curl E, Curl F) - omega^2 * epsilon (E,F)
@@ -365,19 +380,19 @@ int main(int argc, char *argv[])
    a.AddDomainIntegrator(new VectorFEMassIntegrator(restr_omeg),NULL);
 
    int cdim = (dim == 2) ? 1 : dim;
-   PMLMatrixCoefficient pml_c1_Re(cdim,detJ_inv_JT_J_Re, pml);
-   PMLMatrixCoefficient pml_c1_Im(cdim,detJ_inv_JT_J_Im, pml);
-   ScalarMatrixProductCoefficient c1_Re(muinv,pml_c1_Re);
-   ScalarMatrixProductCoefficient c1_Im(muinv,pml_c1_Im);
-   MatrixRestrictedCoefficient restr_c1_Re(c1_Re,attrPML);
-   MatrixRestrictedCoefficient restr_c1_Im(c1_Im,attrPML);
+   PMLDiagMatrixCoefficient pml_c1_Re(cdim,detJ_inv_JT_J_Re, pml);
+   PMLDiagMatrixCoefficient pml_c1_Im(cdim,detJ_inv_JT_J_Im, pml);
+   ScalarVectorProductCoefficient c1_Re(muinv,pml_c1_Re);
+   ScalarVectorProductCoefficient c1_Im(muinv,pml_c1_Im);
+   VectorRestrictedCoefficient restr_c1_Re(c1_Re,attrPML);
+   VectorRestrictedCoefficient restr_c1_Im(c1_Im,attrPML);
 
-   PMLMatrixCoefficient pml_c2_Re(dim, detJ_JT_J_inv_Re,pml);
-   PMLMatrixCoefficient pml_c2_Im(dim, detJ_JT_J_inv_Im,pml);
-   ScalarMatrixProductCoefficient c2_Re(omeg,pml_c2_Re);
-   ScalarMatrixProductCoefficient c2_Im(omeg,pml_c2_Im);
-   MatrixRestrictedCoefficient restr_c2_Re(c2_Re,attrPML);
-   MatrixRestrictedCoefficient restr_c2_Im(c2_Im,attrPML);
+   PMLDiagMatrixCoefficient pml_c2_Re(dim, detJ_JT_J_inv_Re,pml);
+   PMLDiagMatrixCoefficient pml_c2_Im(dim, detJ_JT_J_inv_Im,pml);
+   ScalarVectorProductCoefficient c2_Re(omeg,pml_c2_Re);
+   ScalarVectorProductCoefficient c2_Im(omeg,pml_c2_Im);
+   VectorRestrictedCoefficient restr_c2_Re(c2_Re,attrPML);
+   VectorRestrictedCoefficient restr_c2_Im(c2_Im,attrPML);
 
    // Integrators inside the PML region
    a.AddDomainIntegrator(new CurlCurlIntegrator(restr_c1_Re),
@@ -385,26 +400,30 @@ int main(int argc, char *argv[])
    a.AddDomainIntegrator(new VectorFEMassIntegrator(restr_c2_Re),
                          new VectorFEMassIntegrator(restr_c2_Im));
 
-   // 12. Assemble the bilinear form and the corresponding linear system,
+   // 13. Assemble the bilinear form and the corresponding linear system,
    //     applying any necessary transformations such as: assembly, eliminating
    //     boundary conditions, applying conforming constraints for
    //     non-conforming AMR, etc.
+#ifndef MFEM_USE_SUITESPARSE
+   if (pa) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
+#endif
    a.Assemble(0);
 
    OperatorPtr A;
    Vector B, X;
    a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
 
-   // 13. Solve using a direct or an iterative solver
+   // 14. Solve using a direct or an iterative solver
 #ifdef MFEM_USE_SUITESPARSE
    {
+      if (pa) { cout << "PA not available with MFEM_USE_SUITESPARSE" << endl; }
       ComplexUMFPackSolver csolver(*A.As<ComplexSparseMatrix>());
       csolver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
       csolver.SetPrintLevel(1);
       csolver.Mult(B, X);
    }
 #else
-   // 13a. Set up the Bilinear form a(.,.) for the preconditioner
+   // 14a. Set up the Bilinear form a(.,.) for the preconditioner
    //
    //    In Comp
    //    Domain:   1/mu (Curl E, Curl F) + omega^2 * epsilon (E,F)
@@ -419,50 +438,69 @@ int main(int argc, char *argv[])
       prec.AddDomainIntegrator(new CurlCurlIntegrator(restr_muinv));
       prec.AddDomainIntegrator(new VectorFEMassIntegrator(restr_absomeg));
 
-      PMLMatrixCoefficient pml_c1_abs(cdim,detJ_inv_JT_J_abs, pml);
-      ScalarMatrixProductCoefficient c1_abs(muinv,pml_c1_abs);
-      MatrixRestrictedCoefficient restr_c1_abs(c1_abs,attrPML);
+      PMLDiagMatrixCoefficient pml_c1_abs(cdim,detJ_inv_JT_J_abs, pml);
+      ScalarVectorProductCoefficient c1_abs(muinv,pml_c1_abs);
+      VectorRestrictedCoefficient restr_c1_abs(c1_abs,attrPML);
 
-      PMLMatrixCoefficient pml_c2_abs(dim, detJ_JT_J_inv_abs,pml);
-      ScalarMatrixProductCoefficient c2_abs(absomeg,pml_c2_abs);
-      MatrixRestrictedCoefficient restr_c2_abs(c2_abs,attrPML);
+      PMLDiagMatrixCoefficient pml_c2_abs(dim, detJ_JT_J_inv_abs,pml);
+      ScalarVectorProductCoefficient c2_abs(absomeg,pml_c2_abs);
+      VectorRestrictedCoefficient restr_c2_abs(c2_abs,attrPML);
 
       prec.AddDomainIntegrator(new CurlCurlIntegrator(restr_c1_abs));
       prec.AddDomainIntegrator(new VectorFEMassIntegrator(restr_c2_abs));
 
+      if (pa) { prec.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
       prec.Assemble();
 
-      OperatorPtr PCOpAh;
-      prec.FormSystemMatrix(ess_tdof_list, PCOpAh);
-
-      // 13b. Define and apply a GMRES solver for AU=B with a block diagonal
-      //      preconditioner based on the Gauss-Seidel sparse smoother.
+      // 14b. Define and apply a GMRES solver for AU=B with a block diagonal
+      //      preconditioner based on the Gauss-Seidel or Jacobi sparse smoother.
       Array<int> offsets(3);
       offsets[0] = 0;
       offsets[1] = fespace->GetTrueVSize();
       offsets[2] = fespace->GetTrueVSize();
       offsets.PartialSum();
 
-      GSSmoother gs00(*PCOpAh.As<SparseMatrix>());
-      BlockDiagonalPreconditioner BlockGS(offsets);
-      ScaledOperator gs11(&gs00,
-                          (conv == ComplexOperator::HERMITIAN) ? -1.0 : 1.0);
-      BlockGS.SetDiagonalBlock(0,&gs00);
-      BlockGS.SetDiagonalBlock(1,&gs11);
+      Operator *pc_r = nullptr;
+      Operator *pc_i = nullptr;
+      int s = (conv == ComplexOperator::HERMITIAN) ? -1.0 : 1.0;
+      if (pa)
+      {
+         // Jacobi Smoother
+         OperatorJacobiSmoother *d00 = new OperatorJacobiSmoother(prec, ess_tdof_list);
+         ScaledOperator *d11 = new ScaledOperator(d00, s);
+         pc_r = d00;
+         pc_i = d11;
+      }
+      else
+      {
+         OperatorPtr PCOpAh;
+         prec.SetDiagonalPolicy(mfem::Operator::DIAG_ONE);
+         prec.FormSystemMatrix(ess_tdof_list, PCOpAh);
+
+         // Gauss-Seidel Smoother
+         GSSmoother *gs00 = new GSSmoother(*PCOpAh.As<SparseMatrix>());
+         ScaledOperator *gs11 = new ScaledOperator(gs00, s);
+         pc_r = gs00;
+         pc_i = gs11;
+      }
+
+      BlockDiagonalPreconditioner BlockDP(offsets);
+      BlockDP.SetDiagonalBlock(0, pc_r);
+      BlockDP.SetDiagonalBlock(1, pc_i);
 
       GMRESSolver gmres;
       gmres.SetPrintLevel(1);
       gmres.SetKDim(200);
-      gmres.SetMaxIter(2000);
+      gmres.SetMaxIter(pa ? 5000 : 2000);
       gmres.SetRelTol(1e-5);
       gmres.SetAbsTol(0.0);
       gmres.SetOperator(*A);
-      gmres.SetPreconditioner(BlockGS);
+      gmres.SetPreconditioner(BlockDP);
       gmres.Mult(B, X);
    }
 #endif
 
-   // 14. Recover the solution as a finite element grid function and compute the
+   // 15. Recover the solution as a finite element grid function and compute the
    //     errors if the exact solution is known.
    a.RecoverFEMSolution(X, b, x);
 
@@ -499,7 +537,7 @@ int main(int argc, char *argv[])
            << sqrt(L2Error_Re*L2Error_Re + L2Error_Im*L2Error_Im) << "\n\n";
    }
 
-   // 15. Save the refined mesh and the solution. This output can be viewed
+   // 16. Save the refined mesh and the solution. This output can be viewed
    //     later using GLVis: "glvis -m mesh -g sol".
    {
       ofstream mesh_ofs("ex25.mesh");
@@ -514,7 +552,7 @@ int main(int argc, char *argv[])
       x.imag().Save(sol_i_ofs);
    }
 
-   // 16. Send the solution by socket to a GLVis server.
+   // 17. Send the solution by socket to a GLVis server.
    if (visualization)
    {
       // Define visualization keys for GLVis (see GLVis documentation)
@@ -565,7 +603,7 @@ int main(int argc, char *argv[])
       }
    }
 
-   // 17. Free the used memory.
+   // 18. Free the used memory.
    delete pml;
    delete fespace;
    delete fec;
@@ -763,7 +801,7 @@ void E_bdr_data_Im(const Vector &x, Vector &E)
    }
 }
 
-void detJ_JT_J_inv_Re(const Vector &x, CartesianPML * pml, DenseMatrix &M)
+void detJ_JT_J_inv_Re(const Vector &x, CartesianPML * pml, Vector &D)
 {
    vector<complex<double>> dxs(dim);
    complex<double> det(1.0, 0.0);
@@ -774,14 +812,13 @@ void detJ_JT_J_inv_Re(const Vector &x, CartesianPML * pml, DenseMatrix &M)
       det *= dxs[i];
    }
 
-   M = 0.0;
    for (int i = 0; i < dim; ++i)
    {
-      M(i, i) = (det / pow(dxs[i], 2)).real();
+      D(i) = (det / pow(dxs[i], 2)).real();
    }
 }
 
-void detJ_JT_J_inv_Im(const Vector &x, CartesianPML * pml, DenseMatrix &M)
+void detJ_JT_J_inv_Im(const Vector &x, CartesianPML * pml, Vector &D)
 {
    vector<complex<double>> dxs(dim);
    complex<double> det = 1.0;
@@ -792,14 +829,13 @@ void detJ_JT_J_inv_Im(const Vector &x, CartesianPML * pml, DenseMatrix &M)
       det *= dxs[i];
    }
 
-   M = 0.0;
    for (int i = 0; i < dim; ++i)
    {
-      M(i, i) = (det / pow(dxs[i], 2)).imag();
+      D(i) = (det / pow(dxs[i], 2)).imag();
    }
 }
 
-void detJ_JT_J_inv_abs(const Vector &x, CartesianPML * pml, DenseMatrix &M)
+void detJ_JT_J_inv_abs(const Vector &x, CartesianPML * pml, Vector &D)
 {
    vector<complex<double>> dxs(dim);
    complex<double> det = 1.0;
@@ -810,14 +846,13 @@ void detJ_JT_J_inv_abs(const Vector &x, CartesianPML * pml, DenseMatrix &M)
       det *= dxs[i];
    }
 
-   M = 0.0;
    for (int i = 0; i < dim; ++i)
    {
-      M(i, i) = abs(det / pow(dxs[i], 2));
+      D(i) = abs(det / pow(dxs[i], 2));
    }
 }
 
-void detJ_inv_JT_J_Re(const Vector &x, CartesianPML * pml, DenseMatrix &M)
+void detJ_inv_JT_J_Re(const Vector &x, CartesianPML * pml, Vector &D)
 {
    vector<complex<double>> dxs(dim);
    complex<double> det(1.0, 0.0);
@@ -831,19 +866,18 @@ void detJ_inv_JT_J_Re(const Vector &x, CartesianPML * pml, DenseMatrix &M)
    // in the 2D case the coefficient is scalar 1/det(J)
    if (dim == 2)
    {
-      M = (1.0 / det).real();
+      D = (1.0 / det).real();
    }
    else
    {
-      M = 0.0;
       for (int i = 0; i < dim; ++i)
       {
-         M(i, i) = (pow(dxs[i], 2) / det).real();
+         D(i) = (pow(dxs[i], 2) / det).real();
       }
    }
 }
 
-void detJ_inv_JT_J_Im(const Vector &x, CartesianPML * pml, DenseMatrix &M)
+void detJ_inv_JT_J_Im(const Vector &x, CartesianPML * pml, Vector &D)
 {
    vector<complex<double>> dxs(dim);
    complex<double> det = 1.0;
@@ -856,19 +890,18 @@ void detJ_inv_JT_J_Im(const Vector &x, CartesianPML * pml, DenseMatrix &M)
 
    if (dim == 2)
    {
-      M = (1.0 / det).imag();
+      D = (1.0 / det).imag();
    }
    else
    {
-      M = 0.0;
       for (int i = 0; i < dim; ++i)
       {
-         M(i, i) = (pow(dxs[i], 2) / det).imag();
+         D(i) = (pow(dxs[i], 2) / det).imag();
       }
    }
 }
 
-void detJ_inv_JT_J_abs(const Vector &x, CartesianPML * pml, DenseMatrix &M)
+void detJ_inv_JT_J_abs(const Vector &x, CartesianPML * pml, Vector &D)
 {
    vector<complex<double>> dxs(dim);
    complex<double> det = 1.0;
@@ -881,14 +914,13 @@ void detJ_inv_JT_J_abs(const Vector &x, CartesianPML * pml, DenseMatrix &M)
 
    if (dim == 2)
    {
-      M = abs(1.0 / det);
+      D = abs(1.0 / det);
    }
    else
    {
-      M = 0.0;
       for (int i = 0; i < dim; ++i)
       {
-         M(i, i) = abs(pow(dxs[i], 2) / det);
+         D(i) = abs(pow(dxs[i], 2) / det);
       }
    }
 }

examples/ex25p.cpp

@@ -10,6 +10,10 @@
 //               mpirun -np 4 ex25p -o 2 -f 8.0 -rs 2 -rp 2 -prob 4 -m ../data/inline-quad.mesh
 //               mpirun -np 4 ex25p -o 2 -f 2.0 -rs 1 -rp 1 -prob 4 -m ../data/inline-hex.mesh
 //
+// Device sample runs:
+//               mpirun -np 4 ex25p -o 1 -f 3.0 -rs 3 -rp 1 -prob 2 -pa -d cuda
+//               mpirun -np 4 ex25p -o 2 -f 1.0 -rs 1 -rp 1 -prob 3 -pa -d cuda
+//
 // Description:  This example code solves a simple electromagnetic wave
 //               propagation problem corresponding to the second order
 //               indefinite Maxwell equation
@@ -82,24 +86,24 @@ public:
 };
 
 // Class for returning the PML coefficients of the bilinear form
-class PMLMatrixCoefficient : public MatrixCoefficient
+class PMLDiagMatrixCoefficient : public VectorCoefficient
 {
 private:
    CartesianPML * pml = nullptr;
-   void (*Function)(const Vector &, CartesianPML * , DenseMatrix &);
+   void (*Function)(const Vector &, CartesianPML * , Vector &);
 public:
-   PMLMatrixCoefficient(int dim, void(*F)(const Vector &, CartesianPML *,
-                                          DenseMatrix &),
-                        CartesianPML * pml_)
-      : MatrixCoefficient(dim), pml(pml_), Function(F)
+   PMLDiagMatrixCoefficient(int dim, void(*F)(const Vector &, CartesianPML *,
+                                              Vector &),
+                            CartesianPML * pml_)
+      : VectorCoefficient(dim), pml(pml_), Function(F)
    {}
-   virtual void Eval(DenseMatrix &K, ElementTransformation &T,
+   virtual void Eval(Vector &K, ElementTransformation &T,
                      const IntegrationPoint &ip)
    {
       double x[3];
       Vector transip(x, 3);
       T.Transform(ip, transip);
-      K.SetSize(height, width);
+      K.SetSize(vdim);
       (*Function)(transip, pml, K);
    }
 };
@@ -116,13 +120,13 @@ void source(const Vector &x, Vector & f);
 
 // Functions for computing the necessary coefficients after PML stretching.
 // J is the Jacobian matrix of the stretching function
-void detJ_JT_J_inv_Re(const Vector &x, CartesianPML * pml, DenseMatrix &M);
-void detJ_JT_J_inv_Im(const Vector &x, CartesianPML * pml, DenseMatrix &M);
-void detJ_JT_J_inv_abs(const Vector &x, CartesianPML * pml, DenseMatrix &M);
+void detJ_JT_J_inv_Re(const Vector &x, CartesianPML * pml, Vector & D);
+void detJ_JT_J_inv_Im(const Vector &x, CartesianPML * pml, Vector & D);
+void detJ_JT_J_inv_abs(const Vector &x, CartesianPML * pml, Vector & D);
 
-void detJ_inv_JT_J_Re(const Vector &x, CartesianPML * pml, DenseMatrix &M);
-void detJ_inv_JT_J_Im(const Vector &x, CartesianPML * pml, DenseMatrix &M);
-void detJ_inv_JT_J_abs(const Vector &x, CartesianPML * pml, DenseMatrix &M);
+void detJ_inv_JT_J_Re(const Vector &x, CartesianPML * pml, Vector & D);
+void detJ_inv_JT_J_Im(const Vector &x, CartesianPML * pml, Vector & D);
+void detJ_inv_JT_J_abs(const Vector &x, CartesianPML * pml, Vector & D);
 
 Array2D<double> comp_domain_bdr;
 Array2D<double> domain_bdr;
@@ -160,6 +164,8 @@ int main(int argc, char *argv[])
    double freq = 5.0;
    bool herm_conv = true;
    bool visualization = 1;
+   bool pa = false;
+   const char *device_config = "cpu";
 
    OptionsParser args(argc, argv);
    args.AddOption(&mesh_file, "-m", "--mesh",
@@ -183,12 +189,21 @@ int main(int argc, char *argv[])
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                   "--no-visualization",
                   "Enable or disable GLVis visualization.");
+   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.Parse();
 
    if (iprob > 4) { iprob = 4; }
    prob = (prob_type)iprob;
 
-   // 3. Setup the (serial) mesh on all processors.
+   // 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);
+   device.Print();
+
+   // 4. Setup the (serial) mesh on all processors.
    if (!mesh_file)
    {
       exact_known = true;
@@ -236,7 +251,7 @@ int main(int argc, char *argv[])
    // Setup PML length
    Array2D<double> length(dim, 2); length = 0.0;
 
-   // 4. Setup the Cartesian PML region.
+   // 5. Setup the Cartesian PML region.
    switch (prob)
    {
       case disc:
@@ -262,13 +277,13 @@ int main(int argc, char *argv[])
    comp_domain_bdr = pml->GetCompDomainBdr();
    domain_bdr = pml->GetDomainBdr();
 
-   // 5. Refine the serial mesh on all processors to increase the resolution.
+   // 6. Refine the serial mesh on all processors to increase the resolution.
    for (int l = 0; l < ref_levels; l++)
    {
       mesh->UniformRefinement();
    }
 
-   // 6. Define a parallel mesh by a partitioning of the serial mesh.
+   // 7. Define a parallel mesh by a partitioning of the serial mesh.
    ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
    delete mesh;
    {
@@ -278,13 +293,13 @@ int main(int argc, char *argv[])
       }
    }
 
-   // 6a. Reorient mesh in case of a tet mesh
+   // 7a. Reorient mesh in case of a tet mesh
    pmesh->ReorientTetMesh();
 
-   // 7. Set element attributes in order to distinguish elements in the PML
+   // 8. Set element attributes in order to distinguish elements in the PML
    pml->SetAttributes(pmesh);
 
-   // 8. Define a parallel finite element space on the parallel mesh. Here we
+   // 9. 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);
@@ -294,9 +309,9 @@ int main(int argc, char *argv[])
       cout << "Number of finite element unknowns: " << size << endl;
    }
 
-   // 9. Determine the list of true (i.e. parallel conforming) essential
-   //    boundary dofs. In this example, the boundary conditions are defined
-   //    based on the specific mesh and the problem type.
+   // 10. Determine the list of true (i.e. parallel conforming) essential
+   //     boundary dofs. In this example, the boundary conditions are defined
+   //     based on the specific mesh and the problem type.
    Array<int> ess_tdof_list;
    Array<int> ess_bdr;
    if (pmesh->bdr_attributes.Size())
@@ -336,11 +351,11 @@ int main(int argc, char *argv[])
    }
    fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
 
-   // 10. Setup Complex Operator convention
+   // 11. Setup Complex Operator convention
    ComplexOperator::Convention conv =
       herm_conv ? ComplexOperator::HERMITIAN : ComplexOperator::BLOCK_SYMMETRIC;
 
-   // 11. Set up the parallel linear form b(.) which corresponds to the
+   // 12. Set up the parallel linear form b(.) which corresponds to the
    //     right-hand side of the FEM linear system.
    VectorFunctionCoefficient f(dim, source);
    ParComplexLinearForm b(fespace, conv);
@@ -351,7 +366,7 @@ int main(int argc, char *argv[])
    b.Vector::operator=(0.0);
    b.Assemble();
 
-   // 12. Define the solution vector x as a parallel complex finite element grid
+   // 13. Define the solution vector x as a parallel complex finite element grid
    //     function corresponding to fespace.
    ParComplexGridFunction x(fespace);
    x = 0.0;
@@ -359,7 +374,7 @@ int main(int argc, char *argv[])
    VectorFunctionCoefficient E_Im(dim, E_bdr_data_Im);
    x.ProjectBdrCoefficientTangent(E_Re, E_Im, ess_bdr);
 
-   // 13. Set up the parallel sesquilinear form a(.,.)
+   // 14. Set up the parallel sesquilinear form a(.,.)
    //
    //     In Comp
    //     Domain:   1/mu (Curl E, Curl F) - omega^2 * epsilon (E,F)
@@ -393,19 +408,19 @@ int main(int argc, char *argv[])
    a.AddDomainIntegrator(new VectorFEMassIntegrator(restr_omeg),NULL);
 
    int cdim = (dim == 2) ? 1 : dim;
-   PMLMatrixCoefficient pml_c1_Re(cdim,detJ_inv_JT_J_Re, pml);
-   PMLMatrixCoefficient pml_c1_Im(cdim,detJ_inv_JT_J_Im, pml);
-   ScalarMatrixProductCoefficient c1_Re(muinv,pml_c1_Re);
-   ScalarMatrixProductCoefficient c1_Im(muinv,pml_c1_Im);
-   MatrixRestrictedCoefficient restr_c1_Re(c1_Re,attrPML);
-   MatrixRestrictedCoefficient restr_c1_Im(c1_Im,attrPML);
+   PMLDiagMatrixCoefficient pml_c1_Re(cdim,detJ_inv_JT_J_Re, pml);
+   PMLDiagMatrixCoefficient pml_c1_Im(cdim,detJ_inv_JT_J_Im, pml);
+   ScalarVectorProductCoefficient c1_Re(muinv,pml_c1_Re);
+   ScalarVectorProductCoefficient c1_Im(muinv,pml_c1_Im);
+   VectorRestrictedCoefficient restr_c1_Re(c1_Re,attrPML);
+   VectorRestrictedCoefficient restr_c1_Im(c1_Im,attrPML);
 
-   PMLMatrixCoefficient pml_c2_Re(dim, detJ_JT_J_inv_Re,pml);
-   PMLMatrixCoefficient pml_c2_Im(dim, detJ_JT_J_inv_Im,pml);
-   ScalarMatrixProductCoefficient c2_Re(omeg,pml_c2_Re);
-   ScalarMatrixProductCoefficient c2_Im(omeg,pml_c2_Im);
-   MatrixRestrictedCoefficient restr_c2_Re(c2_Re,attrPML);
-   MatrixRestrictedCoefficient restr_c2_Im(c2_Im,attrPML);
+   PMLDiagMatrixCoefficient pml_c2_Re(dim, detJ_JT_J_inv_Re,pml);
+   PMLDiagMatrixCoefficient pml_c2_Im(dim, detJ_JT_J_inv_Im,pml);
+   ScalarVectorProductCoefficient c2_Re(omeg,pml_c2_Re);
+   ScalarVectorProductCoefficient c2_Im(omeg,pml_c2_Im);
+   VectorRestrictedCoefficient restr_c2_Re(c2_Re,attrPML);
+   VectorRestrictedCoefficient restr_c2_Im(c2_Im,attrPML);
 
    // Integrators inside the PML region
    a.AddDomainIntegrator(new CurlCurlIntegrator(restr_c1_Re),
@@ -413,19 +428,23 @@ int main(int argc, char *argv[])
    a.AddDomainIntegrator(new VectorFEMassIntegrator(restr_c2_Re),
                          new VectorFEMassIntegrator(restr_c2_Im));
 
-   // 14. Assemble the parallel bilinear form and the corresponding linear
+   // 15. 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, etc.
+#ifndef MFEM_USE_SUPERLU
+   if (pa) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
+#endif
    a.Assemble();
 
    OperatorPtr Ah;
    Vector B, X;
    a.FormLinearSystem(ess_tdof_list, x, b, Ah, X, B);
 
-   // 15. Solve using a direct or an iterative solver
+   // 16. Solve using a direct or an iterative solver
 #ifdef MFEM_USE_SUPERLU
    {
+      if (pa) { cout << "PA not available with MFEM_USE_SUPERLU" << endl; }
       // Transform to monolithic HypreParMatrix
       HypreParMatrix *A = Ah.As<ComplexHypreParMatrix>()->GetSystemMatrix();
       SuperLURowLocMatrix SA(*A);
@@ -453,22 +472,20 @@ int main(int argc, char *argv[])
       prec.AddDomainIntegrator(new CurlCurlIntegrator(restr_muinv));
       prec.AddDomainIntegrator(new VectorFEMassIntegrator(restr_absomeg));
 
-      PMLMatrixCoefficient pml_c1_abs(cdim,detJ_inv_JT_J_abs, pml);
-      ScalarMatrixProductCoefficient c1_abs(muinv,pml_c1_abs);
-      MatrixRestrictedCoefficient restr_c1_abs(c1_abs,attrPML);
+      PMLDiagMatrixCoefficient pml_c1_abs(cdim,detJ_inv_JT_J_abs, pml);
+      ScalarVectorProductCoefficient c1_abs(muinv,pml_c1_abs);
+      VectorRestrictedCoefficient restr_c1_abs(c1_abs,attrPML);
 
-      PMLMatrixCoefficient pml_c2_abs(dim, detJ_JT_J_inv_abs,pml);
-      ScalarMatrixProductCoefficient c2_abs(absomeg,pml_c2_abs);
-      MatrixRestrictedCoefficient restr_c2_abs(c2_abs,attrPML);
+      PMLDiagMatrixCoefficient pml_c2_abs(dim, detJ_JT_J_inv_abs,pml);
+      ScalarVectorProductCoefficient c2_abs(absomeg,pml_c2_abs);
+      VectorRestrictedCoefficient restr_c2_abs(c2_abs,attrPML);
 
       prec.AddDomainIntegrator(new CurlCurlIntegrator(restr_c1_abs));
       prec.AddDomainIntegrator(new VectorFEMassIntegrator(restr_c2_abs));
 
+      if (pa) { prec.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
       prec.Assemble();
 
-      OperatorPtr PCOpAh;
-      prec.FormSystemMatrix(ess_tdof_list, PCOpAh);
-
       // 16b. Define and apply a parallel GMRES solver for AU=B with a block
       //      diagonal preconditioner based on hypre's AMS preconditioner.
       Array<int> offsets(3);
@@ -477,21 +494,41 @@ int main(int argc, char *argv[])
       offsets[2] = fespace->GetTrueVSize();
       offsets.PartialSum();
 
-      HypreAMS ams00(*PCOpAh.As<HypreParMatrix>(),fespace);
-      BlockDiagonalPreconditioner BlockAMS(offsets);
-      ScaledOperator ams11(&ams00,
-                           (conv == ComplexOperator::HERMITIAN) ? -1.0 : 1.0);
-      BlockAMS.SetDiagonalBlock(0,&ams00);
-      BlockAMS.SetDiagonalBlock(1,&ams11);
+      Operator *pc_r = nullptr;
+      Operator *pc_i = nullptr;
+      int s = (conv == ComplexOperator::HERMITIAN) ? -1.0 : 1.0;
+      if (pa)
+      {
+         // Jacobi Smoother
+         OperatorJacobiSmoother *d00 = new OperatorJacobiSmoother(prec, ess_tdof_list);
+         ScaledOperator *d11 = new ScaledOperator(d00, s);
+         pc_r = d00;
+         pc_i = d11;
+      }
+      else
+      {
+         OperatorPtr PCOpAh;
+         prec.FormSystemMatrix(ess_tdof_list, PCOpAh);
+
+         // Hypre AMS
+         HypreAMS *ams00 = new HypreAMS(*PCOpAh.As<HypreParMatrix>(), fespace);
+         ScaledOperator *ams11 = new ScaledOperator(ams00, s);
+         pc_r = ams00;
+         pc_i = ams11;
+      }
+
+      BlockDiagonalPreconditioner BlockDP(offsets);
+      BlockDP.SetDiagonalBlock(0, pc_r);
+      BlockDP.SetDiagonalBlock(1, pc_i);
 
       GMRESSolver gmres(MPI_COMM_WORLD);
       gmres.SetPrintLevel(1);
       gmres.SetKDim(200);
-      gmres.SetMaxIter(2000);
+      gmres.SetMaxIter(pa ? 5000 : 2000);
       gmres.SetRelTol(1e-5);
       gmres.SetAbsTol(0.0);
       gmres.SetOperator(*Ah);
-      gmres.SetPreconditioner(BlockAMS);
+      gmres.SetPreconditioner(BlockDP);
       gmres.Mult(B, X);
    }
 #endif
@@ -819,7 +856,7 @@ void E_bdr_data_Im(const Vector &x, Vector &E)
    }
 }
 
-void detJ_JT_J_inv_Re(const Vector &x, CartesianPML * pml, DenseMatrix &M)
+void detJ_JT_J_inv_Re(const Vector &x, CartesianPML * pml, Vector & D)
 {
    vector<complex<double>> dxs(dim);
    complex<double> det(1.0, 0.0);
@@ -830,14 +867,13 @@ void detJ_JT_J_inv_Re(const Vector &x, CartesianPML * pml, DenseMatrix &M)
       det *= dxs[i];
    }
 
-   M = 0.0;
    for (int i = 0; i < dim; ++i)
    {
-      M(i, i) = (det / pow(dxs[i], 2)).real();
+      D(i) = (det / pow(dxs[i], 2)).real();
    }
 }
 
-void detJ_JT_J_inv_Im(const Vector &x, CartesianPML * pml, DenseMatrix &M)
+void detJ_JT_J_inv_Im(const Vector &x, CartesianPML * pml, Vector & D)
 {
    vector<complex<double>> dxs(dim);
    complex<double> det = 1.0;
@@ -848,14 +884,13 @@ void detJ_JT_J_inv_Im(const Vector &x, CartesianPML * pml, DenseMatrix &M)
       det *= dxs[i];
    }
 
-   M = 0.0;
    for (int i = 0; i < dim; ++i)
    {
-      M(i, i) = (det / pow(dxs[i], 2)).imag();
+      D(i) = (det / pow(dxs[i], 2)).imag();
    }
 }
 
-void detJ_JT_J_inv_abs(const Vector &x, CartesianPML * pml, DenseMatrix &M)
+void detJ_JT_J_inv_abs(const Vector &x, CartesianPML * pml, Vector & D)
 {
    vector<complex<double>> dxs(dim);
    complex<double> det = 1.0;
@@ -866,14 +901,13 @@ void detJ_JT_J_inv_abs(const Vector &x, CartesianPML * pml, DenseMatrix &M)
       det *= dxs[i];
    }
 
-   M = 0.0;
    for (int i = 0; i < dim; ++i)
    {
-      M(i, i) = abs(det / pow(dxs[i], 2));
+      D(i) = abs(det / pow(dxs[i], 2));
    }
 }
 
-void detJ_inv_JT_J_Re(const Vector &x, CartesianPML * pml, DenseMatrix &M)
+void detJ_inv_JT_J_Re(const Vector &x, CartesianPML * pml, Vector & D)
 {
    vector<complex<double>> dxs(dim);
    complex<double> det(1.0, 0.0);
@@ -887,19 +921,18 @@ void detJ_inv_JT_J_Re(const Vector &x, CartesianPML * pml, DenseMatrix &M)
    // in the 2D case the coefficient is scalar 1/det(J)
    if (dim == 2)
    {
-      M = (1.0 / det).real();
+      D = (1.0 / det).real();
    }
    else
    {
-      M = 0.0;
       for (int i = 0; i < dim; ++i)
       {
-         M(i, i) = (pow(dxs[i], 2) / det).real();
+         D(i) = (pow(dxs[i], 2) / det).real();
       }
    }
 }
 
-void detJ_inv_JT_J_Im(const Vector &x, CartesianPML * pml, DenseMatrix &M)
+void detJ_inv_JT_J_Im(const Vector &x, CartesianPML * pml, Vector & D)
 {
    vector<complex<double>> dxs(dim);
    complex<double> det = 1.0;
@@ -912,19 +945,18 @@ void detJ_inv_JT_J_Im(const Vector &x, CartesianPML * pml, DenseMatrix &M)
 
    if (dim == 2)
    {
-      M = (1.0 / det).imag();
+      D = (1.0 / det).imag();
    }
    else
    {
-      M = 0.0;
       for (int i = 0; i < dim; ++i)
       {
-         M(i, i) = (pow(dxs[i], 2) / det).imag();
+         D(i) = (pow(dxs[i], 2) / det).imag();
       }
    }
 }
 
-void detJ_inv_JT_J_abs(const Vector &x, CartesianPML * pml, DenseMatrix &M)
+void detJ_inv_JT_J_abs(const Vector &x, CartesianPML * pml, Vector & D)
 {
    vector<complex<double>> dxs(dim);
    complex<double> det = 1.0;
@@ -937,14 +969,13 @@ void detJ_inv_JT_J_abs(const Vector &x, CartesianPML * pml, DenseMatrix &M)
 
    if (dim == 2)
    {
-      M = abs(1.0 / det);
+      D = abs(1.0 / det);
    }
    else
    {
-      M = 0.0;
       for (int i = 0; i < dim; ++i)
       {
-         M(i, i) = abs(pow(dxs[i], 2) / det);
+         D(i) = abs(pow(dxs[i], 2) / det);
       }
    }
 }

fem/bilininteg.cpp

@@ -1522,6 +1522,7 @@ void CurlCurlIntegrator::AssembleElementMatrix
    double w;
 
 #ifdef MFEM_THREAD_SAFE
+   Vector D;
    DenseMatrix curlshape(nd,dimc), curlshape_dFt(nd,dimc), M;
 #else
    curlshape.SetSize(nd,dimc);
@@ -1529,6 +1530,7 @@ void CurlCurlIntegrator::AssembleElementMatrix
 #endif
    elmat.SetSize(nd);
    if (MQ) { M.SetSize(dimc); }
+   if (DQ) { D.SetSize(dimc); }
 
    const IntegrationRule *ir = IntRule;
    if (ir == NULL)
@@ -1572,6 +1574,12 @@ void CurlCurlIntegrator::AssembleElementMatrix
          Mult(curlshape_dFt, M, curlshape);
          AddMultABt(curlshape, curlshape_dFt, elmat);
       }
+      else if (DQ)
+      {
+         DQ->Eval(D, Trans, ip);
+         D *= w;
+         AddMultADAt(curlshape_dFt, D, elmat);
+      }
       else if (Q)
       {
          w *= Q->Eval(Trans, ip);

fem/bilininteg.hpp

@@ -20,6 +20,13 @@
 namespace mfem
 {
 
+// Local maximum size of dofs and quads in 1D
+constexpr int HCURL_MAX_D1D = 5;
+constexpr int HCURL_MAX_Q1D = 6;
+
+constexpr int HDIV_MAX_D1D = 5;
+constexpr int HDIV_MAX_Q1D = 6;
+
 /// Abstract base class BilinearFormIntegrator
 class BilinearFormIntegrator : public NonlinearFormIntegrator
 {
@@ -2293,12 +2300,14 @@ class CurlCurlIntegrator: public BilinearFormIntegrator
 private:
    Vector vec, pointflux;
 #ifndef MFEM_THREAD_SAFE
+   Vector D;
    DenseMatrix curlshape, curlshape_dFt, M;
    DenseMatrix vshape, projcurl;
 #endif
 
 protected:
    Coefficient *Q;
+   VectorCoefficient *DQ;
    MatrixCoefficient *MQ;
 
    // PA extension
@@ -2307,12 +2316,17 @@ protected:
    const DofToQuad *mapsC;         ///< Not owned. DOF-to-quad map, closed.
    const GeometricFactors *geom;   ///< Not owned
    int dim, ne, nq, dofs1D, quad1D;
+   bool symmetric = true; ///< False if using a nonsymmetric matrix coefficient
 
 public:
-   CurlCurlIntegrator() { Q = NULL; MQ = NULL; }
+   CurlCurlIntegrator() { Q = NULL; DQ = NULL; MQ = NULL; }
    /// Construct a bilinear form integrator for Nedelec elements
-   CurlCurlIntegrator(Coefficient &q) : Q(&q) { MQ = NULL; }
-   CurlCurlIntegrator(MatrixCoefficient &m) : MQ(&m) { Q = NULL; }
+   CurlCurlIntegrator(Coefficient &q, const IntegrationRule *ir = NULL) :
+      BilinearFormIntegrator(ir), Q(&q) { DQ = NULL; MQ = NULL; }
+   CurlCurlIntegrator(VectorCoefficient &dq, const IntegrationRule *ir = NULL) :
+      BilinearFormIntegrator(ir), DQ(&dq) { Q = NULL; MQ = NULL; }
+   CurlCurlIntegrator(MatrixCoefficient &mq, const IntegrationRule *ir = NULL) :
+      BilinearFormIntegrator(ir), MQ(&mq) { Q = NULL; DQ = NULL; }
 
    /* Given a particular Finite Element, compute the
       element curl-curl matrix elmat */
@@ -2390,8 +2404,11 @@ protected:
    Vector pa_data;
    const DofToQuad *mapsO;         ///< Not owned. DOF-to-quad map, open.
    const DofToQuad *mapsC;         ///< Not owned. DOF-to-quad map, closed.
+   const DofToQuad *mapsOtest;     ///< Not owned. DOF-to-quad map, open.
+   const DofToQuad *mapsCtest;     ///< Not owned. DOF-to-quad map, closed.
    const GeometricFactors *geom;   ///< Not owned
-   int dim, ne, nq, dofs1D, quad1D, fetype;
+   int dim, ne, nq, dofs1D, dofs1Dtest, quad1D, trial_fetype, test_fetype;
+   bool symmetric = true; ///< False if using a nonsymmetric matrix coefficient
 
 public:
    VectorFEMassIntegrator() { Init(NULL, NULL, NULL); }
@@ -2412,6 +2429,8 @@ public:
 
    using BilinearFormIntegrator::AssemblePA;
    virtual void AssemblePA(const FiniteElementSpace &fes);
+   virtual void AssemblePA(const FiniteElementSpace &trial_fes,
+                           const FiniteElementSpace &test_fes);
    virtual void AddMultPA(const Vector &x, Vector &y) const;
    virtual void AssembleDiagonalPA(Vector& diag);
 };

fem/bilininteg_hcurl.cpp

@@ -19,10 +19,6 @@ using namespace std;
 namespace mfem
 {
 
-// Local maximum size of dofs and quads in 1D
-constexpr int HCURL_MAX_D1D = 5;
-constexpr int HCURL_MAX_Q1D = 6;
-
 // PA H(curl) Mass Assemble 2D kernel
 void PAHcurlSetup2D(const int Q1D,
                     const int coeffDim,
@@ -33,11 +29,11 @@ void PAHcurlSetup2D(const int Q1D,
                     Vector &op)
 {
    const int NQ = Q1D*Q1D;
+   const bool symmetric = (coeffDim != 4);
    auto W = w.Read();
-
    auto J = Reshape(j.Read(), NQ, 2, 2, NE);
    auto coeff = Reshape(_coeff.Read(), coeffDim, NQ, NE);
-   auto y = Reshape(op.Write(), NQ, 3, NE);
+   auto y = Reshape(op.Write(), NQ, symmetric ? 3 : 4, NE);
 
    MFEM_FORALL(e, NE,
    {
@@ -47,12 +43,39 @@ void PAHcurlSetup2D(const int Q1D,
          const double J21 = J(q,1,0,e);
          const double J12 = J(q,0,1,e);
          const double J22 = J(q,1,1,e);
-         const double c_detJ1 = W[q] * coeff(0, q, e) / ((J11*J22)-(J21*J12));
-         const double c_detJ2 = coeffDim == 2 ? W[q] * coeff(1, q, e)
-         / ((J11*J22)-(J21*J12)) : c_detJ1;
-         y(q,0,e) =  (c_detJ2*J12*J12 + c_detJ1*J22*J22); // 1,1
-         y(q,1,e) = -(c_detJ2*J12*J11 + c_detJ1*J22*J21); // 1,2
-         y(q,2,e) =  (c_detJ2*J11*J11 + c_detJ1*J21*J21); // 2,2
+
+         if (coeffDim == 3 || coeffDim == 4) // Matrix coefficient version
+         {
+            // First compute entries of R = MJ^{-T}, without det J factor.
+            const double M11 = coeff(0, q, e);
+            const double M12 = coeff(1, q, e);
+            const double M21 = symmetric ? M12 : coeff(2, q, e);
+            const double M22 = symmetric ? coeff(2, q, e) : coeff(3, q, e);
+            const double R11 = M11*J22 - M12*J12;
+            const double R21 = M21*J22 - M22*J12;
+            const double R12 = -M11*J21 + M12*J11;
+            const double R22 = -M21*J21 + M22*J11;
+
+            // Now set y to J^{-1}R.
+            const double w_detJ = W[q] / ((J11*J22)-(J21*J12));
+            y(q,0,e) = w_detJ * ( J22*R11 - J12*R21); // 1,1
+            y(q,1,e) = w_detJ * (-J21*R11 + J11*R21); // 2,1
+            y(q,2,e) = w_detJ * (symmetric ? (-J21*R12 + J11*R22) :
+            (J22*R12 - J12*R22)); // 2,2 or 1,2
+            if (!symmetric)
+            {
+               y(q,3,e) = w_detJ * (-J21*R12 + J11*R22); // 2,2
+            }
+         }
+         else  // Vector or scalar coefficient version
+         {
+            const double c_detJ1 = W[q] * coeff(0, q, e) / ((J11*J22)-(J21*J12));
+            const double c_detJ2 = (coeffDim == 2) ? W[q] * coeff(1, q, e)
+                                   / ((J11*J22)-(J21*J12)) : c_detJ1;
+            y(q,0,e) =  (c_detJ2*J12*J12 + c_detJ1*J22*J22); // 1,1
+            y(q,1,e) = -(c_detJ2*J12*J11 + c_detJ1*J22*J21); // 1,2
+            y(q,2,e) =  (c_detJ2*J11*J11 + c_detJ1*J21*J21); // 2,2
+         }
       }
    });
 }
@@ -67,10 +90,11 @@ void PAHcurlSetup3D(const int Q1D,
                     Vector &op)
 {
    const int NQ = Q1D*Q1D*Q1D;
+   const bool symmetric = (coeffDim != 9);
    auto W = w.Read();
    auto J = Reshape(j.Read(), NQ, 3, 3, NE);
    auto coeff = Reshape(_coeff.Read(), coeffDim, NQ, NE);
-   auto y = Reshape(op.Write(), NQ, 6, NE);
+   auto y = Reshape(op.Write(), NQ, symmetric ? 6 : 9, NE);
 
    MFEM_FORALL(e, NE,
    {
@@ -89,9 +113,6 @@ void PAHcurlSetup3D(const int Q1D,
          /* */               J21 * (J12 * J33 - J32 * J13) +
          /* */               J31 * (J12 * J23 - J22 * J13);
          const double w_detJ = W[q] / detJ;
-         const double D1 = coeff(0, q, e);
-         const double D2 = coeffDim == 3 ? coeff(1, q, e) : D1;
-         const double D3 = coeffDim == 3 ? coeff(2, q, e) : D1;
          // adj(J)
          const double A11 = (J22 * J33) - (J23 * J32);
          const double A12 = (J32 * J13) - (J12 * J33);
@@ -102,13 +123,66 @@ void PAHcurlSetup3D(const int Q1D,
          const double A31 = (J21 * J32) - (J31 * J22);
          const double A32 = (J31 * J12) - (J11 * J32);
          const double A33 = (J11 * J22) - (J12 * J21);
-         // detJ J^{-1} J^{-T} = (1/detJ) adj(J) D adj(J)^T
-         y(q,0,e) = w_detJ * (D1*A11*A11 + D2*A12*A12 + D3*A13*A13); // 1,1
-         y(q,1,e) = w_detJ * (D1*A11*A21 + D2*A12*A22 + D3*A13*A23); // 2,1
-         y(q,2,e) = w_detJ * (D1*A11*A31 + D2*A12*A32 + D3*A13*A33); // 3,1
-         y(q,3,e) = w_detJ * (D1*A21*A21 + D2*A22*A22 + D3*A23*A23); // 2,2
-         y(q,4,e) = w_detJ * (D1*A21*A31 + D2*A22*A32 + D3*A23*A33); // 3,2
-         y(q,5,e) = w_detJ * (D1*A31*A31 + D2*A32*A32 + D3*A33*A33); // 3,3
+
+         if (coeffDim == 6 || coeffDim == 9) // Matrix coefficient version
+         {
+            // First compute entries of R = MJ^{-T} = M adj(J)^T, without det J factor.
+            const double M11 = coeff(0, q, e);
+            const double M12 = coeff(1, q, e);
+            const double M13 = coeff(2, q, e);
+            const double M21 = (!symmetric) ? coeff(3, q, e) : M12;
+            const double M22 = (!symmetric) ? coeff(4, q, e) : coeff(3, q, e);
+            const double M23 = (!symmetric) ? coeff(5, q, e) : coeff(4, q, e);
+            const double M31 = (!symmetric) ? coeff(6, q, e) : M13;
+            const double M32 = (!symmetric) ? coeff(7, q, e) : M23;
+            const double M33 = (!symmetric) ? coeff(8, q, e) : coeff(5, q, e);
+
+            const double R11 = M11*A11 + M12*A12 + M13*A13;
+            const double R12 = M11*A21 + M12*A22 + M13*A23;
+            const double R13 = M11*A31 + M12*A32 + M13*A33;
+            const double R21 = M21*A11 + M22*A12 + M23*A13;
+            const double R22 = M21*A21 + M22*A22 + M23*A23;
+            const double R23 = M21*A31 + M22*A32 + M23*A33;
+            const double R31 = M31*A11 + M32*A12 + M33*A13;
+            const double R32 = M31*A21 + M32*A22 + M33*A23;
+            const double R33 = M31*A31 + M32*A32 + M33*A33;
+
+            // Now set y to J^{-1} R = adj(J) R
+            y(q,0,e) = w_detJ * (A11*R11 + A12*R21 + A13*R31); // 1,1
+            const double Y12 = w_detJ * (A11*R12 + A12*R22 + A13*R32);
+            y(q,1,e) = Y12; // 1,2
+            y(q,2,e) = w_detJ * (A11*R13 + A12*R23 + A13*R33); // 1,3
+
+            const double Y21 = w_detJ * (A21*R11 + A22*R21 + A23*R31);
+            const double Y22 = w_detJ * (A21*R12 + A22*R22 + A23*R32);
+            const double Y23 = w_detJ * (A21*R13 + A22*R23 + A23*R33);
+
+            const double Y33 = w_detJ * (A31*R13 + A32*R23 + A33*R33);
+
+            y(q,3,e) = symmetric ? Y22 : Y21; // 2,2 or 2,1
+            y(q,4,e) = symmetric ? Y23 : Y22; // 2,3 or 2,2
+            y(q,5,e) = symmetric ? Y33 : Y23; // 3,3 or 2,3
+
+            if (!symmetric)
+            {
+               y(q,6,e) = w_detJ * (A31*R11 + A32*R21 + A33*R31); // 3,1
+               y(q,7,e) = w_detJ * (A31*R12 + A32*R22 + A33*R32); // 3,2
+               y(q,8,e) = Y33; // 3,3
+            }
+         }
+         else  // Vector or scalar coefficient version
+         {
+            const double D1 = coeff(0, q, e);
+            const double D2 = coeffDim == 3 ? coeff(1, q, e) : D1;
+            const double D3 = coeffDim == 3 ? coeff(2, q, e) : D1;
+            // detJ J^{-1} D J^{-T} = (1/detJ) adj(J) D adj(J)^T
+            y(q,0,e) = w_detJ * (D1*A11*A11 + D2*A12*A12 + D3*A13*A13); // 1,1
+            y(q,1,e) = w_detJ * (D1*A11*A21 + D2*A12*A22 + D3*A13*A23); // 2,1
+            y(q,2,e) = w_detJ * (D1*A11*A31 + D2*A12*A32 + D3*A13*A33); // 3,1
+            y(q,3,e) = w_detJ * (D1*A21*A21 + D2*A22*A22 + D3*A23*A23); // 2,2
+            y(q,4,e) = w_detJ * (D1*A21*A31 + D2*A22*A32 + D3*A23*A33); // 3,2
+            y(q,5,e) = w_detJ * (D1*A31*A31 + D2*A32*A32 + D3*A33*A33); // 3,3
+         }
       }
    });
 }
@@ -116,6 +190,7 @@ void PAHcurlSetup3D(const int Q1D,
 void PAHcurlMassApply2D(const int D1D,
                         const int Q1D,
                         const int NE,
+                        const bool symmetric,
                         const Array<double> &_Bo,
                         const Array<double> &_Bc,
                         const Array<double> &_Bot,
@@ -132,7 +207,7 @@ void PAHcurlMassApply2D(const int D1D,
    auto Bc = Reshape(_Bc.Read(), Q1D, D1D);
    auto Bot = Reshape(_Bot.Read(), D1D-1, Q1D);
    auto Bct = Reshape(_Bct.Read(), D1D, Q1D);
-   auto op = Reshape(_op.Read(), Q1D, Q1D, 3, NE);
+   auto op = Reshape(_op.Read(), Q1D, Q1D, symmetric ? 3 : 4, NE);
    auto x = Reshape(_x.Read(), 2*(D1D-1)*D1D, NE);
    auto y = Reshape(_y.ReadWrite(), 2*(D1D-1)*D1D, NE);
 
@@ -194,12 +269,13 @@ void PAHcurlMassApply2D(const int D1D,
          for (int qx = 0; qx < Q1D; ++qx)
          {
             const double O11 = op(qx,qy,0,e);
-            const double O12 = op(qx,qy,1,e);
-            const double O22 = op(qx,qy,2,e);
+            const double O21 = op(qx,qy,1,e);
+            const double O12 = symmetric ? O21 : op(qx,qy,2,e);
+            const double O22 = symmetric ? op(qx,qy,2,e) : op(qx,qy,3,e);
             const double massX = mass[qy][qx][0];
             const double massY = mass[qy][qx][1];
             mass[qy][qx][0] = (O11*massX)+(O12*massY);
-            mass[qy][qx][1] = (O12*massX)+(O22*massY);
+            mass[qy][qx][1] = (O21*massX)+(O22*massY);
          }
       }
 
@@ -215,7 +291,7 @@ void PAHcurlMassApply2D(const int D1D,
             double massX[MAX_D1D];
             for (int dx = 0; dx < D1Dx; ++dx)
             {
-               massX[dx] = 0;
+               massX[dx] = 0.0;
             }
             for (int qx = 0; qx < Q1D; ++qx)
             {
@@ -244,6 +320,7 @@ void PAHcurlMassApply2D(const int D1D,
 void PAHcurlMassAssembleDiagonal2D(const int D1D,
                                    const int Q1D,
                                    const int NE,
+                                   const bool symmetric,
                                    const Array<double> &_Bo,
                                    const Array<double> &_Bc,
                                    const Vector &_op,
@@ -254,7 +331,7 @@ void PAHcurlMassAssembleDiagonal2D(const int D1D,
 
    auto Bo = Reshape(_Bo.Read(), Q1D, D1D-1);
    auto Bc = Reshape(_Bc.Read(), Q1D, D1D);
-   auto op = Reshape(_op.Read(), Q1D, Q1D, 3, NE);
+   auto op = Reshape(_op.Read(), Q1D, Q1D, symmetric ? 3 : 4, NE);
    auto diag = Reshape(_diag.ReadWrite(), 2*(D1D-1)*D1D, NE);
 
    MFEM_FORALL(e, NE,
@@ -277,7 +354,8 @@ void PAHcurlMassAssembleDiagonal2D(const int D1D,
                {
                   const double wy = (c == 1) ? Bo(qy,dy) : Bc(qy,dy);
 
-                  mass[qx] += wy * wy * ((c == 0) ? op(qx,qy,0,e) : op(qx,qy,2,e));
+                  mass[qx] += wy * wy * ((c == 0) ? op(qx,qy,0,e) :
+                  op(qx,qy,symmetric ? 2 : 3, e));
                }
             }
 
@@ -299,6 +377,7 @@ void PAHcurlMassAssembleDiagonal2D(const int D1D,
 void PAHcurlMassAssembleDiagonal3D(const int D1D,
                                    const int Q1D,
                                    const int NE,
+                                   const bool symmetric,
                                    const Array<double> &_Bo,
                                    const Array<double> &_Bc,
                                    const Vector &_op,
@@ -313,7 +392,7 @@ void PAHcurlMassAssembleDiagonal3D(const int D1D,
 
    auto Bo = Reshape(_Bo.Read(), Q1D, D1D-1);
    auto Bc = Reshape(_Bc.Read(), Q1D, D1D);
-   auto op = Reshape(_op.Read(), Q1D, Q1D, Q1D, 6, NE);
+   auto op = Reshape(_op.Read(), Q1D, Q1D, Q1D, symmetric ? 6 : 9, NE);
    auto diag = Reshape(_diag.ReadWrite(), 3*(D1D-1)*D1D*D1D, NE);
 
    MFEM_FORALL(e, NE,
@@ -326,7 +405,8 @@ void PAHcurlMassAssembleDiagonal3D(const int D1D,
          const int D1Dy = (c == 1) ? D1D - 1 : D1D;
          const int D1Dx = (c == 0) ? D1D - 1 : D1D;
 
-         const int opc = (c == 0) ? 0 : ((c == 1) ? 3 : 5);
+         const int opc = (c == 0) ? 0 : ((c == 1) ? (symmetric ? 3 : 4) :
+         (symmetric ? 5 : 8));
 
          double mass[MAX_Q1D];
 
@@ -369,6 +449,7 @@ void PAHcurlMassAssembleDiagonal3D(const int D1D,
 void PAHcurlMassApply3D(const int D1D,
                         const int Q1D,
                         const int NE,
+                        const bool symmetric,
                         const Array<double> &_Bo,
                         const Array<double> &_Bc,
                         const Array<double> &_Bot,
@@ -388,7 +469,7 @@ void PAHcurlMassApply3D(const int D1D,
    auto Bc = Reshape(_Bc.Read(), Q1D, D1D);
    auto Bot = Reshape(_Bot.Read(), D1D-1, Q1D);
    auto Bct = Reshape(_Bct.Read(), D1D, Q1D);
-   auto op = Reshape(_op.Read(), Q1D, Q1D, Q1D, 6, NE);
+   auto op = Reshape(_op.Read(), Q1D, Q1D, Q1D, symmetric ? 6 : 9, NE);
    auto x = Reshape(_x.Read(), 3*(D1D-1)*D1D*D1D, NE);
    auto y = Reshape(_y.ReadWrite(), 3*(D1D-1)*D1D*D1D, NE);
 
@@ -483,15 +564,18 @@ void PAHcurlMassApply3D(const int D1D,
                const double O11 = op(qx,qy,qz,0,e);
                const double O12 = op(qx,qy,qz,1,e);
                const double O13 = op(qx,qy,qz,2,e);
-               const double O22 = op(qx,qy,qz,3,e);
-               const double O23 = op(qx,qy,qz,4,e);
-               const double O33 = op(qx,qy,qz,5,e);
+               const double O21 = symmetric ? O12 : op(qx,qy,qz,3,e);
+               const double O22 = symmetric ? op(qx,qy,qz,3,e) : op(qx,qy,qz,4,e);
+               const double O23 = symmetric ? op(qx,qy,qz,4,e) : op(qx,qy,qz,5,e);
+               const double O31 = symmetric ? O13 : op(qx,qy,qz,6,e);
+               const double O32 = symmetric ? O23 : op(qx,qy,qz,7,e);
+               const double O33 = symmetric ? op(qx,qy,qz,5,e) : op(qx,qy,qz,8,e);
                const double massX = mass[qz][qy][qx][0];
                const double massY = mass[qz][qy][qx][1];
                const double massZ = mass[qz][qy][qx][2];
                mass[qz][qy][qx][0] = (O11*massX)+(O12*massY)+(O13*massZ);
-               mass[qz][qy][qx][1] = (O12*massX)+(O22*massY)+(O23*massZ);
-               mass[qz][qy][qx][2] = (O13*massX)+(O23*massY)+(O33*massZ);
+               mass[qz][qy][qx][1] = (O21*massX)+(O22*massY)+(O23*massZ);
+               mass[qz][qy][qx][2] = (O31*massX)+(O32*massY)+(O33*massZ);
             }
          }
       }
@@ -512,7 +596,7 @@ void PAHcurlMassApply3D(const int D1D,
             {
                for (int dx = 0; dx < D1Dx; ++dx)
                {
-                  massXY[dy][dx] = 0;
+                  massXY[dy][dx] = 0.0;
                }
             }
             for (int qy = 0; qy < Q1D; ++qy)
@@ -594,10 +678,12 @@ static void PACurlCurlSetup3D(const int Q1D,
                               Vector &op)
 {
    const int NQ = Q1D*Q1D*Q1D;
+   const bool symmetric = (coeffDim != 9);
    auto W = w.Read();
    auto J = Reshape(j.Read(), NQ, 3, 3, NE);
    auto coeff = Reshape(_coeff.Read(), coeffDim, NQ, NE);
-   auto y = Reshape(op.Write(), NQ, 6, NE);
+   auto y = Reshape(op.Write(), NQ, symmetric ? 6 : 9, NE);
+
    MFEM_FORALL(e, NE,
    {
       for (int q = 0; q < NQ; ++q)
@@ -615,19 +701,69 @@ static void PACurlCurlSetup3D(const int Q1D,
          /* */               J21 * (J12 * J33 - J32 * J13) +
          /* */               J31 * (J12 * J23 - J22 * J13);
 
-         const double D1 = coeff(0, q, e);
-         const double D2 = coeffDim == 3 ? coeff(1, q, e) : D1;
-         const double D3 = coeffDim == 3 ? coeff(2, q, e) : D1;
-
-         // set y to the 6 entries of J^T D J / det^2
          const double c_detJ = W[q] / detJ;
 
-         y(q,0,e) = c_detJ * (D1*J11*J11 + D2*J21*J21 + D3*J31*J31); // 1,1
-         y(q,1,e) = c_detJ * (D1*J11*J12 + D2*J21*J22 + D3*J31*J32); // 1,2
-         y(q,2,e) = c_detJ * (D1*J11*J13 + D2*J21*J23 + D3*J31*J33); // 1,3
-         y(q,3,e) = c_detJ * (D1*J12*J12 + D2*J22*J22 + D3*J32*J32); // 2,2
-         y(q,4,e) = c_detJ * (D1*J12*J13 + D2*J22*J23 + D3*J32*J33); // 2,3
-         y(q,5,e) = c_detJ * (D1*J13*J13 + D2*J23*J23 + D3*J33*J33); // 3,3
+         if (coeffDim == 6 || coeffDim == 9) // Matrix coefficient version
+         {
+            // Set y to the 6 or 9 entries of J^T M J / det
+            const double M11 = coeff(0, q, e);
+            const double M12 = coeff(1, q, e);
+            const double M13 = coeff(2, q, e);
+            const double M21 = (!symmetric) ? coeff(3, q, e) : M12;
+            const double M22 = (!symmetric) ? coeff(4, q, e) : coeff(3, q, e);
+            const double M23 = (!symmetric) ? coeff(5, q, e) : coeff(4, q, e);
+            const double M31 = (!symmetric) ? coeff(6, q, e) : M13;
+            const double M32 = (!symmetric) ? coeff(7, q, e) : M23;
+            const double M33 = (!symmetric) ? coeff(8, q, e) : coeff(5, q, e);
+
+            // First compute R = MJ
+            const double R11 = M11*J11 + M12*J21 + M13*J31;
+            const double R12 = M11*J12 + M12*J22 + M13*J32;
+            const double R13 = M11*J13 + M12*J23 + M13*J33;
+            const double R21 = M21*J11 + M22*J21 + M23*J31;
+            const double R22 = M21*J12 + M22*J22 + M23*J32;
+            const double R23 = M21*J13 + M22*J23 + M23*J33;
+            const double R31 = M31*J11 + M32*J21 + M33*J31;
+            const double R32 = M31*J12 + M32*J22 + M33*J32;
+            const double R33 = M31*J13 + M32*J23 + M33*J33;
+
+            // Now set y to J^T R / det
+            y(q,0,e) = c_detJ * (J11*R11 + J21*R21 + J31*R31); // 1,1
+            const double Y12 = c_detJ * (J11*R12 + J21*R22 + J31*R32);
+            y(q,1,e) = Y12; // 1,2
+            y(q,2,e) = c_detJ * (J11*R13 + J21*R23 + J31*R33); // 1,3
+
+            const double Y21 = c_detJ * (J12*R11 + J22*R21 + J32*R31);
+            const double Y22 = c_detJ * (J12*R12 + J22*R22 + J32*R32);
+            const double Y23 = c_detJ * (J12*R13 + J22*R23 + J32*R33);
+
+            const double Y33 = c_detJ * (J13*R13 + J23*R23 + J33*R33);
+
+            y(q,3,e) = symmetric ? Y22 : Y21; // 2,2 or 2,1
+            y(q,4,e) = symmetric ? Y23 : Y22; // 2,3 or 2,2
+            y(q,5,e) = symmetric ? Y33 : Y23; // 3,3 or 2,3
+
+            if (!symmetric)
+            {
+               y(q,6,e) = c_detJ * (J13*R11 + J23*R21 + J33*R31); // 3,1
+               y(q,7,e) = c_detJ * (J13*R12 + J23*R22 + J33*R32); // 3,2
+               y(q,8,e) = Y33; // 3,3
+            }
+         }
+         else  // Vector or scalar coefficient version
+         {
+            // Set y to the 6 entries of J^T D J / det^2
+            const double D1 = coeff(0, q, e);
+            const double D2 = coeffDim == 3 ? coeff(1, q, e) : D1;
+            const double D3 = coeffDim == 3 ? coeff(2, q, e) : D1;
+
+            y(q,0,e) = c_detJ * (D1*J11*J11 + D2*J21*J21 + D3*J31*J31); // 1,1
+            y(q,1,e) = c_detJ * (D1*J11*J12 + D2*J21*J22 + D3*J31*J32); // 1,2
+            y(q,2,e) = c_detJ * (D1*J11*J13 + D2*J21*J23 + D3*J31*J33); // 1,3
+            y(q,3,e) = c_detJ * (D1*J12*J12 + D2*J22*J22 + D3*J32*J32); // 2,2
+            y(q,4,e) = c_detJ * (D1*J12*J13 + D2*J22*J23 + D3*J32*J33); // 2,3
+            y(q,5,e) = c_detJ * (D1*J13*J13 + D2*J23*J23 + D3*J33*J33); // 3,3
+         }
       }
    });
 }
@@ -652,6 +788,8 @@ void CurlCurlIntegrator::AssemblePA(const FiniteElementSpace &fes)
    dim = mesh->Dimension();
    MFEM_VERIFY(dim == 2 || dim == 3, "");
 
+   const int dimc = (dim == 3) ? 3 : 1;
+
    ne = fes.GetNE();
    geom = mesh->GetGeometricFactors(*ir, GeometricFactors::JACOBIANS);
    mapsC = &el->GetDofToQuad(*ir, DofToQuad::TENSOR);
@@ -661,36 +799,103 @@ void CurlCurlIntegrator::AssemblePA(const FiniteElementSpace &fes)
 
    MFEM_VERIFY(dofs1D == mapsO->ndof + 1 && quad1D == mapsO->nqpt, "");
 
-   const int ndata = (dim == 2) ? 1 : 6;
+   const int MQsymmDim = MQ ? (MQ->GetWidth() * (MQ->GetWidth() + 1)) / 2 : 0;
+   const int MQfullDim = MQ ? (MQ->GetHeight() * MQ->GetWidth()) : 0;
+   const int MQdim = MQ ? (MQ->IsSymmetric() ? MQsymmDim : MQfullDim) : 0;
+   const int coeffDim = MQ ? MQdim : (DQ ? DQ->GetVDim() : 1);
+
+   symmetric = MQ ? MQ->IsSymmetric() : true;
+
+   const int symmDims = (dims * (dims + 1)) / 2; // 1x1: 1, 2x2: 3, 3x3: 6
+   const int ndata = (dim == 2) ? 1 : (symmetric ? symmDims : MQfullDim);
    pa_data.SetSize(ndata * nq * ne, Device::GetMemoryType());
 
-   Vector coeff(ne * nq);
+   Vector coeff(coeffDim * ne * nq);
    coeff = 1.0;
-   if (Q)
+   auto coeffh = Reshape(coeff.HostWrite(), coeffDim, nq, ne);
+   if (Q || DQ || MQ)
    {
+      Vector D(DQ ? coeffDim : 0);
+      DenseMatrix M;
+      Vector Msymm;
+      if (MQ)
+      {
+         if (symmetric)
+         {
+            Msymm.SetSize(MQsymmDim);
+         }
+         else
+         {
+            M.SetSize(dimc);
+         }
+      }
+
+      if (DQ)
+      {
+         MFEM_VERIFY(coeffDim == dimc, "");
+      }
+      if (MQ)
+      {
+         MFEM_VERIFY(coeffDim == MQdim, "");
+         MFEM_VERIFY(MQ->GetHeight() == dimc && MQ->GetWidth() == dimc, "");
+      }
+
       for (int e=0; e<ne; ++e)
       {
          ElementTransformation *tr = mesh->GetElementTransformation(e);
          for (int p=0; p<nq; ++p)
          {
-            coeff[p + (e * nq)] = Q->Eval(*tr, ir->IntPoint(p));
+            if (MQ)
+            {
+               if (MQ->IsSymmetric())
+               {
+                  MQ->EvalSymmetric(Msymm, *tr, ir->IntPoint(p));
+
+                  for (int i=0; i<MQsymmDim; ++i)
+                  {
+                     coeffh(i, p, e) = Msymm[i];
+                  }
+               }
+               else
+               {
+                  MQ->Eval(M, *tr, ir->IntPoint(p));
+
+                  for (int i=0; i<dimc; ++i)
+                     for (int j=0; j<dimc; ++j)
+                     {
+                        coeffh(j+(i*dimc), p, e) = M(i,j);
+                     }
+               }
+            }
+            else if (DQ)
+            {
+               DQ->Eval(D, *tr, ir->IntPoint(p));
+               for (int i=0; i<coeffDim; ++i)
+               {
+                  coeffh(i, p, e) = D[i];
+               }
+            }
+            else
+            {
+               coeffh(0, p, e) = Q->Eval(*tr, ir->IntPoint(p));
+            }
          }
       }
    }
 
-   if (el->GetDerivType() == mfem::FiniteElement::CURL && dim == 3)
+   if (el->GetDerivType() != mfem::FiniteElement::CURL)
    {
-      PACurlCurlSetup3D(quad1D, 1, ne, ir->GetWeights(), geom->J,
-                        coeff, pa_data);
+      MFEM_ABORT("Unknown kernel.");
    }
-   else if (el->GetDerivType() == mfem::FiniteElement::CURL && dim == 2)
+
+   if (dim == 3)
    {
-      PACurlCurlSetup2D(quad1D, ne, ir->GetWeights(), geom->J,
-                        coeff, pa_data);
+      PACurlCurlSetup3D(quad1D, coeffDim, ne, ir->GetWeights(), geom->J, coeff,
+                        pa_data);
    }
    else
    {
-      MFEM_ABORT("Unknown kernel.");
+      PACurlCurlSetup2D(quad1D, ne, ir->GetWeights(), geom->J, coeff, pa_data);
    }
 }
 
@@ -727,7 +932,7 @@ static void PACurlCurlApply2D(const int D1D,
       {
          for (int qx = 0; qx < Q1D; ++qx)
          {
-            curl[qy][qx] = 0;
+            curl[qy][qx] = 0.0;
          }
       }
 
@@ -789,7 +994,7 @@ static void PACurlCurlApply2D(const int D1D,
             double gradX[MAX_D1D];
             for (int dx = 0; dx < D1Dx; ++dx)
             {
-               gradX[dx] = 0;
+               gradX[dx] = 0.0;
             }
             for (int qx = 0; qx < Q1D; ++qx)
             {
@@ -817,6 +1022,7 @@ static void PACurlCurlApply2D(const int D1D,
 template<int MAX_D1D = HCURL_MAX_D1D, int MAX_Q1D = HCURL_MAX_Q1D>
 static void PACurlCurlApply3D(const int D1D,
                               const int Q1D,
+                              const bool symmetric,
                               const int NE,
                               const Array<double> &_Bo,
                               const Array<double> &_Bc,
@@ -844,7 +1050,7 @@ static void PACurlCurlApply3D(const int D1D,
    auto Bct = Reshape(_Bct.Read(), D1D, Q1D);
    auto Gc = Reshape(_Gc.Read(), Q1D, D1D);
    auto Gct = Reshape(_Gct.Read(), D1D, Q1D);
-   auto op = Reshape(_op.Read(), Q1D, Q1D, Q1D, 6, NE);
+   auto op = Reshape(_op.Read(), Q1D, Q1D, Q1D, (symmetric ? 6 : 9), NE);
    auto x = Reshape(_x.Read(), 3*(D1D-1)*D1D*D1D, NE);
    auto y = Reshape(_y.ReadWrite(), 3*(D1D-1)*D1D*D1D, NE);
 
@@ -1087,15 +1293,18 @@ static void PACurlCurlApply3D(const int D1D,
                const double O11 = op(qx,qy,qz,0,e);
                const double O12 = op(qx,qy,qz,1,e);
                const double O13 = op(qx,qy,qz,2,e);
-               const double O22 = op(qx,qy,qz,3,e);
-               const double O23 = op(qx,qy,qz,4,e);
-               const double O33 = op(qx,qy,qz,5,e);
+               const double O21 = symmetric ? O12 : op(qx,qy,qz,3,e);
+               const double O22 = symmetric ? op(qx,qy,qz,3,e) : op(qx,qy,qz,4,e);
+               const double O23 = symmetric ? op(qx,qy,qz,4,e) : op(qx,qy,qz,5,e);
+               const double O31 = symmetric ? O13 : op(qx,qy,qz,6,e);
+               const double O32 = symmetric ? O23 : op(qx,qy,qz,7,e);
+               const double O33 = symmetric ? op(qx,qy,qz,5,e) : op(qx,qy,qz,8,e);
 
                const double c1 = (O11 * curl[qz][qy][qx][0]) + (O12 * curl[qz][qy][qx][1]) +
                                  (O13 * curl[qz][qy][qx][2]);
-               const double c2 = (O12 * curl[qz][qy][qx][0]) + (O22 * curl[qz][qy][qx][1]) +
+               const double c2 = (O21 * curl[qz][qy][qx][0]) + (O22 * curl[qz][qy][qx][1]) +
                                  (O23 * curl[qz][qy][qx][2]);
-               const double c3 = (O13 * curl[qz][qy][qx][0]) + (O23 * curl[qz][qy][qx][1]) +
+               const double c3 = (O31 * curl[qz][qy][qx][0]) + (O32 * curl[qz][qy][qx][1]) +
                                  (O33 * curl[qz][qy][qx][2]);
 
                curl[qz][qy][qx][0] = c1;
@@ -1326,7 +1535,7 @@ void CurlCurlIntegrator::AddMultPA(const Vector &x, Vector &y) const
 {
    if (dim == 3)
    {
-      PACurlCurlApply3D(dofs1D, quad1D, ne, mapsO->B, mapsC->B, mapsO->Bt,
+      PACurlCurlApply3D(dofs1D, quad1D, symmetric, ne, mapsO->B, mapsC->B, mapsO->Bt,
                         mapsC->Bt, mapsC->G, mapsC->Gt, pa_data, x, y);
    }
    else if (dim == 2)
@@ -1757,7 +1966,7 @@ void PAHcurlH1Apply3D(const int D1D,
             {
                for (int dx = 0; dx < D1Dx; ++dx)
                {
-                  massXY[dy][dx] = 0;
+                  massXY[dy][dx] = 0.0;
                }
             }
             for (int qy = 0; qy < Q1D; ++qy)
@@ -1900,7 +2109,7 @@ void PAHcurlH1Apply2D(const int D1D,
             double massX[MAX_D1D];
             for (int dx = 0; dx < D1Dx; ++dx)
             {
-               massX[dx] = 0;
+               massX[dx] = 0.0;
             }
             for (int qx = 0; qx < Q1D; ++qx)
             {
@@ -1926,15 +2135,14 @@ void PAHcurlH1Apply2D(const int D1D,
    }); // end of element loop
 }
 
-// PA H(curl) Mass Assemble 3D kernel
-static void PAHcurlL2Setup3D(const int Q1D,
-                             const int coeffDim,
-                             const int NE,
-                             const Array<double> &w,
-                             Vector &_coeff,
-                             Vector &op)
+// PA H(curl) assemble kernel
+void PAHcurlL2Setup(const int NQ,
+                    const int coeffDim,
+                    const int NE,
+                    const Array<double> &w,
+                    Vector &_coeff,
+                    Vector &op)
 {
-   const int NQ = Q1D*Q1D*Q1D;
    auto W = w.Read();
    auto coeff = Reshape(_coeff.Read(), coeffDim, NQ, NE);
    auto y = Reshape(op.Write(), coeffDim, NQ, NE);
@@ -2035,7 +2243,7 @@ void MixedVectorCurlIntegrator::AssemblePA(const FiniteElementSpace &trial_fes,
    if (testType == mfem::FiniteElement::CURL &&
        trialType == mfem::FiniteElement::CURL && dim == 3)
    {
-      PAHcurlL2Setup3D(quad1D, coeffDim, ne, ir->GetWeights(), coeff, pa_data);
+      PAHcurlL2Setup(nq, coeffDim, ne, ir->GetWeights(), coeff, pa_data);
    }
    else if (testType == mfem::FiniteElement::DIV &&
             trialType == mfem::FiniteElement::CURL && dim == 3 &&
@@ -2346,7 +2554,7 @@ static void PAHcurlL2Apply3D(const int D1D,
             {
                for (int dx = 0; dx < D1Dx; ++dx)
                {
-                  massXY[dy][dx] = 0;
+                  massXY[dy][dx] = 0.0;
                }
             }
             for (int qy = 0; qy < Q1D; ++qy)
@@ -2354,7 +2562,7 @@ static void PAHcurlL2Apply3D(const int D1D,
                double massX[MAX_D1D];
                for (int dx = 0; dx < D1Dx; ++dx)
                {
-                  massX[dx] = 0;
+                  massX[dx] = 0.0;
                }
                for (int qx = 0; qx < Q1D; ++qx)
                {
@@ -2425,7 +2633,7 @@ static void PAHcurlHdivApply3D(const int D1D,
    auto Gc = Reshape(_Gc.Read(), Q1D, D1D);
    auto op = Reshape(_op.Read(), Q1D, Q1D, Q1D, 6, NE);
    auto x = Reshape(_x.Read(), 3*(D1D-1)*D1D*D1D, NE);
-   auto y = Reshape(_y.ReadWrite(), 3*(D1Dtest-1)*(D1Dtest-1)*D1D, NE);
+   auto y = Reshape(_y.ReadWrite(), 3*(D1Dtest-1)*(D1Dtest-1)*D1Dtest, NE);
 
    MFEM_FORALL(e, NE,
    {
@@ -2700,7 +2908,7 @@ static void PAHcurlHdivApply3D(const int D1D,
             {
                for (int dx = 0; dx < D1Dx; ++dx)
                {
-                  massXY[dy][dx] = 0;
+                  massXY[dy][dx] = 0.0;
                }
             }
             for (int qy = 0; qy < Q1D; ++qy)
@@ -2708,7 +2916,7 @@ static void PAHcurlHdivApply3D(const int D1D,
                double massX[HCURL_MAX_D1D];
                for (int dx = 0; dx < D1Dx; ++dx)
                {
-                  massX[dx] = 0;
+                  massX[dx] = 0.0;
                }
                for (int qx = 0; qx < Q1D; ++qx)
                {
@@ -2844,7 +3052,7 @@ void MixedVectorWeakCurlIntegrator::AssemblePA(const FiniteElementSpace
 
    if (trialType == mfem::FiniteElement::CURL && dim == 3)
    {
-      PAHcurlL2Setup3D(quad1D, coeffDim, ne, ir->GetWeights(), coeff, pa_data);
+      PAHcurlL2Setup(nq, coeffDim, ne, ir->GetWeights(), coeff, pa_data);
    }
    else
    {

fem/bilininteg_hdiv.cpp

@@ -23,11 +23,6 @@ using namespace std;
 namespace mfem
 {
 
-// Local maximum size of dofs and quads in 1D
-constexpr int HDIV_MAX_D1D = 5;
-constexpr int HDIV_MAX_Q1D = 6;
-
-
 // PA H(div) Mass Assemble 2D kernel
 void PAHdivSetup2D(const int Q1D,
                    const int NE,

fem/bilininteg_vectorfe.cpp

@@ -34,6 +34,7 @@ void PAHcurlSetup3D(const int Q1D,
 void PAHcurlMassAssembleDiagonal2D(const int D1D,
                                    const int Q1D,
                                    const int NE,
+                                   const bool symmetric,
                                    const Array<double> &_Bo,
                                    const Array<double> &_Bc,
                                    const Vector &_op,
@@ -42,6 +43,7 @@ void PAHcurlMassAssembleDiagonal2D(const int D1D,
 void PAHcurlMassAssembleDiagonal3D(const int D1D,
                                    const int Q1D,
                                    const int NE,
+                                   const bool symmetric,
                                    const Array<double> &_Bo,
                                    const Array<double> &_Bc,
                                    const Vector &_op,
@@ -50,6 +52,7 @@ void PAHcurlMassAssembleDiagonal3D(const int D1D,
 void PAHcurlMassApply2D(const int D1D,
                         const int Q1D,
                         const int NE,
+                        const bool symmetric,
                         const Array<double> &_Bo,
                         const Array<double> &_Bc,
                         const Array<double> &_Bot,
@@ -61,6 +64,7 @@ void PAHcurlMassApply2D(const int D1D,
 void PAHcurlMassApply3D(const int D1D,
                         const int Q1D,
                         const int NE,
+                        const bool symmetric,
                         const Array<double> &_Bo,
                         const Array<double> &_Bc,
                         const Array<double> &_Bot,
@@ -143,20 +147,573 @@ void PAHdivMassApply3D(const int D1D,
                        const Vector &_x,
                        Vector &_y);
 
+void PAHcurlL2Setup(const int NQ,
+                    const int coeffDim,
+                    const int NE,
+                    const Array<double> &w,
+                    Vector &_coeff,
+                    Vector &op);
+
+// PA H(curl) x H(div) mass assemble 3D kernel, with factor
+// dF^{-1} C dF for a vector or matrix coefficient C.
+// If transpose, use dF^T C dF^{-T} for H(div) x H(curl).
+void PAHcurlHdivSetup3D(const int Q1D,
+                        const int coeffDim,
+                        const int NE,
+                        const bool transpose,
+                        const Array<double> &_w,
+                        const Vector &j,
+                        Vector &_coeff,
+                        Vector &op)
+{
+   const int NQ = Q1D*Q1D*Q1D;
+   const bool symmetric = (coeffDim != 9);
+   auto W = _w.Read();
+   auto J = Reshape(j.Read(), NQ, 3, 3, NE);
+   auto coeff = Reshape(_coeff.Read(), coeffDim, NQ, NE);
+   auto y = Reshape(op.Write(), 9, NQ, NE);
+
+   const int i11 = 0;
+   const int i12 = transpose ? 3 : 1;
+   const int i13 = transpose ? 6 : 2;
+   const int i21 = transpose ? 1 : 3;
+   const int i22 = 4;
+   const int i23 = transpose ? 7 : 5;
+   const int i31 = transpose ? 2 : 6;
+   const int i32 = transpose ? 5 : 7;
+   const int i33 = 8;
+
+   MFEM_FORALL(e, NE,
+   {
+      for (int q = 0; q < NQ; ++q)
+      {
+         const double J11 = J(q,0,0,e);
+         const double J21 = J(q,1,0,e);
+         const double J31 = J(q,2,0,e);
+         const double J12 = J(q,0,1,e);
+         const double J22 = J(q,1,1,e);
+         const double J32 = J(q,2,1,e);
+         const double J13 = J(q,0,2,e);
+         const double J23 = J(q,1,2,e);
+         const double J33 = J(q,2,2,e);
+         const double detJ = J11 * (J22 * J33 - J32 * J23) -
+         /* */               J21 * (J12 * J33 - J32 * J13) +
+         /* */               J31 * (J12 * J23 - J22 * J13);
+         const double w_detJ = W[q] / detJ;
+         // adj(J)
+         const double A11 = (J22 * J33) - (J23 * J32);
+         const double A12 = (J32 * J13) - (J12 * J33);
+         const double A13 = (J12 * J23) - (J22 * J13);
+         const double A21 = (J31 * J23) - (J21 * J33);
+         const double A22 = (J11 * J33) - (J13 * J31);
+         const double A23 = (J21 * J13) - (J11 * J23);
+         const double A31 = (J21 * J32) - (J31 * J22);
+         const double A32 = (J31 * J12) - (J11 * J32);
+         const double A33 = (J11 * J22) - (J12 * J21);
+
+         if (coeffDim == 6 || coeffDim == 9) // Matrix coefficient version
+         {
+            // First compute entries of R = MJ
+            const double M11 = (!symmetric) ? coeff(i11, q, e) : coeff(0, q, e);
+            const double M12 = (!symmetric) ? coeff(i12, q, e) : coeff(1, q, e);
+            const double M13 = (!symmetric) ? coeff(i13, q, e) : coeff(2, q, e);
+            const double M21 = (!symmetric) ? coeff(i21, q, e) : M12;
+            const double M22 = (!symmetric) ? coeff(i22, q, e) : coeff(3, q, e);
+            const double M23 = (!symmetric) ? coeff(i23, q, e) : coeff(4, q, e);
+            const double M31 = (!symmetric) ? coeff(i31, q, e) : M13;
+            const double M32 = (!symmetric) ? coeff(i32, q, e) : M23;
+            const double M33 = (!symmetric) ? coeff(i33, q, e) : coeff(5, q, e);
+
+            const double R11 = M11*J11 + M12*J12 + M13*J13;
+            const double R12 = M11*J21 + M12*J22 + M13*J23;
+            const double R13 = M11*J31 + M12*J32 + M13*J33;
+            const double R21 = M21*J11 + M22*J12 + M23*J13;
+            const double R22 = M21*J21 + M22*J22 + M23*J23;
+            const double R23 = M21*J31 + M22*J32 + M23*J33;
+            const double R31 = M31*J11 + M32*J12 + M33*J13;
+            const double R32 = M31*J21 + M32*J22 + M33*J23;
+            const double R33 = M31*J31 + M32*J32 + M33*J33;
+
+            // Now set y to detJ J^{-1} R = adj(J) R
+            y(i11,q,e) = w_detJ * (A11*R11 + A12*R21 + A13*R31); // 1,1
+            y(i12,q,e) = w_detJ * (A11*R12 + A12*R22 + A13*R32); // 1,2
+            y(i13,q,e) = w_detJ * (A11*R13 + A12*R23 + A13*R33); // 1,3
+            y(i21,q,e) = w_detJ * (A21*R11 + A22*R21 + A23*R31); // 2,1
+            y(i22,q,e) = w_detJ * (A21*R12 + A22*R22 + A23*R32); // 2,2
+            y(i23,q,e) = w_detJ * (A21*R13 + A22*R23 + A23*R33); // 2,3
+            y(i31,q,e) = w_detJ * (A31*R11 + A32*R21 + A33*R31); // 3,1
+            y(i32,q,e) = w_detJ * (A31*R12 + A32*R22 + A33*R32); // 3,2
+            y(i33,q,e) = w_detJ * (A31*R13 + A32*R23 + A33*R33); // 3,3
+         }
+         else if (coeffDim == 3)  // Vector coefficient version
+         {
+            const double D1 = coeff(0, q, e);
+            const double D2 = coeff(1, q, e);
+            const double D3 = coeff(2, q, e);
+            // detJ J^{-1} DJ = adj(J) DJ
+            y(i11,q,e) = w_detJ * (D1*A11*J11 + D2*A12*J21 + D3*A13*J31); // 1,1
+            y(i12,q,e) = w_detJ * (D1*A11*J12 + D2*A12*J22 + D3*A13*J32); // 1,2
+            y(i13,q,e) = w_detJ * (D1*A11*J13 + D2*A12*J23 + D3*A13*J33); // 1,3
+            y(i21,q,e) = w_detJ * (D1*A21*J11 + D2*A22*J21 + D3*A23*J31); // 2,1
+            y(i22,q,e) = w_detJ * (D1*A21*J12 + D2*A22*J22 + D3*A23*J32); // 2,2
+            y(i23,q,e) = w_detJ * (D1*A21*J13 + D2*A22*J23 + D3*A23*J33); // 2,3
+            y(i31,q,e) = w_detJ * (D1*A31*J11 + D2*A32*J21 + D3*A33*J31); // 3,1
+            y(i32,q,e) = w_detJ * (D1*A31*J12 + D2*A32*J22 + D3*A33*J32); // 3,2
+            y(i33,q,e) = w_detJ * (D1*A31*J13 + D2*A32*J23 + D3*A33*J33); // 3,3
+         }
+      }
+   });
+}
+
+// PA H(curl) x H(div) mass assemble 2D kernel, with factor
+// dF^{-1} C dF for a vector or matrix coefficient C.
+// If transpose, use dF^T C dF^{-T} for H(div) x H(curl).
+void PAHcurlHdivSetup2D(const int Q1D,
+                        const int coeffDim,
+                        const int NE,
+                        const bool transpose,
+                        const Array<double> &_w,
+                        const Vector &j,
+                        Vector &_coeff,
+                        Vector &op)
+{
+   const int NQ = Q1D*Q1D;
+   const bool symmetric = (coeffDim != 4);
+   auto W = _w.Read();
+   auto J = Reshape(j.Read(), NQ, 2, 2, NE);
+   auto coeff = Reshape(_coeff.Read(), coeffDim, NQ, NE);
+   auto y = Reshape(op.Write(), 4, NQ, NE);
+
+   const int i11 = 0;
+   const int i12 = transpose ? 2 : 1;
+   const int i21 = transpose ? 1 : 2;
+   const int i22 = 3;
+
+   MFEM_FORALL(e, NE,
+   {
+      for (int q = 0; q < NQ; ++q)
+      {
+         const double J11 = J(q,0,0,e);
+         const double J21 = J(q,1,0,e);
+         const double J12 = J(q,0,1,e);
+         const double J22 = J(q,1,1,e);
+         const double w_detJ = W[q] / (J11*J22) - (J21*J12);
+
+         if (coeffDim == 3 || coeffDim == 4) // Matrix coefficient version
+         {
+            // First compute entries of R = MJ
+            const double M11 = coeff(i11, q, e);
+            const double M12 = (!symmetric) ? coeff(i12, q, e) : coeff(1, q, e);
+            const double M21 = (!symmetric) ? coeff(i21, q, e) : M12;
+            const double M22 = (!symmetric) ? coeff(i22, q, e) : coeff(2, q, e);
+
+            const double R11 = M11*J11 + M12*J21;
+            const double R12 = M11*J12 + M12*J22;
+            const double R21 = M21*J11 + M22*J21;
+            const double R22 = M21*J12 + M22*J22;
+
+            // Now set y to J^{-1} R
+            y(i11,q,e) = w_detJ * ( J22*R11 - J12*R21); // 1,1
+            y(i12,q,e) = w_detJ * ( J22*R12 - J12*R22); // 1,2
+            y(i21,q,e) = w_detJ * (-J21*R11 + J11*R21); // 2,1
+            y(i22,q,e) = w_detJ * (-J21*R12 + J11*R22); // 2,2
+         }
+         else if (coeffDim == 2) // Vector coefficient version
+         {
+            const double D1 = coeff(0, q, e);
+            const double D2 = coeff(1, q, e);
+            const double R11 = D1*J11;
+            const double R12 = D1*J12;
+            const double R21 = D2*J21;
+            const double R22 = D2*J22;
+            y(i11,q,e) = w_detJ * ( J22*R11 - J12*R21); // 1,1
+            y(i12,q,e) = w_detJ * ( J22*R12 - J12*R22); // 1,2
+            y(i21,q,e) = w_detJ * (-J21*R11 + J11*R21); // 2,1
+            y(i22,q,e) = w_detJ * (-J21*R12 + J11*R22); // 2,2
+         }
+      }
+   });
+}
+
+// Mass operator for H(curl) and H(div) functions, using Piola transformations
+// u = dF^{-T} \hat{u} in H(curl), v = (1 / det dF) dF \hat{v} in H(div).
+void PAHcurlHdivMassApply3D(const int D1D,
+                            const int D1Dtest,
+                            const int Q1D,
+                            const int NE,
+                            const bool scalarCoeff,
+                            const bool trialHcurl,
+                            const Array<double> &_Bo,
+                            const Array<double> &_Bc,
+                            const Array<double> &_Bot,
+                            const Array<double> &_Bct,
+                            const Vector &_op,
+                            const Vector &_x,
+                            Vector &_y)
+{
+   constexpr static int MAX_D1D = HCURL_MAX_D1D;
+   constexpr static int MAX_Q1D = HCURL_MAX_Q1D;
+
+   MFEM_VERIFY(D1D <= MAX_D1D, "Error: D1D > MAX_D1D");
+   MFEM_VERIFY(Q1D <= MAX_Q1D, "Error: Q1D > MAX_Q1D");
+   constexpr static int VDIM = 3;
+
+   auto Bo = Reshape(_Bo.Read(), Q1D, D1D-1);
+   auto Bc = Reshape(_Bc.Read(), Q1D, D1D);
+   auto Bot = Reshape(_Bot.Read(), D1Dtest-1, Q1D);
+   auto Bct = Reshape(_Bct.Read(), D1Dtest, Q1D);
+   auto op = Reshape(_op.Read(), scalarCoeff ? 1 : 9, Q1D, Q1D, Q1D, NE);
+   auto x = Reshape(_x.Read(), 3*(D1D-1)*D1D*(trialHcurl ? D1D : D1D-1), NE);
+   auto y = Reshape(_y.ReadWrite(), 3*(D1Dtest-1)*D1Dtest*
+                    (trialHcurl ? D1Dtest-1 : D1Dtest), NE);
+
+   MFEM_FORALL(e, NE,
+   {
+      double mass[MAX_Q1D][MAX_Q1D][MAX_Q1D][VDIM];
+
+      for (int qz = 0; qz < Q1D; ++qz)
+      {
+         for (int qy = 0; qy < Q1D; ++qy)
+         {
+            for (int qx = 0; qx < Q1D; ++qx)
+            {
+               for (int c = 0; c < VDIM; ++c)
+               {
+                  mass[qz][qy][qx][c] = 0.0;
+               }
+            }
+         }
+      }
+
+      int osc = 0;
+      for (int c = 0; c < VDIM; ++c)  // loop over x, y, z trial components
+      {
+         const int D1Dz = trialHcurl ? ((c == 2) ? D1D - 1 : D1D) :
+                          ((c == 2) ? D1D : D1D - 1);
+         const int D1Dy = trialHcurl ? ((c == 1) ? D1D - 1 : D1D) :
+                          ((c == 1) ? D1D : D1D - 1);
+         const int D1Dx = trialHcurl ? ((c == 0) ? D1D - 1 : D1D) :
+                          ((c == 0) ? D1D : D1D - 1);
+
+         for (int dz = 0; dz < D1Dz; ++dz)
+         {
+            double massXY[MAX_Q1D][MAX_Q1D];
+            for (int qy = 0; qy < Q1D; ++qy)
+            {
+               for (int qx = 0; qx < Q1D; ++qx)
+               {
+                  massXY[qy][qx] = 0.0;
+               }
+            }
+
+            for (int dy = 0; dy < D1Dy; ++dy)
+            {
+               double massX[MAX_Q1D];
+               for (int qx = 0; qx < Q1D; ++qx)
+               {
+                  massX[qx] = 0.0;
+               }
+
+               for (int dx = 0; dx < D1Dx; ++dx)
+               {
+                  const double t = x(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc, e);
+                  for (int qx = 0; qx < Q1D; ++qx)
+                  {
+                     massX[qx] += t * (trialHcurl ? ((c == 0) ? Bo(qx,dx) : Bc(qx,dx)) :
+                                       ((c == 0) ? Bc(qx,dx) : Bo(qx,dx)));
+                  }
+               }
+
+               for (int qy = 0; qy < Q1D; ++qy)
+               {
+                  const double wy = trialHcurl ? ((c == 1) ? Bo(qy,dy) : Bc(qy,dy)) :
+                                    ((c == 1) ? Bc(qy,dy) : Bo(qy,dy));
+                  for (int qx = 0; qx < Q1D; ++qx)
+                  {
+                     const double wx = massX[qx];
+                     massXY[qy][qx] += wx * wy;
+                  }
+               }
+            }
+
+            for (int qz = 0; qz < Q1D; ++qz)
+            {
+               const double wz = trialHcurl ? ((c == 2) ? Bo(qz,dz) : Bc(qz,dz)) :
+                                 ((c == 2) ? Bc(qz,dz) : Bo(qz,dz));
+               for (int qy = 0; qy < Q1D; ++qy)
+               {
+                  for (int qx = 0; qx < Q1D; ++qx)
+                  {
+                     mass[qz][qy][qx][c] += massXY[qy][qx] * wz;
+                  }
+               }
+            }
+         }
+
+         osc += D1Dx * D1Dy * D1Dz;
+      }  // loop (c) over components
+
+      // Apply D operator.
+      for (int qz = 0; qz < Q1D; ++qz)
+      {
+         for (int qy = 0; qy < Q1D; ++qy)
+         {
+            for (int qx = 0; qx < Q1D; ++qx)
+            {
+               const double O11 = op(0,qx,qy,qz,e);
+               const double O12 = scalarCoeff ? 0.0 : op(1,qx,qy,qz,e);
+               const double O13 = scalarCoeff ? 0.0 : op(2,qx,qy,qz,e);
+               const double O21 = scalarCoeff ? 0.0 : op(3,qx,qy,qz,e);
+               const double O22 = scalarCoeff ? O11 : op(4,qx,qy,qz,e);
+               const double O23 = scalarCoeff ? 0.0 : op(5,qx,qy,qz,e);
+               const double O31 = scalarCoeff ? 0.0 : op(6,qx,qy,qz,e);
+               const double O32 = scalarCoeff ? 0.0 : op(7,qx,qy,qz,e);
+               const double O33 = scalarCoeff ? O11 : op(8,qx,qy,qz,e);
+               const double massX = mass[qz][qy][qx][0];
+               const double massY = mass[qz][qy][qx][1];
+               const double massZ = mass[qz][qy][qx][2];
+               mass[qz][qy][qx][0] = (O11*massX)+(O12*massY)+(O13*massZ);
+               mass[qz][qy][qx][1] = (O21*massX)+(O22*massY)+(O23*massZ);
+               mass[qz][qy][qx][2] = (O31*massX)+(O32*massY)+(O33*massZ);
+            }
+         }
+      }
+
+      for (int qz = 0; qz < Q1D; ++qz)
+      {
+         double massXY[HDIV_MAX_D1D][HDIV_MAX_D1D];
+
+         osc = 0;
+         for (int c = 0; c < VDIM; ++c)  // loop over x, y, z test components
+         {
+            const int D1Dz = trialHcurl ? ((c == 2) ? D1Dtest : D1Dtest - 1) :
+                             ((c == 2) ? D1Dtest - 1 : D1Dtest);
+            const int D1Dy = trialHcurl ? ((c == 1) ? D1Dtest : D1Dtest - 1) :
+                             ((c == 1) ? D1Dtest - 1 : D1Dtest);
+            const int D1Dx = trialHcurl ? ((c == 0) ? D1Dtest : D1Dtest - 1) :
+                             ((c == 0) ? D1Dtest - 1 : D1Dtest);
+
+            for (int dy = 0; dy < D1Dy; ++dy)
+            {
+               for (int dx = 0; dx < D1Dx; ++dx)
+               {
+                  massXY[dy][dx] = 0.0;
+               }
+            }
+            for (int qy = 0; qy < Q1D; ++qy)
+            {
+               double massX[HDIV_MAX_D1D];
+               for (int dx = 0; dx < D1Dx; ++dx)
+               {
+                  massX[dx] = 0.0;
+               }
+               for (int qx = 0; qx < Q1D; ++qx)
+               {
+                  for (int dx = 0; dx < D1Dx; ++dx)
+                  {
+                     massX[dx] += mass[qz][qy][qx][c] * (trialHcurl ?
+                                                         ((c == 0) ? Bct(dx,qx) : Bot(dx,qx)) :
+                                                         ((c == 0) ? Bot(dx,qx) : Bct(dx,qx)));
+                  }
+               }
+               for (int dy = 0; dy < D1Dy; ++dy)
+               {
+                  const double wy = trialHcurl ? ((c == 1) ? Bct(dy,qy) : Bot(dy,qy)) :
+                                    ((c == 1) ? Bot(dy,qy) : Bct(dy,qy));
+                  for (int dx = 0; dx < D1Dx; ++dx)
+                  {
+                     massXY[dy][dx] += massX[dx] * wy;
+                  }
+               }
+            }
+
+            for (int dz = 0; dz < D1Dz; ++dz)
+            {
+               const double wz = trialHcurl ? ((c == 2) ? Bct(dz,qz) : Bot(dz,qz)) :
+                                 ((c == 2) ? Bot(dz,qz) : Bct(dz,qz));
+               for (int dy = 0; dy < D1Dy; ++dy)
+               {
+                  for (int dx = 0; dx < D1Dx; ++dx)
+                  {
+                     y(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc, e) +=
+                        massXY[dy][dx] * wz;
+                  }
+               }
+            }
+
+            osc += D1Dx * D1Dy * D1Dz;
+         }  // loop c
+      }  // loop qz
+   }); // end of element loop
+}
+
+// Mass operator for H(curl) and H(div) functions, using Piola transformations
+// u = dF^{-T} \hat{u} in H(curl), v = (1 / det dF) dF \hat{v} in H(div).
+void PAHcurlHdivMassApply2D(const int D1D,
+                            const int D1Dtest,
+                            const int Q1D,
+                            const int NE,
+                            const bool scalarCoeff,
+                            const bool trialHcurl,
+                            const Array<double> &_Bo,
+                            const Array<double> &_Bc,
+                            const Array<double> &_Bot,
+                            const Array<double> &_Bct,
+                            const Vector &_op,
+                            const Vector &_x,
+                            Vector &_y)
+{
+   constexpr static int MAX_D1D = HCURL_MAX_D1D;
+   constexpr static int MAX_Q1D = HCURL_MAX_Q1D;
+
+   MFEM_VERIFY(D1D <= MAX_D1D, "Error: D1D > MAX_D1D");
+   MFEM_VERIFY(Q1D <= MAX_Q1D, "Error: Q1D > MAX_Q1D");
+   constexpr static int VDIM = 2;
+
+   auto Bo = Reshape(_Bo.Read(), Q1D, D1D-1);
+   auto Bc = Reshape(_Bc.Read(), Q1D, D1D);
+   auto Bot = Reshape(_Bot.Read(), D1Dtest-1, Q1D);
+   auto Bct = Reshape(_Bct.Read(), D1Dtest, Q1D);
+   auto op = Reshape(_op.Read(), scalarCoeff ? 1 : 4, Q1D, Q1D, NE);
+   auto x = Reshape(_x.Read(), 2*(D1D-1)*D1D, NE);
+   auto y = Reshape(_y.ReadWrite(), 2*(D1Dtest-1)*D1Dtest, NE);
+
+   MFEM_FORALL(e, NE,
+   {
+      double mass[MAX_Q1D][MAX_Q1D][VDIM];
+
+      for (int qy = 0; qy < Q1D; ++qy)
+      {
+         for (int qx = 0; qx < Q1D; ++qx)
+         {
+            for (int c = 0; c < VDIM; ++c)
+            {
+               mass[qy][qx][c] = 0.0;
+            }
+         }
+      }
+
+      int osc = 0;
+      for (int c = 0; c < VDIM; ++c)  // loop over x, y trial components
+      {
+         const int D1Dy = trialHcurl ? ((c == 1) ? D1D - 1 : D1D) :
+                          ((c == 1) ? D1D : D1D - 1);
+         const int D1Dx = trialHcurl ? ((c == 0) ? D1D - 1 : D1D) :
+                          ((c == 0) ? D1D : D1D - 1);
+
+         for (int dy = 0; dy < D1Dy; ++dy)
+         {
+            double massX[MAX_Q1D];
+            for (int qx = 0; qx < Q1D; ++qx)
+            {
+               massX[qx] = 0.0;
+            }
+
+            for (int dx = 0; dx < D1Dx; ++dx)
+            {
+               const double t = x(dx + (dy * D1Dx) + osc, e);
+               for (int qx = 0; qx < Q1D; ++qx)
+               {
+                  massX[qx] += t * (trialHcurl ? ((c == 0) ? Bo(qx,dx) : Bc(qx,dx)) :
+                                    ((c == 0) ? Bc(qx,dx) : Bo(qx,dx)));
+               }
+            }
+
+            for (int qy = 0; qy < Q1D; ++qy)
+            {
+               const double wy = trialHcurl ? ((c == 1) ? Bo(qy,dy) : Bc(qy,dy)) :
+                                 ((c == 1) ? Bc(qy,dy) : Bo(qy,dy));
+               for (int qx = 0; qx < Q1D; ++qx)
+               {
+                  mass[qy][qx][c] += massX[qx] * wy;
+               }
+            }
+         }
+
+         osc += D1Dx * D1Dy;
+      }  // loop (c) over components
+
+      // Apply D operator.
+      for (int qy = 0; qy < Q1D; ++qy)
+      {
+         for (int qx = 0; qx < Q1D; ++qx)
+         {
+            const double O11 = op(0,qx,qy,e);
+            const double O12 = scalarCoeff ? 0.0 : op(1,qx,qy,e);
+            const double O21 = scalarCoeff ? 0.0 : op(2,qx,qy,e);
+            const double O22 = scalarCoeff ? O11 : op(3,qx,qy,e);
+            const double massX = mass[qy][qx][0];
+            const double massY = mass[qy][qx][1];
+            mass[qy][qx][0] = (O11*massX)+(O12*massY);
+            mass[qy][qx][1] = (O21*massX)+(O22*massY);
+         }
+      }
+
+      osc = 0;
+      for (int c = 0; c < VDIM; ++c)  // loop over x, y test components
+      {
+         const int D1Dy = trialHcurl ? ((c == 1) ? D1Dtest : D1Dtest - 1) :
+                          ((c == 1) ? D1Dtest - 1 : D1Dtest);
+         const int D1Dx = trialHcurl ? ((c == 0) ? D1Dtest : D1Dtest - 1) :
+                          ((c == 0) ? D1Dtest - 1 : D1Dtest);
+
+         for (int qy = 0; qy < Q1D; ++qy)
+         {
+            double massX[HDIV_MAX_D1D];
+            for (int dx = 0; dx < D1Dx; ++dx)
+            {
+               massX[dx] = 0.0;
+            }
+            for (int qx = 0; qx < Q1D; ++qx)
+            {
+               for (int dx = 0; dx < D1Dx; ++dx)
+               {
+                  massX[dx] += mass[qy][qx][c] * (trialHcurl ?
+                                                  ((c == 0) ? Bct(dx,qx) : Bot(dx,qx)) :
+                                                  ((c == 0) ? Bot(dx,qx) : Bct(dx,qx)));
+               }
+            }
+            for (int dy = 0; dy < D1Dy; ++dy)
+            {
+               const double wy = trialHcurl ? ((c == 1) ? Bct(dy,qy) : Bot(dy,qy)) :
+                                 ((c == 1) ? Bot(dy,qy) : Bct(dy,qy));
+               for (int dx = 0; dx < D1Dx; ++dx)
+               {
+                  y(dx + (dy * D1Dx) + osc, e) += massX[dx] * wy;
+               }
+            }
+         }
+
+         osc += D1Dx * D1Dy;
+      }  // loop c
+   }); // end of element loop
+}
+
 void VectorFEMassIntegrator::AssemblePA(const FiniteElementSpace &fes)
 {
-   // Assumes tensor-product elements
-   Mesh *mesh = fes.GetMesh();
-   const FiniteElement *fel = fes.GetFE(0);
+   AssemblePA(fes, fes);
+}
 
-   const VectorTensorFiniteElement *el =
-      dynamic_cast<const VectorTensorFiniteElement*>(fel);
-   MFEM_VERIFY(el != NULL, "Only VectorTensorFiniteElement is supported!");
+void VectorFEMassIntegrator::AssemblePA(const FiniteElementSpace &trial_fes,
+                                        const FiniteElementSpace &test_fes)
+{
+   // Assumes tensor-product elements
+   Mesh *mesh = trial_fes.GetMesh();
+
+   const FiniteElement *trial_fel = trial_fes.GetFE(0);
+   const VectorTensorFiniteElement *trial_el =
+      dynamic_cast<const VectorTensorFiniteElement*>(trial_fel);
+   MFEM_VERIFY(trial_el != NULL, "Only VectorTensorFiniteElement is supported!");
+
+   const FiniteElement *test_fel = test_fes.GetFE(0);
+   const VectorTensorFiniteElement *test_el =
+      dynamic_cast<const VectorTensorFiniteElement*>(test_fel);
+   MFEM_VERIFY(test_el != NULL, "Only VectorTensorFiniteElement is supported!");
 
    const IntegrationRule *ir
-      = IntRule ? IntRule : &MassIntegrator::GetRule(*el, *el,
+      = IntRule ? IntRule : &MassIntegrator::GetRule(*trial_el, *trial_el,
                                                      *mesh->GetElementTransformation(0));
-   const int dims = el->GetDim();
+   const int dims = trial_el->GetDim();
    MFEM_VERIFY(dims == 2 || dims == 3, "");
 
    const int symmDims = (dims * (dims + 1)) / 2; // 1x1: 1, 2x2: 3, 3x3: 6
@@ -164,36 +721,99 @@ void VectorFEMassIntegrator::AssemblePA(const FiniteElementSpace &fes)
    dim = mesh->Dimension();
    MFEM_VERIFY(dim == 2 || dim == 3, "");
 
-   ne = fes.GetNE();
+   ne = trial_fes.GetNE();
+   MFEM_VERIFY(ne == test_fes.GetNE(),
+               "Different meshes for test and trial spaces");
    geom = mesh->GetGeometricFactors(*ir, GeometricFactors::JACOBIANS);
-   mapsC = &el->GetDofToQuad(*ir, DofToQuad::TENSOR);
-   mapsO = &el->GetDofToQuadOpen(*ir, DofToQuad::TENSOR);
+   mapsC = &trial_el->GetDofToQuad(*ir, DofToQuad::TENSOR);
+   mapsO = &trial_el->GetDofToQuadOpen(*ir, DofToQuad::TENSOR);
    dofs1D = mapsC->ndof;
    quad1D = mapsC->nqpt;
 
+   mapsCtest = &test_el->GetDofToQuad(*ir, DofToQuad::TENSOR);
+   mapsOtest = &test_el->GetDofToQuadOpen(*ir, DofToQuad::TENSOR);
+   dofs1Dtest = mapsCtest->ndof;
+
    MFEM_VERIFY(dofs1D == mapsO->ndof + 1 && quad1D == mapsO->nqpt, "");
 
-   pa_data.SetSize(symmDims * nq * ne, Device::GetMemoryType());
+   trial_fetype = trial_el->GetDerivType();
+   test_fetype = test_el->GetDerivType();
 
-   const int coeffDim = VQ ? VQ->GetVDim() : 1;
+   const int MQsymmDim = MQ ? (MQ->GetWidth() * (MQ->GetWidth() + 1)) / 2 : 0;
+   const int MQfullDim = MQ ? (MQ->GetHeight() * MQ->GetWidth()) : 0;
+   const int MQdim = MQ ? (MQ->IsSymmetric() ? MQsymmDim : MQfullDim) : 0;
+   const int coeffDim = MQ ? MQdim : (VQ ? VQ->GetVDim() : 1);
+
+   symmetric = MQ ? MQ->IsSymmetric() : true;
+
+   if ((trial_fetype == mfem::FiniteElement::CURL &&
+        test_fetype == mfem::FiniteElement::DIV) ||
+       (trial_fetype == mfem::FiniteElement::DIV &&
+        test_fetype == mfem::FiniteElement::CURL))
+      pa_data.SetSize((coeffDim == 1 ? 1 : dim*dim) * nq * ne,
+                      Device::GetMemoryType());
+   else
+      pa_data.SetSize((symmetric ? symmDims : MQfullDim) * nq * ne,
+                      Device::GetMemoryType());
 
    Vector coeff(coeffDim * ne * nq);
    coeff = 1.0;
    auto coeffh = Reshape(coeff.HostWrite(), coeffDim, nq, ne);
-   if (Q || VQ)
+   if (Q || VQ || MQ)
    {
       Vector D(VQ ? coeffDim : 0);
+      DenseMatrix M;
+      Vector Msymm;
+      if (MQ)
+      {
+         if (symmetric)
+         {
+            Msymm.SetSize(MQsymmDim);
+         }
+         else
+         {
+            M.SetSize(dim);
+         }
+      }
+
       if (VQ)
       {
          MFEM_VERIFY(coeffDim == dim, "");
       }
+      if (MQ)
+      {
+         MFEM_VERIFY(coeffDim == MQdim, "");
+         MFEM_VERIFY(MQ->GetHeight() == dim && MQ->GetWidth() == dim, "");
+      }
 
       for (int e=0; e<ne; ++e)
       {
          ElementTransformation *tr = mesh->GetElementTransformation(e);
          for (int p=0; p<nq; ++p)
          {
-            if (VQ)
+            if (MQ)
+            {
+               if (MQ->IsSymmetric())
+               {
+                  MQ->EvalSymmetric(Msymm, *tr, ir->IntPoint(p));
+
+                  for (int i=0; i<MQsymmDim; ++i)
+                  {
+                     coeffh(i, p, e) = Msymm[i];
+                  }
+               }
+               else
+               {
+                  MQ->Eval(M, *tr, ir->IntPoint(p));
+
+                  for (int i=0; i<dim; ++i)
+                     for (int j=0; j<dim; ++j)
+                     {
+                        coeffh(j+(i*dim), p, e) = M(i,j);
+                     }
+               }
+            }
+            else if (VQ)
             {
                VQ->Eval(D, *tr, ir->IntPoint(p));
                for (int i=0; i<coeffDim; ++i)
@@ -209,28 +829,52 @@ void VectorFEMassIntegrator::AssemblePA(const FiniteElementSpace &fes)
       }
    }
 
-   fetype = el->GetDerivType();
-
-   if (el->GetDerivType() == mfem::FiniteElement::CURL && dim == 3)
+   if (trial_fetype == mfem::FiniteElement::CURL && test_fetype == trial_fetype
+       && dim == 3)
    {
       PAHcurlSetup3D(quad1D, coeffDim, ne, ir->GetWeights(), geom->J,
                      coeff, pa_data);
    }
-   else if (el->GetDerivType() == mfem::FiniteElement::CURL && dim == 2)
+   else if (trial_fetype == mfem::FiniteElement::CURL
+            && test_fetype == trial_fetype && dim == 2)
    {
       PAHcurlSetup2D(quad1D, coeffDim, ne, ir->GetWeights(), geom->J,
                      coeff, pa_data);
    }
-   else if (el->GetDerivType() == mfem::FiniteElement::DIV && dim == 3)
+   else if (trial_fetype == mfem::FiniteElement::DIV
+            && test_fetype == trial_fetype && dim == 3)
    {
       PAHdivSetup3D(quad1D, ne, ir->GetWeights(), geom->J,
                     coeff, pa_data);
    }
-   else if (el->GetDerivType() == mfem::FiniteElement::DIV && dim == 2)
+   else if (trial_fetype == mfem::FiniteElement::DIV
+            && test_fetype == trial_fetype && dim == 2)
    {
       PAHdivSetup2D(quad1D, ne, ir->GetWeights(), geom->J,
                     coeff, pa_data);
    }
+   else if (((trial_fetype == mfem::FiniteElement::CURL &&
+              test_fetype == mfem::FiniteElement::DIV) ||
+             (trial_fetype == mfem::FiniteElement::DIV &&
+              test_fetype == mfem::FiniteElement::CURL)) &&
+            test_fel->GetOrder() == trial_fel->GetOrder())
+   {
+      if (coeffDim == 1)
+      {
+         PAHcurlL2Setup(nq, coeffDim, ne, ir->GetWeights(), coeff, pa_data);
+      }
+      else
+      {
+         const bool tr = (trial_fetype == mfem::FiniteElement::DIV &&
+                          test_fetype == mfem::FiniteElement::CURL);
+         if (dim == 3)
+            PAHcurlHdivSetup3D(quad1D, coeffDim, ne, tr, ir->GetWeights(),
+                               geom->J, coeff, pa_data);
+         else
+            PAHcurlHdivSetup2D(quad1D, coeffDim, ne, tr, ir->GetWeights(),
+                               geom->J, coeff, pa_data);
+      }
+   }
    else
    {
       MFEM_ABORT("Unknown kernel.");
@@ -241,12 +885,13 @@ void VectorFEMassIntegrator::AssembleDiagonalPA(Vector& diag)
 {
    if (dim == 3)
    {
-      if (fetype == mfem::FiniteElement::CURL)
+      if (trial_fetype == mfem::FiniteElement::CURL && test_fetype == trial_fetype)
       {
-         PAHcurlMassAssembleDiagonal3D(dofs1D, quad1D, ne,
+         PAHcurlMassAssembleDiagonal3D(dofs1D, quad1D, ne, symmetric,
                                        mapsO->B, mapsC->B, pa_data, diag);
       }
-      else if (fetype == mfem::FiniteElement::DIV)
+      else if (trial_fetype == mfem::FiniteElement::DIV &&
+               test_fetype == trial_fetype)
       {
          PAHdivMassAssembleDiagonal3D(dofs1D, quad1D, ne,
                                       mapsO->B, mapsC->B, pa_data, diag);
@@ -258,12 +903,13 @@ void VectorFEMassIntegrator::AssembleDiagonalPA(Vector& diag)
    }
    else
    {
-      if (fetype == mfem::FiniteElement::CURL)
+      if (trial_fetype == mfem::FiniteElement::CURL && test_fetype == trial_fetype)
       {
-         PAHcurlMassAssembleDiagonal2D(dofs1D, quad1D, ne,
+         PAHcurlMassAssembleDiagonal2D(dofs1D, quad1D, ne, symmetric,
                                        mapsO->B, mapsC->B, pa_data, diag);
       }
-      else if (fetype == mfem::FiniteElement::DIV)
+      else if (trial_fetype == mfem::FiniteElement::DIV &&
+               test_fetype == trial_fetype)
       {
          PAHdivMassAssembleDiagonal2D(dofs1D, quad1D, ne,
                                       mapsO->B, mapsC->B, pa_data, diag);
@@ -279,16 +925,33 @@ void VectorFEMassIntegrator::AddMultPA(const Vector &x, Vector &y) const
 {
    if (dim == 3)
    {
-      if (fetype == mfem::FiniteElement::CURL)
+      if (trial_fetype == mfem::FiniteElement::CURL && test_fetype == trial_fetype)
       {
-         PAHcurlMassApply3D(dofs1D, quad1D, ne, mapsO->B, mapsC->B, mapsO->Bt,
-                            mapsC->Bt, pa_data, x, y);
+         PAHcurlMassApply3D(dofs1D, quad1D, ne, symmetric, mapsO->B, mapsC->B,
+                            mapsO->Bt, mapsC->Bt, pa_data, x, y);
       }
-      else if (fetype == mfem::FiniteElement::DIV)
+      else if (trial_fetype == mfem::FiniteElement::DIV &&
+               test_fetype == trial_fetype)
       {
          PAHdivMassApply3D(dofs1D, quad1D, ne, mapsO->B, mapsC->B, mapsO->Bt,
                            mapsC->Bt, pa_data, x, y);
       }
+      else if (trial_fetype == mfem::FiniteElement::CURL &&
+               test_fetype == mfem::FiniteElement::DIV)
+      {
+         const bool scalarCoeff = !(VQ || MQ);
+         PAHcurlHdivMassApply3D(dofs1D, dofs1Dtest, quad1D, ne, scalarCoeff,
+                                true, mapsO->B, mapsC->B, mapsOtest->Bt,
+                                mapsCtest->Bt, pa_data, x, y);
+      }
+      else if (trial_fetype == mfem::FiniteElement::DIV &&
+               test_fetype == mfem::FiniteElement::CURL)
+      {
+         const bool scalarCoeff = !(VQ || MQ);
+         PAHcurlHdivMassApply3D(dofs1D, dofs1Dtest, quad1D, ne, scalarCoeff,
+                                false, mapsO->B, mapsC->B, mapsOtest->Bt,
+                                mapsCtest->Bt, pa_data, x, y);
+      }
       else
       {
          MFEM_ABORT("Unknown kernel.");
@@ -296,16 +959,28 @@ void VectorFEMassIntegrator::AddMultPA(const Vector &x, Vector &y) const
    }
    else
    {
-      if (fetype == mfem::FiniteElement::CURL)
+      if (trial_fetype == mfem::FiniteElement::CURL && test_fetype == trial_fetype)
       {
-         PAHcurlMassApply2D(dofs1D, quad1D, ne, mapsO->B, mapsC->B, mapsO->Bt,
-                            mapsC->Bt, pa_data, x, y);
+         PAHcurlMassApply2D(dofs1D, quad1D, ne, symmetric, mapsO->B, mapsC->B,
+                            mapsO->Bt, mapsC->Bt, pa_data, x, y);
       }
-      else if (fetype == mfem::FiniteElement::DIV)
+      else if (trial_fetype == mfem::FiniteElement::DIV &&
+               test_fetype == trial_fetype)
       {
          PAHdivMassApply2D(dofs1D, quad1D, ne, mapsO->B, mapsC->B, mapsO->Bt,
                            mapsC->Bt, pa_data, x, y);
       }
+      else if ((trial_fetype == mfem::FiniteElement::CURL &&
+                test_fetype == mfem::FiniteElement::DIV) ||
+               (trial_fetype == mfem::FiniteElement::DIV &&
+                test_fetype == mfem::FiniteElement::CURL))
+      {
+         const bool scalarCoeff = !(VQ || MQ);
+         const bool trialHcurl = (trial_fetype == mfem::FiniteElement::CURL);
+         PAHcurlHdivMassApply2D(dofs1D, dofs1Dtest, quad1D, ne, scalarCoeff,
+                                trialHcurl, mapsO->B, mapsC->B, mapsOtest->Bt,
+                                mapsCtest->Bt, pa_data, x, y);
+      }
       else
       {
          MFEM_ABORT("Unknown kernel.");

fem/coefficient.cpp

@@ -319,6 +319,31 @@ void MatrixFunctionCoefficient::Eval(DenseMatrix &K, ElementTransformation &T,
    }
 }
 
+void MatrixFunctionCoefficient::EvalSymmetric(Vector &K,
+                                              ElementTransformation &T,
+                                              const IntegrationPoint &ip)
+{
+   MFEM_VERIFY(symmetric && height == width && height < 4 && SymmFunction,
+               "MatrixFunctionCoefficient is not symmetric");
+
+   double x[3];
+   Vector transip(x, 3);
+
+   T.Transform(ip, transip);
+
+   K.SetSize((width * (width + 1)) / 2); // 1x1: 1, 2x2: 3, 3x3: 6
+
+   if (SymmFunction)
+   {
+      (*SymmFunction)(transip, K);
+   }
+
+   if (Q)
+   {
+      K *= Q->Eval(T, ip, GetTime());
+   }
+}
+
 MatrixArrayCoefficient::MatrixArrayCoefficient (int dim)
    : MatrixCoefficient (dim)
 {

fem/coefficient.hpp

@@ -695,13 +695,16 @@ class MatrixCoefficient
 protected:
    int height, width;
    double time;
+   bool symmetric;
 
 public:
    /// Construct a dim x dim matrix coefficient.
-   explicit MatrixCoefficient(int dim) { height = width = dim; time = 0.; }
+   explicit MatrixCoefficient(int dim, bool symm=false)
+   { height = width = dim; time = 0.; symmetric = symm; }
 
    /// Construct a h x w matrix coefficient.
-   MatrixCoefficient(int h, int w) : height(h), width(w), time(0.) { }
+   MatrixCoefficient(int h, int w, bool symm=false) :
+      height(h), width(w), time(0.), symmetric(symm) { }
 
    /// Set the time for time dependent coefficients
    void SetTime(double t) { time = t; }
@@ -718,6 +721,9 @@ public:
    /// For backward compatibility get the width of the matrix.
    int GetVDim() const { return width; }
 
+   void SetSymmetric(bool s) { symmetric = s; }
+   bool IsSymmetric() const { return symmetric; }
+
    /** @brief Evaluate the matrix coefficient in the element described by @a T
        at the point @a ip, storing the result in @a K. */
    /** @note When this method is called, the caller must make sure that the
@@ -726,6 +732,15 @@ public:
    virtual void Eval(DenseMatrix &K, ElementTransformation &T,
                      const IntegrationPoint &ip) = 0;
 
+   /** @brief Evaluate the upper triangular entries of the matrix coefficient
+       in the symmetric case, similarly to Eval. Matrix entry (i,j) is stored
+       in K[j - i + os_i] for 0 <= i <= j < width, os_0 = 0,
+       os_{i+1} = os_i + width - i. That is, K = {M(0,0), ..., M(0,w-1),
+       M(1,1), ..., M(1,w-1), ..., M(w-1,w-1) with w = width. */
+   virtual void EvalSymmetric(Vector &K, ElementTransformation &T,
+                              const IntegrationPoint &ip)
+   { mfem_error("MatrixCoefficient::EvalSymmetric"); }
+
    virtual ~MatrixCoefficient() { }
 };
 
@@ -753,6 +768,7 @@ class MatrixFunctionCoefficient : public MatrixCoefficient
 {
 private:
    void (*Function)(const Vector &, DenseMatrix &);
+   void (*SymmFunction)(const Vector &, Vector &);
    void (*TDFunction)(const Vector &, double, DenseMatrix &);
    Coefficient *Q;
    DenseMatrix mat;
@@ -790,10 +806,26 @@ public:
       mat.SetSize(0);
    }
 
+   /// Construct a symmetric square matrix coefficient from a C-function
+   /// defining a vector function used by EvalSymmetric
+   MatrixFunctionCoefficient(int dim, void (*F)(const Vector &, Vector &),
+                             Coefficient *q = NULL)
+      : MatrixCoefficient(dim, true), Q(q)
+   {
+      SymmFunction = F;
+      Function = NULL;
+      TDFunction = NULL;
+      mat.SetSize(0);
+   }
+
    /// Evaluate the matrix coefficient at @a ip.
    virtual void Eval(DenseMatrix &K, ElementTransformation &T,
                      const IntegrationPoint &ip);
 
+   /// Evaluate the symmetric matrix coefficient at @a ip.
+   virtual void EvalSymmetric(Vector &K, ElementTransformation &T,
+                              const IntegrationPoint &ip);
+
    virtual ~MatrixFunctionCoefficient() { }
 };
 

fem/complex_fem.cpp

@@ -10,6 +10,7 @@
 // CONTRIBUTING.md for details.
 
 #include "complex_fem.hpp"
+#include "../general/forall.hpp"
 
 using namespace std;
 
@@ -19,16 +20,21 @@ namespace mfem
 ComplexGridFunction::ComplexGridFunction(FiniteElementSpace *fes)
    : Vector(2*(fes->GetVSize()))
 {
-   gfr = new GridFunction(fes, data);
-   gfi = new GridFunction(fes, &data[fes->GetVSize()]);
+   UseDevice(true);
+   this->Vector::operator=(0.0);
+
+   gfr = new GridFunction();
+   gfr->MakeRef(fes, *this, 0);
+
+   gfi = new GridFunction();
+   gfi->MakeRef(fes, *this, fes->GetVSize());
 }
 
 void
 ComplexGridFunction::Update()
 {
-   FiniteElementSpace * fes = gfr->FESpace();
-
-   int vsize = fes->GetVSize();
+   FiniteElementSpace *fes = gfr->FESpace();
+   const int vsize = fes->GetVSize();
 
    const Operator *T = fes->GetUpdateOperator();
    if (T)
@@ -40,30 +46,36 @@ ComplexGridFunction::Update()
 
       // Our data array now contains old data as well as being the wrong size so
       // reallocate it.
+      UseDevice(true);
       this->SetSize(2 * vsize);
+      this->Vector::operator=(0.0);
 
       // Create temporary vectors which point to the new data array
-      Vector gf_r(data, vsize);
-      Vector gf_i((data) ? &data[vsize] : data, vsize);
+      Vector gf_r; gf_r.MakeRef(*this, 0, vsize);
+      Vector gf_i; gf_i.MakeRef(*this, vsize, vsize);
 
       // Copy the updated GridFunctions into the new data array
       gf_r = *gfr;
       gf_i = *gfi;
+      gf_r.SyncAliasMemory(*this);
+      gf_i.SyncAliasMemory(*this);
 
       // Replace the individual data arrays with pointers into the new data
       // array
-      gfr->NewDataAndSize(data, vsize);
-      gfi->NewDataAndSize((data) ? &data[vsize] : data, vsize);
+      gfr->MakeRef(*this, 0, vsize);
+      gfi->MakeRef(*this, vsize, vsize);
    }
    else
    {
       // The existing data will not be transferred to the new GridFunctions so
-      // delete it a allocate a new array
+      // delete it and allocate a new array
+      UseDevice(true);
       this->SetSize(2 * vsize);
+      this->Vector::operator=(0.0);
 
       // Point the individual GridFunctions to the new data array
-      gfr->NewDataAndSize(data, vsize);
-      gfi->NewDataAndSize((data) ? &data[vsize] : data, vsize);
+      gfr->MakeRef(*this, 0, vsize);
+      gfi->MakeRef(*this, vsize, vsize);
 
       // These updates will only set the proper 'sequence' value within the
       // individual GridFunction objects because their sizes are already correct
@@ -76,16 +88,24 @@ void
 ComplexGridFunction::ProjectCoefficient(Coefficient &real_coeff,
                                         Coefficient &imag_coeff)
 {
+   gfr->SyncMemory(*this);
+   gfi->SyncMemory(*this);
    gfr->ProjectCoefficient(real_coeff);
    gfi->ProjectCoefficient(imag_coeff);
+   gfr->SyncAliasMemory(*this);
+   gfi->SyncAliasMemory(*this);
 }
 
 void
 ComplexGridFunction::ProjectCoefficient(VectorCoefficient &real_vcoeff,
                                         VectorCoefficient &imag_vcoeff)
 {
+   gfr->SyncMemory(*this);
+   gfi->SyncMemory(*this);
    gfr->ProjectCoefficient(real_vcoeff);
    gfi->ProjectCoefficient(imag_vcoeff);
+   gfr->SyncAliasMemory(*this);
+   gfi->SyncAliasMemory(*this);
 }
 
 void
@@ -93,8 +113,12 @@ ComplexGridFunction::ProjectBdrCoefficient(Coefficient &real_coeff,
                                            Coefficient &imag_coeff,
                                            Array<int> &attr)
 {
+   gfr->SyncMemory(*this);
+   gfi->SyncMemory(*this);
    gfr->ProjectBdrCoefficient(real_coeff, attr);
    gfi->ProjectBdrCoefficient(imag_coeff, attr);
+   gfr->SyncAliasMemory(*this);
+   gfi->SyncAliasMemory(*this);
 }
 
 void
@@ -102,8 +126,12 @@ ComplexGridFunction::ProjectBdrCoefficientNormal(VectorCoefficient &real_vcoeff,
                                                  VectorCoefficient &imag_vcoeff,
                                                  Array<int> &attr)
 {
+   gfr->SyncMemory(*this);
+   gfi->SyncMemory(*this);
    gfr->ProjectBdrCoefficientNormal(real_vcoeff, attr);
    gfi->ProjectBdrCoefficientNormal(imag_vcoeff, attr);
+   gfr->SyncAliasMemory(*this);
+   gfi->SyncAliasMemory(*this);
 }
 
 void
@@ -113,18 +141,28 @@ ComplexGridFunction::ProjectBdrCoefficientTangent(VectorCoefficient
                                                   &imag_vcoeff,
                                                   Array<int> &attr)
 {
+   gfr->SyncMemory(*this);
+   gfi->SyncMemory(*this);
    gfr->ProjectBdrCoefficientTangent(real_vcoeff, attr);
    gfi->ProjectBdrCoefficientTangent(imag_vcoeff, attr);
+   gfr->SyncAliasMemory(*this);
+   gfi->SyncAliasMemory(*this);
 }
 
 
-ComplexLinearForm::ComplexLinearForm(FiniteElementSpace *f,
+ComplexLinearForm::ComplexLinearForm(FiniteElementSpace *fes,
                                      ComplexOperator::Convention convention)
-   : Vector(2*(f->GetVSize())),
+   : Vector(2*(fes->GetVSize())),
      conv(convention)
 {
-   lfr = new LinearForm(f, data);
-   lfi = new LinearForm(f, &data[f->GetVSize()]);
+   UseDevice(true);
+   this->Vector::operator=(0.0);
+
+   lfr = new LinearForm();
+   lfr->MakeRef(fes, *this, 0);
+
+   lfi = new LinearForm();
+   lfi->MakeRef(fes, *this, fes->GetVSize());
 }
 
 ComplexLinearForm::ComplexLinearForm(FiniteElementSpace *fes,
@@ -133,8 +171,14 @@ ComplexLinearForm::ComplexLinearForm(FiniteElementSpace *fes,
    : Vector(2*(fes->GetVSize())),
      conv(convention)
 {
-   lfr = new LinearForm(fes, lf_r);  lfr->SetData(data);
-   lfi = new LinearForm(fes, lf_i);  lfi->SetData(&data[fes->GetVSize()]);
+   UseDevice(true);
+   this->Vector::operator=(0.0);
+
+   lfr = new LinearForm(fes, lf_r);
+   lfi = new LinearForm(fes, lf_i);
+
+   lfr->MakeRef(fes, *this, 0);
+   lfi->MakeRef(fes, *this, fes->GetVSize());
 }
 
 ComplexLinearForm::~ComplexLinearForm()
@@ -189,42 +233,43 @@ void
 ComplexLinearForm::Update()
 {
    FiniteElementSpace *fes = lfr->FESpace();
-
    this->Update(fes);
 }
 
 void
 ComplexLinearForm::Update(FiniteElementSpace *fes)
 {
-   int vsize = fes->GetVSize();
-   SetSize(2 * vsize);
+   UseDevice(true);
+   SetSize(2 * fes->GetVSize());
+   this->Vector::operator=(0.0);
 
-   Vector vlfr(data, vsize);
-   Vector vlfi((data) ? &data[vsize] : data, vsize);
-
-   lfr->Update(fes, vlfr, 0);
-   lfi->Update(fes, vlfi, 0);
+   lfr->MakeRef(fes, *this, 0);
+   lfi->MakeRef(fes, *this, fes->GetVSize());
 }
 
 void
 ComplexLinearForm::Assemble()
 {
+   lfr->SyncMemory(*this);
+   lfi->SyncMemory(*this);
    lfr->Assemble();
    lfi->Assemble();
-   if (conv == ComplexOperator::BLOCK_SYMMETRIC)
-   {
-      *lfi *= -1.0;
-   }
+   if (conv == ComplexOperator::BLOCK_SYMMETRIC) { *lfi *= -1.0; }
+   lfr->SyncAliasMemory(*this);
+   lfi->SyncAliasMemory(*this);
 }
 
 complex<double>
 ComplexLinearForm::operator()(const ComplexGridFunction &gf) const
 {
-   double s = (conv == ComplexOperator::HERMITIAN)?1.0:-1.0;
+   double s = (conv == ComplexOperator::HERMITIAN) ? 1.0 : -1.0;
+   lfr->SyncMemory(*this);
+   lfi->SyncMemory(*this);
    return complex<double>((*lfr)(gf.real()) - s * (*lfi)(gf.imag()),
                           (*lfr)(gf.imag()) + s * (*lfi)(gf.real()));
 }
 
+
 bool SesquilinearForm::RealInteg()
 {
    int nint = blfr->GetFBFI()->Size() + blfr->GetDBFI()->Size() +
@@ -341,34 +386,45 @@ SesquilinearForm::FormLinearSystem(const Array<int> &ess_tdof_list,
                                    Vector &X, Vector &B,
                                    int ci)
 {
-   FiniteElementSpace * fes = blfr->FESpace();
-   int vsize  = fes->GetVSize();
+   FiniteElementSpace *fes = blfr->FESpace();
+   const int vsize = fes->GetVSize();
 
-   // Allocate temporary vectors
-   Vector b_0(vsize); b_0 = 0.0;
+   // Allocate temporary vector
+   Vector b_0;
+   b_0.UseDevice(true);
+   b_0.SetSize(vsize);
+   b_0 = 0.0;
 
    // Extract the real and imaginary parts of the input vectors
    MFEM_ASSERT(x.Size() == 2 * vsize, "Input GridFunction of incorrect size!");
-   Vector x_r(x.GetData(), vsize);
-   Vector x_i(&(x.GetData())[vsize], vsize);
+   x.Read();
+   Vector x_r; x_r.MakeRef(x, 0, vsize);
+   Vector x_i; x_i.MakeRef(x, vsize, vsize);
 
    MFEM_ASSERT(b.Size() == 2 * vsize, "Input LinearForm of incorrect size!");
-   Vector b_r(b.GetData(), vsize);
-   Vector b_i(&(b.GetData())[vsize], vsize);
+   b.Read();
+   Vector b_r; b_r.MakeRef(b, 0, vsize);
+   Vector b_i; b_i.MakeRef(b, vsize, vsize);
 
    if (conv == ComplexOperator::BLOCK_SYMMETRIC) { b_i *= -1.0; }
 
-   int tvsize = fes->GetTrueVSize();
+   const int tvsize = fes->GetTrueVSize();
    OperatorHandle A_r, A_i;
 
+   X.UseDevice(true);
    X.SetSize(2 * tvsize);
-   B.SetSize(2 * tvsize);
+   X = 0.0;
 
-   Vector X_0(tvsize), B_0(tvsize);
-   Vector X_r(X.GetData(),tvsize);
-   Vector X_i(&(X.GetData())[tvsize], tvsize);
-   Vector B_r(B.GetData(), tvsize);
-   Vector B_i(&(B.GetData())[tvsize], tvsize);
+   B.UseDevice(true);
+   B.SetSize(2 * tvsize);
+   B = 0.0;
+
+   Vector X_r; X_r.MakeRef(X, 0, tvsize);
+   Vector X_i; X_i.MakeRef(X, tvsize, tvsize);
+   Vector B_r; B_r.MakeRef(B, 0, tvsize);
+   Vector B_i; B_i.MakeRef(B, tvsize, tvsize);
+
+   Vector X_0, B_0;
 
    if (RealInteg())
    {
@@ -418,13 +474,18 @@ SesquilinearForm::FormLinearSystem(const Array<int> &ess_tdof_list,
       // conform with standard essential BC treatment
       if (A_i.Is<ConstrainedOperator>())
       {
-         int n = ess_tdof_list.Size();
-         for (int k = 0; k < n; k++)
+         const int n = ess_tdof_list.Size();
+         auto d_B_r = B_r.Write();
+         auto d_B_i = B_i.Write();
+         auto d_X_r = X_r.Read();
+         auto d_X_i = X_i.Read();
+         auto d_idx = ess_tdof_list.Read();
+         MFEM_FORALL(i, n,
          {
-            int j = ess_tdof_list[k];
-            B_r(j) = X_r(j);
-            B_i(j) = X_i(j);
-         }
+            const int j = d_idx[i];
+            d_B_r[j] = d_X_r[j];
+            d_B_i[j] = d_X_i[j];
+         });
          A_i.As<ConstrainedOperator>()->SetDiagonalPolicy
          (mfem::Operator::DiagonalPolicy::DIAG_ZERO);
       }
@@ -436,6 +497,16 @@ SesquilinearForm::FormLinearSystem(const Array<int> &ess_tdof_list,
       b_i *= -1.0;
    }
 
+   x_r.SyncAliasMemory(x);
+   x_i.SyncAliasMemory(x);
+   b_r.SyncAliasMemory(b);
+   b_i.SyncAliasMemory(b);
+
+   X_r.SyncAliasMemory(X);
+   X_i.SyncAliasMemory(X);
+   B_r.SyncAliasMemory(B);
+   B_i.SyncAliasMemory(B);
+
    // A = A_r + i A_i
    A.Clear();
    if ( A_r.Type() == Operator::MFEM_SPARSEMAT ||
@@ -528,29 +599,32 @@ void
 SesquilinearForm::RecoverFEMSolution(const Vector &X, const Vector &b,
                                      Vector &x)
 {
-   FiniteElementSpace * fes = blfr->FESpace();
+   FiniteElementSpace *fes = blfr->FESpace();
 
    const SparseMatrix *P = fes->GetConformingProlongation();
-
-   int vsize  = fes->GetVSize();
-   int tvsize = X.Size() / 2;
-
-   Vector X_r(X.GetData(), tvsize);
-   Vector X_i(&(X.GetData())[tvsize], tvsize);
-
-   Vector x_r(x.GetData(), vsize);
-   Vector x_i(&(x.GetData())[vsize], vsize);
-
    if (!P)
    {
       x = X;
+      return;
    }
-   else
-   {
-      // Apply conforming prolongation
-      P->Mult(X_r, x_r);
-      P->Mult(X_i, x_i);
-   }
+
+   const int vsize  = fes->GetVSize();
+   const int tvsize = X.Size() / 2;
+
+   X.Read();
+   Vector X_r; X_r.MakeRef(const_cast<Vector&>(X), 0, tvsize);
+   Vector X_i; X_i.MakeRef(const_cast<Vector&>(X), tvsize, tvsize);
+
+   x.Write();
+   Vector x_r; x_r.MakeRef(x, 0, vsize);
+   Vector x_i; x_i.MakeRef(x, vsize, vsize);
+
+   // Apply conforming prolongation
+   P->Mult(X_r, x_r);
+   P->Mult(X_i, x_i);
+
+   x_r.SyncAliasMemory(x);
+   x_i.SyncAliasMemory(x);
 }
 
 void
@@ -566,16 +640,21 @@ SesquilinearForm::Update(FiniteElementSpace *nfes)
 ParComplexGridFunction::ParComplexGridFunction(ParFiniteElementSpace *pfes)
    : Vector(2*(pfes->GetVSize()))
 {
-   pgfr = new ParGridFunction(pfes, data);
-   pgfi = new ParGridFunction(pfes, (data) ? &data[pfes->GetVSize()]:data);
+   UseDevice(true);
+   this->Vector::operator=(0.0);
+
+   pgfr = new ParGridFunction();
+   pgfr->MakeRef(pfes, *this, 0);
+
+   pgfi = new ParGridFunction();
+   pgfi->MakeRef(pfes, *this, pfes->GetVSize());
 }
 
 void
 ParComplexGridFunction::Update()
 {
-   ParFiniteElementSpace * pfes = pgfr->ParFESpace();
-
-   int vsize = pfes->GetVSize();
+   ParFiniteElementSpace *pfes = pgfr->ParFESpace();
+   const int vsize = pfes->GetVSize();
 
    const Operator *T = pfes->GetUpdateOperator();
    if (T)
@@ -587,30 +666,34 @@ ParComplexGridFunction::Update()
 
       // Our data array now contains old data as well as being the wrong size so
       // reallocate it.
+      UseDevice(true);
       this->SetSize(2 * vsize);
+      this->Vector::operator=(0.0);
 
       // Create temporary vectors which point to the new data array
-      Vector gf_r(data, vsize);
-      Vector gf_i((data) ? &data[vsize] : data, vsize);
+      Vector gf_r; gf_r.MakeRef(*this, 0, vsize);
+      Vector gf_i; gf_i.MakeRef(*this, vsize, vsize);
 
       // Copy the updated GridFunctions into the new data array
-      gf_r = *pgfr;
-      gf_i = *pgfi;
+      gf_r = *pgfr; gf_r.SyncAliasMemory(*this);
+      gf_i = *pgfi; gf_i.SyncAliasMemory(*this);
 
       // Replace the individual data arrays with pointers into the new data
       // array
-      pgfr->NewDataAndSize(data, vsize);
-      pgfi->NewDataAndSize((data) ? &data[vsize] : data, vsize);
+      pgfr->MakeRef(*this, 0, vsize);
+      pgfi->MakeRef(*this, vsize, vsize);
    }
    else
    {
       // The existing data will not be transferred to the new GridFunctions so
-      // delete it a allocate a new array
+      // delete it and allocate a new array
+      UseDevice(true);
       this->SetSize(2 * vsize);
+      this->Vector::operator=(0.0);
 
       // Point the individual GridFunctions to the new data array
-      pgfr->NewDataAndSize(data, vsize);
-      pgfi->NewDataAndSize((data) ? &data[vsize] : data, vsize);
+      pgfr->MakeRef(*this, 0, vsize);
+      pgfi->MakeRef(*this, vsize, vsize);
 
       // These updates will only set the proper 'sequence' value within the
       // individual GridFunction objects because their sizes are already correct
@@ -623,16 +706,24 @@ void
 ParComplexGridFunction::ProjectCoefficient(Coefficient &real_coeff,
                                            Coefficient &imag_coeff)
 {
+   pgfr->SyncMemory(*this);
+   pgfi->SyncMemory(*this);
    pgfr->ProjectCoefficient(real_coeff);
    pgfi->ProjectCoefficient(imag_coeff);
+   pgfr->SyncAliasMemory(*this);
+   pgfi->SyncAliasMemory(*this);
 }
 
 void
 ParComplexGridFunction::ProjectCoefficient(VectorCoefficient &real_vcoeff,
                                            VectorCoefficient &imag_vcoeff)
 {
+   pgfr->SyncMemory(*this);
+   pgfi->SyncMemory(*this);
    pgfr->ProjectCoefficient(real_vcoeff);
    pgfi->ProjectCoefficient(imag_vcoeff);
+   pgfr->SyncAliasMemory(*this);
+   pgfi->SyncAliasMemory(*this);
 }
 
 void
@@ -640,8 +731,12 @@ ParComplexGridFunction::ProjectBdrCoefficient(Coefficient &real_coeff,
                                               Coefficient &imag_coeff,
                                               Array<int> &attr)
 {
+   pgfr->SyncMemory(*this);
+   pgfi->SyncMemory(*this);
    pgfr->ProjectBdrCoefficient(real_coeff, attr);
    pgfi->ProjectBdrCoefficient(imag_coeff, attr);
+   pgfr->SyncAliasMemory(*this);
+   pgfi->SyncAliasMemory(*this);
 }
 
 void
@@ -651,8 +746,12 @@ ParComplexGridFunction::ProjectBdrCoefficientNormal(VectorCoefficient
                                                     &imag_vcoeff,
                                                     Array<int> &attr)
 {
+   pgfr->SyncMemory(*this);
+   pgfi->SyncMemory(*this);
    pgfr->ProjectBdrCoefficientNormal(real_vcoeff, attr);
    pgfi->ProjectBdrCoefficientNormal(imag_vcoeff, attr);
+   pgfr->SyncAliasMemory(*this);
+   pgfi->SyncAliasMemory(*this);
 }
 
 void
@@ -662,36 +761,51 @@ ParComplexGridFunction::ProjectBdrCoefficientTangent(VectorCoefficient
                                                      &imag_vcoeff,
                                                      Array<int> &attr)
 {
+   pgfr->SyncMemory(*this);
+   pgfi->SyncMemory(*this);
    pgfr->ProjectBdrCoefficientTangent(real_vcoeff, attr);
    pgfi->ProjectBdrCoefficientTangent(imag_vcoeff, attr);
+   pgfr->SyncAliasMemory(*this);
+   pgfi->SyncAliasMemory(*this);
 }
 
 void
 ParComplexGridFunction::Distribute(const Vector *tv)
 {
-   ParFiniteElementSpace * pfes = pgfr->ParFESpace();
-   HYPRE_Int size = pfes->GetTrueVSize();
+   ParFiniteElementSpace *pfes = pgfr->ParFESpace();
+   const int tvsize = pfes->GetTrueVSize();
 
-   double * tvd = tv->GetData();
-   Vector tvr(tvd, size);
-   Vector tvi((tvd) ? &tvd[size] : tvd, size);
+   tv->Read();
+   Vector tvr; tvr.MakeRef(const_cast<Vector&>(*tv), 0, tvsize);
+   Vector tvi; tvi.MakeRef(const_cast<Vector&>(*tv), tvsize, tvsize);
 
+   pgfr->SyncMemory(*this);
+   pgfi->SyncMemory(*this);
    pgfr->Distribute(tvr);
    pgfi->Distribute(tvi);
+   pgfr->SyncAliasMemory(*this);
+   pgfi->SyncAliasMemory(*this);
 }
 
 void
 ParComplexGridFunction::ParallelProject(Vector &tv) const
 {
-   ParFiniteElementSpace * pfes = pgfr->ParFESpace();
-   HYPRE_Int size = pfes->GetTrueVSize();
+   ParFiniteElementSpace *pfes = pgfr->ParFESpace();
+   const int tvsize = pfes->GetTrueVSize();
 
-   double * tvd = tv.GetData();
-   Vector tvr(tvd, size);
-   Vector tvi((tvd) ? &tvd[size] : tvd, size);
+   tv.Write();
+   Vector tvr; tvr.MakeRef(tv, 0, tvsize);
+   Vector tvi; tvi.MakeRef(tv, tvsize, tvsize);
 
+   pgfr->SyncMemory(*this);
+   pgfi->SyncMemory(*this);
    pgfr->ParallelProject(tvr);
    pgfi->ParallelProject(tvi);
+   pgfr->SyncAliasMemory(*this);
+   pgfi->SyncAliasMemory(*this);
+
+   tvr.SyncAliasMemory(tv);
+   tvi.SyncAliasMemory(tv);
 }
 
 
@@ -701,10 +815,16 @@ ParComplexLinearForm::ParComplexLinearForm(ParFiniteElementSpace *pfes,
    : Vector(2*(pfes->GetVSize())),
      conv(convention)
 {
-   plfr = new ParLinearForm(pfes, data);
-   plfi = new ParLinearForm(pfes, (data) ? &data[pfes->GetVSize()]:data);
+   UseDevice(true);
+   this->Vector::operator=(0.0);
 
-   HYPRE_Int * tdof_offsets_fes = pfes->GetTrueDofOffsets();
+   plfr = new ParLinearForm();
+   plfr->MakeRef(pfes, *this, 0);
+
+   plfi = new ParLinearForm();
+   plfi->MakeRef(pfes, *this, pfes->GetVSize());
+
+   HYPRE_Int *tdof_offsets_fes = pfes->GetTrueDofOffsets();
 
    int n = (HYPRE_AssumedPartitionCheck()) ? 2 : pfes->GetNRanks();
    tdof_offsets = new HYPRE_Int[n+1];
@@ -724,12 +844,16 @@ ParComplexLinearForm::ParComplexLinearForm(ParFiniteElementSpace *pfes,
    : Vector(2*(pfes->GetVSize())),
      conv(convention)
 {
-   plfr = new ParLinearForm(pfes, plf_r);
-   plfr->SetData(data);
-   plfi = new ParLinearForm(pfes, plf_i);
-   plfi->SetData((data) ? &data[pfes->GetVSize()]:data);
+   UseDevice(true);
+   this->Vector::operator=(0.0);
 
-   HYPRE_Int * tdof_offsets_fes = pfes->GetTrueDofOffsets();
+   plfr = new ParLinearForm(pfes, plf_r);
+   plfi = new ParLinearForm(pfes, plf_i);
+
+   plfr->MakeRef(pfes, *this, 0);
+   plfi->MakeRef(pfes, *this, pfes->GetVSize());
+
+   HYPRE_Int *tdof_offsets_fes = pfes->GetTrueDofOffsets();
 
    int n = (HYPRE_AssumedPartitionCheck()) ? 2 : pfes->GetNRanks();
    tdof_offsets = new HYPRE_Int[n+1];
@@ -792,58 +916,71 @@ ParComplexLinearForm::AddBdrFaceIntegrator(LinearFormIntegrator *lfi_real,
 void
 ParComplexLinearForm::Update(ParFiniteElementSpace *pf)
 {
-   ParFiniteElementSpace *pfes = (pf!=NULL)?pf:plfr->ParFESpace();
-   int vsize = pfes->GetVSize();
-   SetSize(2 * vsize);
+   ParFiniteElementSpace *pfes = (pf != NULL) ? pf : plfr->ParFESpace();
 
-   Vector vplfr(data, vsize);
-   Vector vplfi((data) ? &data[vsize] : data, vsize);
+   UseDevice(true);
+   SetSize(2 * pfes->GetVSize());
+   this->Vector::operator=(0.0);
 
-   plfr->Update(pfes, vplfr, 0);
-   plfi->Update(pfes, vplfi, 0);
+   plfr->MakeRef(pfes, *this, 0);
+   plfi->MakeRef(pfes, *this, pfes->GetVSize());
 }
 
 void
 ParComplexLinearForm::Assemble()
 {
+   plfr->SyncMemory(*this);
+   plfi->SyncMemory(*this);
    plfr->Assemble();
    plfi->Assemble();
-   if (conv == ComplexOperator::BLOCK_SYMMETRIC)
-   {
-      *plfi *= -1.0;
-   }
+   if (conv == ComplexOperator::BLOCK_SYMMETRIC) { *plfi *= -1.0; }
+   plfr->SyncAliasMemory(*this);
+   plfi->SyncAliasMemory(*this);
 }
 
 void
 ParComplexLinearForm::ParallelAssemble(Vector &tv)
 {
-   HYPRE_Int size = plfr->ParFESpace()->GetTrueVSize();
+   const int tvsize = plfr->ParFESpace()->GetTrueVSize();
 
-   double * tvd = tv.GetData();
-   Vector tvr(tvd, size);
-   Vector tvi((tvd) ? &tvd[size] : tvd, size);
+   tv.Write();
+   Vector tvr; tvr.MakeRef(tv, 0, tvsize);
+   Vector tvi; tvi.MakeRef(tv, tvsize, tvsize);
 
+   plfr->SyncMemory(*this);
+   plfi->SyncMemory(*this);
    plfr->ParallelAssemble(tvr);
    plfi->ParallelAssemble(tvi);
+   plfr->SyncAliasMemory(*this);
+   plfi->SyncAliasMemory(*this);
+
+   tvr.SyncAliasMemory(tv);
+   tvi.SyncAliasMemory(tv);
 }
 
 HypreParVector *
 ParComplexLinearForm::ParallelAssemble()
 {
-   const ParFiniteElementSpace * pfes = plfr->ParFESpace();
+   const ParFiniteElementSpace *pfes = plfr->ParFESpace();
+   const int tvsize = pfes->GetTrueVSize();
 
-   HypreParVector * tv = new HypreParVector(pfes->GetComm(),
-                                            2*(pfes->GlobalTrueVSize()),
-                                            tdof_offsets);
+   HypreParVector *tv = new HypreParVector(pfes->GetComm(),
+                                           2*(pfes->GlobalTrueVSize()),
+                                           tdof_offsets);
 
-   HYPRE_Int size = pfes->GetTrueVSize();
-
-   double * tvd = tv->GetData();
-   Vector tvr(tvd, size);
-   Vector tvi((tvd) ? &tvd[size] : tvd, size);
+   tv->Write();
+   Vector tvr; tvr.MakeRef(*tv, 0, tvsize);
+   Vector tvi; tvi.MakeRef(*tv, tvsize, tvsize);
 
+   plfr->SyncMemory(*this);
+   plfi->SyncMemory(*this);
    plfr->ParallelAssemble(tvr);
    plfi->ParallelAssemble(tvi);
+   plfr->SyncAliasMemory(*this);
+   plfi->SyncAliasMemory(*this);
+
+   tvr.SyncAliasMemory(*tv);
+   tvi.SyncAliasMemory(*tv);
 
    return tv;
 }
@@ -851,13 +988,14 @@ ParComplexLinearForm::ParallelAssemble()
 complex<double>
 ParComplexLinearForm::operator()(const ParComplexGridFunction &gf) const
 {
-   double s = (conv == ComplexOperator::HERMITIAN)?1.0:-1.0;
+   plfr->SyncMemory(*this);
+   plfi->SyncMemory(*this);
+   double s = (conv == ComplexOperator::HERMITIAN) ? 1.0 : -1.0;
    return complex<double>((*plfr)(gf.real()) - s * (*plfi)(gf.imag()),
                           (*plfr)(gf.imag()) + s * (*plfi)(gf.real()));
 }
 
 
-
 bool ParSesquilinearForm::RealInteg()
 {
    int nint = pblfr->GetFBFI()->Size() + pblfr->GetDBFI()->Size() +
@@ -964,7 +1102,6 @@ ParSesquilinearForm::ParallelAssemble()
    return new ComplexHypreParMatrix(pblfr->ParallelAssemble(),
                                     pblfi->ParallelAssemble(),
                                     true, true, conv);
-
 }
 
 void
@@ -974,35 +1111,45 @@ ParSesquilinearForm::FormLinearSystem(const Array<int> &ess_tdof_list,
                                       Vector &X, Vector &B,
                                       int ci)
 {
-   ParFiniteElementSpace * pfes = pblfr->ParFESpace();
-   int vsize = pfes->GetVSize();
+   ParFiniteElementSpace *pfes = pblfr->ParFESpace();
+   const int vsize = pfes->GetVSize();
 
-   // Allocate temporary vectors
-   Vector b_0(vsize); b_0 = 0.0;
+   // Allocate temporary vector
+   Vector b_0;
+   b_0.UseDevice(true);
+   b_0.SetSize(vsize);
+   b_0 = 0.0;
 
    // Extract the real and imaginary parts of the input vectors
    MFEM_ASSERT(x.Size() == 2 * vsize, "Input GridFunction of incorrect size!");
-   Vector x_r(x.GetData(), vsize);
-   Vector x_i(&(x.GetData())[vsize], vsize);
+   x.Read();
+   Vector x_r; x_r.MakeRef(x, 0, vsize);
+   Vector x_i; x_i.MakeRef(x, vsize, vsize);
 
    MFEM_ASSERT(b.Size() == 2 * vsize, "Input LinearForm of incorrect size!");
-   Vector b_r(b.GetData(), vsize);
-   Vector b_i(&(b.GetData())[vsize], vsize);
+   b.Read();
+   Vector b_r; b_r.MakeRef(b, 0, vsize);
+   Vector b_i; b_i.MakeRef(b, vsize, vsize);
 
    if (conv == ComplexOperator::BLOCK_SYMMETRIC) { b_i *= -1.0; }
 
-   int tvsize = pfes->GetTrueVSize();
-
+   const int tvsize = pfes->GetTrueVSize();
    OperatorHandle A_r, A_i;
 
+   X.UseDevice(true);
    X.SetSize(2 * tvsize);
-   B.SetSize(2 * tvsize);
+   X = 0.0;
 
-   Vector X_0(tvsize), B_0(tvsize);
-   Vector X_r(X.GetData(),tvsize);
-   Vector X_i(&(X.GetData())[tvsize], tvsize);
-   Vector B_r(B.GetData(), tvsize);
-   Vector B_i(&(B.GetData())[tvsize], tvsize);
+   B.UseDevice(true);
+   B.SetSize(2 * tvsize);
+   B = 0.0;
+
+   Vector X_r; X_r.MakeRef(X, 0, tvsize);
+   Vector X_i; X_i.MakeRef(X, tvsize, tvsize);
+   Vector B_r; B_r.MakeRef(B, 0, tvsize);
+   Vector B_i; B_i.MakeRef(B, tvsize, tvsize);
+
+   Vector X_0, B_0;
 
    if (RealInteg())
    {
@@ -1042,24 +1189,29 @@ ParSesquilinearForm::FormLinearSystem(const Array<int> &ess_tdof_list,
 
    if (RealInteg() && ImagInteg())
    {
-      int n = ess_tdof_list.Size();
       // Modify RHS to conform with standard essential BC treatment
-      for (int k = 0; k < n; k++)
+      const int n = ess_tdof_list.Size();
+      auto d_B_r = B_r.Write();
+      auto d_B_i = B_i.Write();
+      auto d_X_r = X_r.Read();
+      auto d_X_i = X_i.Read();
+      auto d_idx = ess_tdof_list.Read();
+      MFEM_FORALL(i, n,
       {
-         int j=ess_tdof_list[k];
-         B_r(j) = X_r(j);
-         B_i(j) = X_i(j);
-      }
+         const int j = d_idx[i];
+         d_B_r[j] = d_X_r[j];
+         d_B_i[j] = d_X_i[j];
+      });
       // Modify offdiagonal blocks (imaginary parts of the matrix) to conform
       // with standard essential BC treatment
-      if ( A_i.Type() == Operator::Hypre_ParCSR )
+      if (A_i.Type() == Operator::Hypre_ParCSR)
       {
          HypreParMatrix * Ah;
          A_i.Get(Ah);
          hypre_ParCSRMatrix *Aih = *Ah;
          for (int k = 0; k < n; k++)
          {
-            int j = ess_tdof_list[k];
+            const int j = ess_tdof_list[k];
             Aih->diag->data[Aih->diag->i[j]] = 0.0;
          }
       }
@@ -1076,6 +1228,16 @@ ParSesquilinearForm::FormLinearSystem(const Array<int> &ess_tdof_list,
       b_i *= -1.0;
    }
 
+   x_r.SyncAliasMemory(x);
+   x_i.SyncAliasMemory(x);
+   b_r.SyncAliasMemory(b);
+   b_i.SyncAliasMemory(b);
+
+   X_r.SyncAliasMemory(X);
+   X_i.SyncAliasMemory(X);
+   B_r.SyncAliasMemory(B);
+   B_i.SyncAliasMemory(B);
+
    // A = A_r + i A_i
    A.Clear();
    if ( A_r.Type() == Operator::Hypre_ParCSR ||
@@ -1175,22 +1337,27 @@ void
 ParSesquilinearForm::RecoverFEMSolution(const Vector &X, const Vector &b,
                                         Vector &x)
 {
-   ParFiniteElementSpace * pfes = pblfr->ParFESpace();
+   ParFiniteElementSpace *pfes = pblfr->ParFESpace();
 
    const Operator &P = *pfes->GetProlongationMatrix();
 
-   int vsize  = pfes->GetVSize();
-   int tvsize = X.Size() / 2;
+   const int vsize  = pfes->GetVSize();
+   const int tvsize = X.Size() / 2;
 
-   Vector X_r(X.GetData(), tvsize);
-   Vector X_i(&(X.GetData())[tvsize], tvsize);
+   X.Read();
+   Vector X_r; X_r.MakeRef(const_cast<Vector&>(X), 0, tvsize);
+   Vector X_i; X_i.MakeRef(const_cast<Vector&>(X), tvsize, tvsize);
 
-   Vector x_r(x.GetData(), vsize);
-   Vector x_i(&(x.GetData())[vsize], vsize);
+   x.Write();
+   Vector x_r; x_r.MakeRef(x, 0, vsize);
+   Vector x_i; x_i.MakeRef(x, vsize, vsize);
 
    // Apply conforming prolongation
    P.Mult(X_r, x_r);
    P.Mult(X_i, x_i);
+
+   x_r.SyncAliasMemory(x);
+   x_i.SyncAliasMemory(x);
 }
 
 void

fem/complex_fem.hpp

@@ -99,8 +99,8 @@ public:
                      ComplexOperator::Convention
                      convention = ComplexOperator::HERMITIAN);
 
-   /** @brief Create a ComplexLinearForm on the FiniteElementSpace @a f, using
-       the same integrators as the LinearForms @a lfr (real) and @a lfi (imag) .
+   /** @brief Create a ComplexLinearForm on the FiniteElementSpace @a fes, using
+       the same integrators as the LinearForms @a lf_r (real) and @a lf_i (imag).
 
        The pointer @a fes is not owned by the newly constructed object.
 
@@ -195,8 +195,8 @@ private:
    BilinearForm *blfr;
    BilinearForm *blfi;
 
-   /* These methods check if the real/imag parts of the sesqulinear form are not
-      empty */
+   /* These methods check if the real/imag parts of the sesquilinear form are
+      not empty */
    bool RealInteg();
    bool ImagInteg();
 
@@ -204,7 +204,7 @@ public:
    SesquilinearForm(FiniteElementSpace *fes,
                     ComplexOperator::Convention
                     convention = ComplexOperator::HERMITIAN);
-   /** @brief Create a SesquilinearForm on the FiniteElementSpace @a f, using
+   /** @brief Create a SesquilinearForm on the FiniteElementSpace @a fes, using
        the same integrators as the BilinearForms @a bfr and @a bfi .
 
        The pointer @a fes is not owned by the newly constructed object.
@@ -324,7 +324,7 @@ protected:
 public:
 
    /* @brief Construct a ParComplexGridFunction associated with the
-      ParFiniteElementSpace @a *f. */
+      ParFiniteElementSpace @a *pf. */
    ParComplexGridFunction(ParFiniteElementSpace *pf);
 
    void Update();
@@ -416,8 +416,8 @@ public:
                         convention = ComplexOperator::HERMITIAN);
 
    /** @brief Create a ParComplexLinearForm on the ParFiniteElementSpace @a pf,
-       using the same integrators as the LinearForms @a plfr (real) and @a plfi
-       (imag) .
+       using the same integrators as the LinearForms @a plf_r (real) and
+       @a plf_i (imag).
 
       The pointer @a fes is not owned by the newly constructed object.
 

fem/gridfunc.cpp

@@ -199,8 +199,7 @@ void GridFunction::MakeRef(FiniteElementSpace *f, Vector &v, int v_offset)
    if (f != fes) { Destroy(); }
    fes = f;
    v.UseDevice(true);
-   NewMemoryAndSize(Memory<double>(v.GetMemory(), v_offset, fes->GetVSize()),
-                    fes->GetVSize(), true);
+   this->Vector::MakeRef(v, v_offset, fes->GetVSize());
    sequence = fes->GetSequence();
 }
 

fem/linearform.cpp

@@ -204,6 +204,14 @@ void LinearForm::Update(FiniteElementSpace *f, Vector &v, int v_offset)
    ResetDeltaLocations();
 }
 
+void LinearForm::MakeRef(FiniteElementSpace *f, Vector &v, int v_offset)
+{
+   MFEM_ASSERT(v.Size() >= v_offset + f->GetVSize(), "");
+   fes = f;
+   v.UseDevice(true);
+   this->Vector::MakeRef(v, v_offset, fes->GetVSize());
+}
+
 void LinearForm::AssembleDelta()
 {
    if (dlfi_delta.Size() == 0) { return; }

fem/linearform.hpp

@@ -26,7 +26,7 @@ protected:
    /// FE space on which the LinearForm lives. Not owned.
    FiniteElementSpace *fes;
 
-   /** @brief Indicates the LinerFormIntegrator%s stored in #dlfi, #dlfi_delta,
+   /** @brief Indicates the LinearFormIntegrator%s stored in #dlfi, #dlfi_delta,
        #blfi, and #flfi are owned by another LinearForm. */
    int extern_lfs;
 
@@ -175,6 +175,15 @@ public:
        @note This method does not perform assembly. */
    void Update(FiniteElementSpace *f, Vector &v, int v_offset);
 
+   /** @brief Make the LinearForm reference external data on a new
+       FiniteElementSpace. */
+   /** This method changes the FiniteElementSpace associated with the LinearForm
+       @a *f and sets the data of the Vector @a v (plus the @a v_offset)
+       as external data in the LinearForm.
+       @note This version of the method will also perform bounds checks when
+       the build option MFEM_DEBUG is enabled. */
+   virtual void MakeRef(FiniteElementSpace *f, Vector &v, int v_offset);
+
    /// Return the action of the LinearForm as a linear mapping.
    /** Linear forms are linear functionals which map GridFunctions to
        the real numbers.  This method performs this mapping which in

fem/plinearform.cpp

@@ -21,7 +21,6 @@ namespace mfem
 void ParLinearForm::Update(ParFiniteElementSpace *pf)
 {
    if (pf) { pfes = pf; }
-
    LinearForm::Update(pfes);
 }
 
@@ -31,6 +30,18 @@ void ParLinearForm::Update(ParFiniteElementSpace *pf, Vector &v, int v_offset)
    LinearForm::Update(pf,v,v_offset);
 }
 
+void ParLinearForm::MakeRef(FiniteElementSpace *f, Vector &v, int v_offset)
+{
+   LinearForm::MakeRef(f, v, v_offset);
+   pfes = dynamic_cast<ParFiniteElementSpace*>(f);
+}
+
+void ParLinearForm::MakeRef(ParFiniteElementSpace *pf, Vector &v, int v_offset)
+{
+   LinearForm::MakeRef(pf, v, v_offset);
+   pfes = pf;
+}
+
 void ParLinearForm::ParallelAssemble(Vector &tv)
 {
    const Operator* prolong = pfes->GetProlongationMatrix();

fem/plinearform.hpp

@@ -92,6 +92,25 @@ public:
        @note This method does not perform assembly. */
    void Update(ParFiniteElementSpace *pf, Vector &v, int v_offset);
 
+
+   /** @brief Make the ParLinearForm reference external data on a new
+       FiniteElementSpace. */
+   /** This method changes the FiniteElementSpace associated with the ParLinearForm
+       to @a *f and sets the data of the Vector @a v (plus the @a v_offset) as external
+       data in the ParLinearForm.
+       @note This version of the method will also perform bounds checks when
+       the build option MFEM_DEBUG is enabled. */
+   virtual void MakeRef(FiniteElementSpace *f, Vector &v, int v_offset);
+
+   /** @brief Make the ParLinearForm reference external data on a new
+       ParFiniteElementSpace. */
+   /** This method changes the ParFiniteElementSpace associated with the ParLinearForm
+       to @a *pf and sets the data of the Vector @a v (plus the @a v_offset) as external
+       data in the ParLinearForm.
+       @note This version of the method will also perform bounds checks when
+       the build option MFEM_DEBUG is enabled. */
+   void MakeRef(ParFiniteElementSpace *pf, Vector &v, int v_offset);
+
    /// Assemble the vector on the true dofs, i.e. P^t v.
    void ParallelAssemble(Vector &tv);
 

linalg/complex_operator.cpp

@@ -26,10 +26,10 @@ ComplexOperator::ComplexOperator(Operator * Op_Real, Operator * Op_Imag,
    , ownReal_(ownReal)
    , ownImag_(ownImag)
    , convention_(convention)
-   , x_r_(NULL, width / 2)
-   , x_i_(NULL, width / 2)
-   , y_r_(NULL, height / 2)
-   , y_i_(NULL, height / 2)
+   , x_r_()
+   , x_i_()
+   , y_r_()
+   , y_i_()
    , u_(NULL)
    , v_(NULL)
 {}
@@ -68,14 +68,26 @@ const Operator & ComplexOperator::imag() const
 
 void ComplexOperator::Mult(const Vector &x, Vector &y) const
 {
-   double * x_data = x.GetData();
-   x_r_.SetData(x_data);
-   x_i_.SetData(&x_data[width / 2]);
+   x.Read();
+   y.UseDevice(true); y = 0.0;
 
-   y_r_.SetData(&y[0]);
-   y_i_.SetData(&y[height / 2]);
+   x_r_.MakeRef(const_cast<Vector&>(x), 0, width/2);
+   x_i_.MakeRef(const_cast<Vector&>(x), width/2, width/2);
+
+   y_r_.MakeRef(y, 0, height/2);
+   y_i_.MakeRef(y, height/2, height/2);
 
    this->Mult(x_r_, x_i_, y_r_, y_i_);
+
+   y_r_.SyncAliasMemory(y);
+   y_i_.SyncAliasMemory(y);
+
+   // Destroy alias vectors to prevent dangling aliases when the base vectors
+   // are deleted
+   x_r_.Destroy();
+   x_i_.Destroy();
+   y_r_.Destroy();
+   y_i_.Destroy();
 }
 
 void ComplexOperator::Mult(const Vector &x_r, const Vector &x_i,
@@ -91,31 +103,47 @@ void ComplexOperator::Mult(const Vector &x_r, const Vector &x_i,
       y_r = 0.0;
       y_i = 0.0;
    }
+
    if (Op_Imag_)
    {
-      if (!v_) { v_ = new Vector(Op_Imag_->Height()); }
+      if (!v_) { v_ = new Vector(); }
+      v_->UseDevice(true);
+      v_->SetSize(Op_Imag_->Height());
+
       Op_Imag_->Mult(x_i, *v_);
-      y_r_ -= *v_;
+      y_r.Add(-1.0, *v_);
       Op_Imag_->Mult(x_r, *v_);
-      y_i_ += *v_;
+      y_i.Add(1.0, *v_);
    }
 
    if (convention_ == BLOCK_SYMMETRIC)
    {
-      y_i_ *= -1.0;
+      y_i *= -1.0;
    }
 }
 
 void ComplexOperator::MultTranspose(const Vector &x, Vector &y) const
 {
-   double * x_data = x.GetData();
-   y_r_.SetData(x_data);
-   y_i_.SetData(&x_data[height / 2]);
+   x.Read();
+   y.UseDevice(true); y = 0.0;
 
-   x_r_.SetData(&y[0]);
-   x_i_.SetData(&y[width / 2]);
+   x_r_.MakeRef(const_cast<Vector&>(x), 0, height/2);
+   x_i_.MakeRef(const_cast<Vector&>(x), height/2, height/2);
 
-   this->MultTranspose(y_r_, y_i_, x_r_, x_i_);
+   y_r_.MakeRef(y, 0, width/2);
+   y_i_.MakeRef(y, width/2, width/2);
+
+   this->MultTranspose(x_r_, x_i_, y_r_, y_i_);
+
+   y_r_.SyncAliasMemory(y);
+   y_i_.SyncAliasMemory(y);
+
+   // Destroy alias vectors to prevent dangling aliases when the base vectors
+   // are deleted
+   x_r_.Destroy();
+   x_i_.Destroy();
+   y_r_.Destroy();
+   y_i_.Destroy();
 }
 
 void ComplexOperator::MultTranspose(const Vector &x_r, const Vector &x_i,
@@ -136,13 +164,17 @@ void ComplexOperator::MultTranspose(const Vector &x_r, const Vector &x_i,
       y_r = 0.0;
       y_i = 0.0;
    }
+
    if (Op_Imag_)
    {
-      if (!u_) { u_ = new Vector(Op_Imag_->Width()); }
+      if (!u_) { u_ = new Vector(); }
+      u_->UseDevice(true);
+      u_->SetSize(Op_Imag_->Width());
+
       Op_Imag_->MultTranspose(x_i, *u_);
-      y_r_.Add(convention_ == BLOCK_SYMMETRIC ? -1.0 : 1.0, *u_);
+      y_r.Add(convention_ == BLOCK_SYMMETRIC ? -1.0 : 1.0, *u_);
       Op_Imag_->MultTranspose(x_r, *u_);
-      y_i_ -= *u_;
+      y_i.Add(-1.0, *u_);
    }
 }
 

linalg/solvers.cpp

@@ -134,7 +134,7 @@ OperatorJacobiSmoother::OperatorJacobiSmoother(const BilinearForm &a,
 
 OperatorJacobiSmoother::OperatorJacobiSmoother(const Vector &d,
                                                const Array<int> &ess_tdofs,
-                                               const double dmpng)
+                                               const double dmpng, const bool inverse)
    :
    Solver(d.Size()),
    N(d.Size()),
@@ -143,16 +143,30 @@ OperatorJacobiSmoother::OperatorJacobiSmoother(const Vector &d,
    ess_tdof_list(ess_tdofs),
    residual(N)
 {
-   Setup(d);
+   Setup(d, inverse);
 }
 
-void OperatorJacobiSmoother::Setup(const Vector &diag)
+void OperatorJacobiSmoother::Setup(const Vector &diag, const bool inverse)
 {
    residual.UseDevice(true);
    const double delta = damping;
    auto D = diag.Read();
    auto DI = dinv.Write();
-   MFEM_FORALL(i, N, DI[i] = delta / D[i]; );
+   if (inverse)
+   {
+      if (delta > 0.0)
+      {
+         MFEM_FORALL(i, N, DI[i] = delta * D[i]; );
+      }
+      else
+      {
+         MFEM_FORALL(i, N, DI[i] = D[i]; );
+      }
+   }
+   else
+   {
+      MFEM_FORALL(i, N, DI[i] = delta / D[i]; );
+   }
    auto I = ess_tdof_list.Read();
    MFEM_FORALL(i, ess_tdof_list.Size(), DI[I[i]] = delta; );
 }

linalg/solvers.hpp

@@ -125,13 +125,14 @@ public:
        the matrix-free setting. */
    OperatorJacobiSmoother(const Vector &d,
                           const Array<int> &ess_tdof_list,
-                          const double damping=1.0);
+                          const double damping=1.0, 
+                          const bool inverse=false);
    ~OperatorJacobiSmoother() {}
 
    void Mult(const Vector &x, Vector &y) const;
    void MultTranspose(const Vector &x, Vector &y) const { Mult(x, y); }
    void SetOperator(const Operator &op) { oper = &op; }
-   void Setup(const Vector &diag);
+   void Setup(const Vector &diag, const bool inverse=false);
 
 private:
    const int N;

tests/unit/fem/test_assemblediagonalpa.cpp

@@ -17,9 +17,107 @@ using namespace mfem;
 namespace assemblediagonalpa
 {
 
+int dimension;
+
+double coeffFunction(const Vector& x)
+{
+   if (dimension == 2)
+   {
+      return sin(8.0 * M_PI * x[0]) * cos(6.0 * M_PI * x[1]) + 2.0;
+   }
+   else
+   {
+      return sin(8.0 * M_PI * x[0]) * cos(6.0 * M_PI * x[1]) *
+             sin(4.0 * M_PI * x[2]) +
+             2.0;
+   }
+}
+
+void vectorCoeffFunction(const Vector & x, Vector & f)
+{
+   f = 0.0;
+   if (dimension > 1)
+   {
+      f[0] = sin(M_PI * x[1]);
+      f[1] = sin(2.5 * M_PI * x[0]);
+   }
+   if (dimension == 3)
+   {
+      f[2] = sin(6.1 * M_PI * x[2]);
+   }
+}
+
+void asymmetricMatrixCoeffFunction(const Vector & x, DenseMatrix & f)
+{
+   f = 0.0;
+   if (dimension == 2)
+   {
+      f(0,0) = 1.1 + sin(M_PI * x[1]);  // 1,1
+      f(1,0) = cos(1.3 * M_PI * x[1]);  // 2,1
+      f(0,1) = cos(2.5 * M_PI * x[0]);  // 1,2
+      f(1,1) = 1.1 + sin(4.9 * M_PI * x[0]);  // 2,2
+   }
+   else if (dimension == 3)
+   {
+      f(0,0) = 1.1 + sin(M_PI * x[1]);  // 1,1
+      f(0,1) = cos(2.5 * M_PI * x[0]);  // 1,2
+      f(0,2) = sin(4.9 * M_PI * x[2]);  // 1,3
+      f(1,0) = cos(M_PI * x[0]);  // 2,1
+      f(1,1) = 1.1 + sin(6.1 * M_PI * x[1]);  // 2,2
+      f(1,2) = cos(6.1 * M_PI * x[2]);  // 2,3
+      f(2,0) = sin(1.5 * M_PI * x[1]);  // 3,1
+      f(2,1) = cos(2.9 * M_PI * x[0]);  // 3,2
+      f(2,2) = 1.1 + sin(6.1 * M_PI * x[2]);  // 3,3
+   }
+}
+
+void fullSymmetricMatrixCoeffFunction(const Vector & x, DenseMatrix & f)
+{
+   f = 0.0;
+   if (dimension == 2)
+   {
+      f(0,0) = 1.1 + sin(M_PI * x[1]);  // 1,1
+      f(0,1) = cos(2.5 * M_PI * x[0]);  // 1,2
+      f(1,1) = 1.1 + sin(4.9 * M_PI * x[0]);  // 2,2
+      f(1,0) = f(0,1);
+   }
+   else if (dimension == 3)
+   {
+      f(0,0) = sin(M_PI * x[1]);  // 1,1
+      f(0,1) = cos(2.5 * M_PI * x[0]);  // 1,2
+      f(0,2) = sin(4.9 * M_PI * x[2]);  // 1,3
+      f(1,1) = sin(6.1 * M_PI * x[1]);  // 2,2
+      f(1,2) = cos(6.1 * M_PI * x[2]);  // 2,3
+      f(2,2) = sin(6.1 * M_PI * x[2]);  // 3,3
+      f(1,0) = f(0,1);
+      f(2,0) = f(0,2);
+      f(2,1) = f(1,2);
+   }
+}
+
+void symmetricMatrixCoeffFunction(const Vector & x, Vector & f)
+{
+   f = 0.0;
+   if (dimension == 2)
+   {
+      f[0] = 1.1 + sin(M_PI * x[1]);  // 1,1
+      f[1] = cos(2.5 * M_PI * x[0]);  // 1,2
+      f[2] = 1.1 + sin(4.9 * M_PI * x[0]);  // 2,2
+   }
+   else if (dimension == 3)
+   {
+      f[0] = sin(M_PI * x[1]);  // 1,1
+      f[1] = cos(2.5 * M_PI * x[0]);  // 1,2
+      f[2] = sin(4.9 * M_PI * x[2]);  // 1,3
+      f[3] = sin(6.1 * M_PI * x[1]);  // 2,2
+      f[4] = cos(6.1 * M_PI * x[2]);  // 2,3
+      f[5] = sin(6.1 * M_PI * x[2]);  // 3,3
+   }
+}
+
 TEST_CASE("massdiag")
 {
-   for (int dimension = 2; dimension < 4; ++dimension)
+   for (dimension = 2; dimension < 4; ++dimension)
    {
       for (int ne = 1; ne < 3; ++ne)
       {
@@ -67,7 +165,7 @@ TEST_CASE("massdiag")
 
 TEST_CASE("diffusiondiag")
 {
-   for (int dimension = 2; dimension < 4; ++dimension)
+   for (dimension = 2; dimension < 4; ++dimension)
    {
       for (int ne = 1; ne < 3; ++ne)
       {
@@ -199,81 +297,139 @@ TEST_CASE("Vector Diffusion Diagonal PA",
 
 TEST_CASE("Hcurl/Hdiv diagonal PA")
 {
-   for (int dimension = 2; dimension < 4; ++dimension)
+   for (dimension = 2; dimension < 4; ++dimension)
    {
-      for (int spaceType = 0; spaceType < 2; ++spaceType)
-         for (int integrator = 0; integrator < 2; ++integrator)
+      for (int coeffType = 0; coeffType < 5; ++coeffType)
+      {
+         const int numSpaces = (coeffType == 0) ? 2 : 1;
+         const int numIntegrators = (coeffType == 0) ? 2 : 1;
+
+         Coefficient* coeff = nullptr;
+         VectorCoefficient* vcoeff = nullptr;
+         MatrixCoefficient* mcoeff = nullptr;
+         MatrixCoefficient* smcoeff = nullptr;
+         if (coeffType == 0)
          {
-            for (int ne = 1; ne < 3; ++ne)
+            coeff = new ConstantCoefficient(12.34);
+         }
+         else if (coeffType == 1)
+         {
+            coeff = new FunctionCoefficient(&coeffFunction);
+         }
+         else if (coeffType == 2)
+         {
+            vcoeff = new VectorFunctionCoefficient(dimension, &vectorCoeffFunction);
+         }
+         else if (coeffType == 3)
+         {
+            mcoeff = new MatrixFunctionCoefficient(dimension,
+                                                   &fullSymmetricMatrixCoeffFunction);
+            smcoeff = new MatrixFunctionCoefficient(dimension,
+                                                    &symmetricMatrixCoeffFunction);
+         }
+         else if (coeffType == 4)
+         {
+            mcoeff = new MatrixFunctionCoefficient(dimension,
+                                                   &asymmetricMatrixCoeffFunction);
+            smcoeff = new MatrixFunctionCoefficient(dimension,
+                                                    &asymmetricMatrixCoeffFunction);
+         }
+
+         for (int spaceType = 0; spaceType < numSpaces; ++spaceType)
+         {
+            for (int integrator = 0; integrator < numIntegrators; ++integrator)
             {
-               if (spaceType == 0)
-                  std::cout << "Testing " << dimension <<
-                            "D partial assembly H(curl) diagonal for integrator " << integrator << ": "
-                            << std::pow(ne, dimension) << " elements." << std::endl;
-               else
-                  std::cout << "Testing " << dimension <<
-                            "D partial assembly H(div) diagonal for integrator " << integrator << ": "
-                            << std::pow(ne, dimension) << " elements." << std::endl;
-
-               for (int order = 1; order < 4; ++order)
+               for (int ne = 1; ne < 3; ++ne)
                {
-                  Mesh * mesh;
-                  if (dimension == 2)
-                  {
-                     mesh = new Mesh(ne, ne, Element::QUADRILATERAL, 1, 1.0, 1.0);
-                  }
+                  if (spaceType == 0)
+                     std::cout << "Testing " << dimension <<
+                               "D partial assembly H(curl) diagonal for integrator " << integrator
+                               << " and coeffType " << coeffType << ": "
+                               << std::pow(ne, dimension) << " elements." << std::endl;
                   else
-                  {
-                     mesh = new Mesh(ne, ne, ne, Element::HEXAHEDRON, 1, 1.0, 1.0, 1.0);
-                  }
+                     std::cout << "Testing " << dimension <<
+                               "D partial assembly H(div) diagonal for integrator " << integrator
+                               << " and coeffType " << coeffType << ": "
+                               << std::pow(ne, dimension) << " elements." << std::endl;
 
-                  FiniteElementCollection* fec = (spaceType == 0) ?
-                                                 (FiniteElementCollection*) new ND_FECollection(order, dimension) :
-                                                 (FiniteElementCollection*) new RT_FECollection(order, dimension);
-
-                  FiniteElementSpace fespace(mesh, fec);
-                  BilinearForm paform(&fespace);
-                  BilinearForm faform(&fespace);
-                  ConstantCoefficient one(1.0);
-                  paform.SetAssemblyLevel(AssemblyLevel::PARTIAL);
-                  if (integrator == 0)
+                  for (int order = 1; order < 4; ++order)
                   {
-                     paform.AddDomainIntegrator(new VectorFEMassIntegrator(one));
-                     faform.AddDomainIntegrator(new VectorFEMassIntegrator(one));
-                  }
-                  else
-                  {
-                     if (spaceType == 0)
+                     Mesh * mesh;
+                     if (dimension == 2)
                      {
-                        paform.AddDomainIntegrator(new CurlCurlIntegrator(one));
-                        faform.AddDomainIntegrator(new CurlCurlIntegrator(one));
+                        mesh = new Mesh(ne, ne, Element::QUADRILATERAL, 1, 1.0, 1.0);
                      }
                      else
                      {
-                        paform.AddDomainIntegrator(new DivDivIntegrator(one));
-                        faform.AddDomainIntegrator(new DivDivIntegrator(one));
+                        mesh = new Mesh(ne, ne, ne, Element::HEXAHEDRON, 1, 1.0, 1.0, 1.0);
                      }
+
+                     FiniteElementCollection* fec = (spaceType == 0) ?
+                                                    (FiniteElementCollection*) new ND_FECollection(order, dimension) :
+                                                    (FiniteElementCollection*) new RT_FECollection(order, dimension);
+
+                     FiniteElementSpace fespace(mesh, fec);
+                     BilinearForm paform(&fespace);
+                     BilinearForm faform(&fespace);
+                     paform.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+                     if (integrator == 0)
+                     {
+                        if (coeffType >= 3)
+                        {
+                           paform.AddDomainIntegrator(new VectorFEMassIntegrator(*smcoeff));
+                           faform.AddDomainIntegrator(new VectorFEMassIntegrator(*mcoeff));
+                        }
+                        else if (coeffType == 2)
+                        {
+                           paform.AddDomainIntegrator(new VectorFEMassIntegrator(*vcoeff));
+                           faform.AddDomainIntegrator(new VectorFEMassIntegrator(*vcoeff));
+                        }
+                        else
+                        {
+                           paform.AddDomainIntegrator(new VectorFEMassIntegrator(*coeff));
+                           faform.AddDomainIntegrator(new VectorFEMassIntegrator(*coeff));
+                        }
+                     }
+                     else
+                     {
+                        if (spaceType == 0)
+                        {
+                           paform.AddDomainIntegrator(new CurlCurlIntegrator(*coeff));
+                           faform.AddDomainIntegrator(new CurlCurlIntegrator(*coeff));
+                        }
+                        else
+                        {
+                           paform.AddDomainIntegrator(new DivDivIntegrator(*coeff));
+                           faform.AddDomainIntegrator(new DivDivIntegrator(*coeff));
+                        }
+                     }
+                     paform.Assemble();
+                     Vector pa_diag(fespace.GetVSize());
+                     paform.AssembleDiagonal(pa_diag);
+
+                     faform.Assemble();
+                     faform.Finalize();
+                     Vector assembly_diag(fespace.GetVSize());
+                     faform.SpMat().GetDiag(assembly_diag);
+
+                     assembly_diag -= pa_diag;
+                     double error = assembly_diag.Norml2();
+                     std::cout << "    order: " << order << ", error norm: " << error << std::endl;
+                     REQUIRE(assembly_diag.Norml2() < 1.e-11);
+
+                     delete mesh;
+                     delete fec;
                   }
-                  paform.Assemble();
-                  Vector pa_diag(fespace.GetVSize());
-                  paform.AssembleDiagonal(pa_diag);
+               }  // ne
+            }  // integrator
+         }  // spaceType
 
-                  faform.Assemble();
-                  faform.Finalize();
-                  Vector assembly_diag(fespace.GetVSize());
-                  faform.SpMat().GetDiag(assembly_diag);
-
-                  assembly_diag -= pa_diag;
-                  double error = assembly_diag.Norml2();
-                  std::cout << "    order: " << order << ", error norm: " << error << std::endl;
-                  REQUIRE(assembly_diag.Norml2() < 1.e-12);
-
-                  delete mesh;
-                  delete fec;
-               }
-            }
-         }
-   }
+         delete coeff;
+         delete vcoeff;
+         delete mcoeff;
+         delete smcoeff;
+      }  // coeffType
+   }  // dimension
 }
 
 } // namespace assemblediagonalpa

tests/unit/fem/test_pa_coeff.cpp

@@ -59,6 +59,74 @@ double linearFunction(const Vector & x)
    }
 }
 
+void asymmetricMatrixCoeffFunction(const Vector & x, DenseMatrix & f)
+{
+   f = 0.0;
+   if (dimension == 2)
+   {
+      f(0,0) = 1.1 + sin(M_PI * x[1]);  // 1,1
+      f(1,0) = cos(1.3 * M_PI * x[1]);  // 2,1
+      f(0,1) = cos(2.5 * M_PI * x[0]);  // 1,2
+      f(1,1) = 1.1 + sin(4.9 * M_PI * x[0]);  // 2,2
+   }
+   else if (dimension == 3)
+   {
+      f(0,0) = 1.1 + sin(M_PI * x[1]);  // 1,1
+      f(0,1) = cos(2.5 * M_PI * x[0]);  // 1,2
+      f(0,2) = sin(4.9 * M_PI * x[2]);  // 1,3
+      f(1,0) = cos(M_PI * x[0]);  // 2,1
+      f(1,1) = 1.1 + sin(6.1 * M_PI * x[1]);  // 2,2
+      f(1,2) = cos(6.1 * M_PI * x[2]);  // 2,3
+      f(2,0) = sin(1.5 * M_PI * x[1]);  // 3,1
+      f(2,1) = cos(2.9 * M_PI * x[0]);  // 3,2
+      f(2,2) = 1.1 + sin(6.1 * M_PI * x[2]);  // 3,3
+   }
+}
+
+void fullSymmetricMatrixCoeffFunction(const Vector & x, DenseMatrix & f)
+{
+   f = 0.0;
+   if (dimension == 2)
+   {
+      f(0,0) = 1.1 + sin(M_PI * x[1]);  // 1,1
+      f(0,1) = cos(2.5 * M_PI * x[0]);  // 1,2
+      f(1,1) = 1.1 + sin(4.9 * M_PI * x[0]);  // 2,2
+      f(1,0) = f(0,1);
+   }
+   else if (dimension == 3)
+   {
+      f(0,0) = sin(M_PI * x[1]);  // 1,1
+      f(0,1) = cos(2.5 * M_PI * x[0]);  // 1,2
+      f(0,2) = sin(4.9 * M_PI * x[2]);  // 1,3
+      f(1,1) = sin(6.1 * M_PI * x[1]);  // 2,2
+      f(1,2) = cos(6.1 * M_PI * x[2]);  // 2,3
+      f(2,2) = sin(6.1 * M_PI * x[2]);  // 3,3
+      f(1,0) = f(0,1);
+      f(2,0) = f(0,2);
+      f(2,1) = f(1,2);
+   }
+}
+
+void symmetricMatrixCoeffFunction(const Vector & x, Vector & f)
+{
+   f = 0.0;
+   if (dimension == 2)
+   {
+      f[0] = 1.1 + sin(M_PI * x[1]);  // 1,1
+      f[1] = cos(2.5 * M_PI * x[0]);  // 1,2
+      f[2] = 1.1 + sin(4.9 * M_PI * x[0]);  // 2,2
+   }
+   else if (dimension == 3)
+   {
+      f[0] = sin(M_PI * x[1]);  // 1,1
+      f[1] = cos(2.5 * M_PI * x[0]);  // 1,2
+      f[2] = sin(4.9 * M_PI * x[2]);  // 1,3
+      f[3] = sin(6.1 * M_PI * x[1]);  // 2,2
+      f[4] = cos(6.1 * M_PI * x[2]);  // 2,3
+      f[5] = sin(6.1 * M_PI * x[2]);  // 3,3
+   }
+}
+
 TEST_CASE("H1 pa_coeff")
 {
    for (dimension = 2; dimension < 4; ++dimension)
@@ -185,11 +253,13 @@ TEST_CASE("Hcurl/Hdiv pa_coeff")
          mesh = new Mesh(ne, ne, ne, Element::HEXAHEDRON, 1, 1.0, 1.0, 1.0);
       }
 
-      for (int coeffType = 0; coeffType < 3; ++coeffType)
+      for (int coeffType = 0; coeffType < 5; ++coeffType)
       {
          Coefficient* coeff = nullptr;
          Coefficient* coeff2 = nullptr;
          VectorCoefficient* vcoeff = nullptr;
+         MatrixCoefficient* mcoeff = nullptr;
+         MatrixCoefficient* smcoeff = nullptr;
          if (coeffType == 0)
          {
             coeff = new ConstantCoefficient(12.34);
@@ -205,33 +275,73 @@ TEST_CASE("Hcurl/Hdiv pa_coeff")
             vcoeff = new VectorFunctionCoefficient(dimension, &vectorCoeffFunction);
             coeff2 = new FunctionCoefficient(&linearFunction);
          }
-
-         for (int spaceType = 0; spaceType < 2; ++spaceType)
+         else if (coeffType == 3)
          {
-            if (spaceType == 1 && coeffType == 2)
+            mcoeff = new MatrixFunctionCoefficient(dimension,
+                                                   &fullSymmetricMatrixCoeffFunction);
+            smcoeff = new MatrixFunctionCoefficient(dimension,
+                                                    &symmetricMatrixCoeffFunction);
+            coeff2 = new FunctionCoefficient(&linearFunction);
+         }
+         else if (coeffType == 4)
+         {
+            mcoeff = new MatrixFunctionCoefficient(dimension,
+                                                   &asymmetricMatrixCoeffFunction);
+            smcoeff = new MatrixFunctionCoefficient(dimension,
+                                                    &asymmetricMatrixCoeffFunction);
+            coeff2 = new FunctionCoefficient(&linearFunction);
+         }
+
+         enum MixedSpaces {Hcurl, Hdiv, HcurlHdiv, HdivHcurl, NumSpaceTypes};
+
+         for (int spaceType = 0; spaceType < NumSpaceTypes; ++spaceType)
+         {
+            if (spaceType == Hdiv && coeffType >= 2)
             {
                continue;   // Case not implemented yet
             }
 
-            const int numIntegrators = (coeffType == 2) ? 2 : 3;
+            const int numIntegrators =
+               (spaceType >= HcurlHdiv) ? 1 : ((coeffType == 2) ? 2 : 3);
+
             for (int integrator = 0; integrator < numIntegrators; ++integrator)
             {
-               if (spaceType == 0)
+               if (spaceType == Hcurl)
                   std::cout << "Testing " << dimension
-                            << "D ND partial assembly with " << "coeffType "
-                            << coeffType << " and " << "integrator "
+                            << "D ND partial assembly with coeffType "
+                            << coeffType << " and integrator "
                             << integrator << std::endl;
-               else
+               else if (spaceType == Hdiv)
                   std::cout << "Testing " << dimension
-                            << "D RT partial assembly with " << "coeffType "
-                            << coeffType << " and " << "integrator "
+                            << "D RT partial assembly with coeffType "
+                            << coeffType << " and integrator "
+                            << integrator << std::endl;
+               else if (spaceType == HcurlHdiv)
+                  std::cout << "Testing " << dimension
+                            << "D ND x RT partial assembly with coeffType "
+                            << coeffType << " and integrator "
+                            << integrator << std::endl;
+               else  // HdivHcurl
+                  std::cout << "Testing " << dimension
+                            << "D RT x ND partial assembly with coeffType "
+                            << coeffType << " and integrator "
                             << integrator << std::endl;
 
                for (int order = 1; order < 4; ++order)
                {
-                  FiniteElementCollection* fec = (spaceType == 0) ?
-                                                 (FiniteElementCollection*) new ND_FECollection(order, dimension) :
-                                                 (FiniteElementCollection*) new RT_FECollection(order, dimension);
+                  FiniteElementCollection* fec = nullptr;
+                  if (spaceType == Hcurl || spaceType == HcurlHdiv)
+                  {
+                     fec = (FiniteElementCollection*) new ND_FECollection(order, dimension);
+                  }
+                  else if (spaceType == HdivHcurl)
+                  {
+                     fec = (FiniteElementCollection*) new RT_FECollection(order - 1, dimension);
+                  }
+                  else
+                  {
+                     fec = (FiniteElementCollection*) new RT_FECollection(order, dimension);
+                  }
 
                   FiniteElementSpace fespace(mesh, fec);
 
@@ -270,59 +380,144 @@ TEST_CASE("Hcurl/Hdiv pa_coeff")
                      }
                   }
 
-                  BilinearForm paform(&fespace);
-                  paform.SetAssemblyLevel(AssemblyLevel::PARTIAL);
-                  BilinearForm assemblyform(&fespace);
-                  if (integrator < 2)
+                  Vector xin(fespace.GetTrueVSize());
+                  xin.Randomize();
+
+                  Vector y_mat, y_assembly, y_pa;
+
+                  if (spaceType >= HcurlHdiv)
                   {
-                     if (coeffType == 2)
+                     FiniteElementCollection* fecTest = nullptr;
+                     if (spaceType == HcurlHdiv)
+                     {
+                        fecTest = (FiniteElementCollection*) new RT_FECollection(order - 1, dimension);
+                     }
+                     else
+                     {
+                        fecTest = (FiniteElementCollection*) new ND_FECollection(order, dimension);
+                     }
+
+                     FiniteElementSpace fespaceTest(mesh, fecTest);
+
+                     MixedBilinearForm paform(&fespace, &fespaceTest);
+                     paform.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+                     MixedBilinearForm assemblyform(&fespace, &fespaceTest);
+
+                     const int testSize = fespaceTest.GetTrueVSize();
+                     y_mat.SetSize(testSize);
+                     y_mat = 0.0;
+                     y_assembly.SetSize(testSize);
+                     y_assembly = 0.0;
+                     y_pa.SetSize(testSize);
+                     y_pa = 0.0;
+
+                     if (coeffType >= 3)
+                     {
+                        paform.AddDomainIntegrator(new VectorFEMassIntegrator(*smcoeff));
+                        assemblyform.AddDomainIntegrator(new VectorFEMassIntegrator(*mcoeff));
+                     }
+                     else if (coeffType == 2)
                      {
                         paform.AddDomainIntegrator(new VectorFEMassIntegrator(*vcoeff));
-                        assemblyform.AddDomainIntegrator(
-                           new VectorFEMassIntegrator(*vcoeff));
+                        assemblyform.AddDomainIntegrator(new VectorFEMassIntegrator(*vcoeff));
                      }
                      else
                      {
                         paform.AddDomainIntegrator(new VectorFEMassIntegrator(*coeff));
-                        assemblyform.AddDomainIntegrator(
-                           new VectorFEMassIntegrator(*coeff));
+                        assemblyform.AddDomainIntegrator(new VectorFEMassIntegrator(*coeff));
                      }
+
+                     Array<int> empty_ess; // empty
+
+                     paform.Assemble();
+                     OperatorHandle paopr;
+                     paform.FormRectangularSystemMatrix(ess_tdof_list, empty_ess, paopr);
+
+                     assemblyform.Assemble();
+                     SparseMatrix A_explicit;
+                     assemblyform.FormRectangularSystemMatrix(ess_tdof_list, empty_ess, A_explicit);
+
+                     paopr->Mult(xin, y_pa);
+                     assemblyform.Mult(xin, y_assembly);
+                     A_explicit.Mult(xin, y_mat);
+
+                     delete fecTest;
                   }
-                  if (integrator > 0)
+                  else
                   {
-                     if (spaceType == 0)
+                     BilinearForm paform(&fespace);
+                     paform.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+                     BilinearForm assemblyform(&fespace);
+
+                     y_mat.SetSize(xin.Size());
+                     y_mat = 0.0;
+                     y_assembly.SetSize(xin.Size());
+                     y_assembly = 0.0;
+                     y_pa.SetSize(xin.Size());
+                     y_pa = 0.0;
+
+                     if (integrator < 2)
                      {
-                        paform.AddDomainIntegrator(new CurlCurlIntegrator(*coeff2));
-                        assemblyform.AddDomainIntegrator(new CurlCurlIntegrator(*coeff2));
+                        if (coeffType >= 3)
+                        {
+                           paform.AddDomainIntegrator(new VectorFEMassIntegrator(*smcoeff));
+                           assemblyform.AddDomainIntegrator(new VectorFEMassIntegrator(*mcoeff));
+                        }
+                        else if (coeffType == 2)
+                        {
+                           paform.AddDomainIntegrator(new VectorFEMassIntegrator(*vcoeff));
+                           assemblyform.AddDomainIntegrator(new VectorFEMassIntegrator(*vcoeff));
+
+                        }
+                        else
+                        {
+                           paform.AddDomainIntegrator(new VectorFEMassIntegrator(*coeff));
+                           assemblyform.AddDomainIntegrator(new VectorFEMassIntegrator(*coeff));
+                        }
                      }
-                     else
+                     if (integrator > 0)
                      {
-                        paform.AddDomainIntegrator(new DivDivIntegrator(*coeff2));
-                        assemblyform.AddDomainIntegrator(new DivDivIntegrator(*coeff2));
+                        if (spaceType == Hcurl)
+                        {
+                           const FiniteElement *fel = fespace.GetFE(0);
+                           const IntegrationRule *intRule = &MassIntegrator::GetRule(*fel, *fel,
+                                                                                     *mesh->GetElementTransformation(0));
+
+                           if (coeffType >= 3 && dimension == 3)
+                           {
+                              paform.AddDomainIntegrator(new CurlCurlIntegrator(*smcoeff, intRule));
+                              assemblyform.AddDomainIntegrator(new CurlCurlIntegrator(*mcoeff, intRule));
+                           }
+                           else if (coeffType == 2 && dimension == 3)
+                           {
+                              paform.AddDomainIntegrator(new CurlCurlIntegrator(*vcoeff, intRule));
+                              assemblyform.AddDomainIntegrator(new CurlCurlIntegrator(*vcoeff, intRule));
+                           }
+                           else
+                           {
+                              paform.AddDomainIntegrator(new CurlCurlIntegrator(*coeff2));
+                              assemblyform.AddDomainIntegrator(new CurlCurlIntegrator(*coeff2));
+                           }
+                        }
+                        else
+                        {
+                           paform.AddDomainIntegrator(new DivDivIntegrator(*coeff2));
+                           assemblyform.AddDomainIntegrator(new DivDivIntegrator(*coeff2));
+                        }
                      }
+                     paform.Assemble();
+                     OperatorHandle paopr;
+                     paform.FormSystemMatrix(ess_tdof_list, paopr);
+
+                     assemblyform.SetDiagonalPolicy(Matrix::DIAG_ONE);
+                     assemblyform.Assemble();
+                     SparseMatrix A_explicit;
+                     assemblyform.FormSystemMatrix(ess_tdof_list, A_explicit);
+
+                     paopr->Mult(xin, y_pa);
+                     assemblyform.Mult(xin, y_assembly);
+                     A_explicit.Mult(xin, y_mat);
                   }
-                  paform.Assemble();
-                  OperatorHandle paopr;
-                  paform.FormSystemMatrix(ess_tdof_list, paopr);
-
-                  assemblyform.SetDiagonalPolicy(Matrix::DIAG_ONE);
-                  assemblyform.Assemble();
-                  assemblyform.Finalize();
-                  SparseMatrix A_explicit;
-                  assemblyform.FormSystemMatrix(ess_tdof_list, A_explicit);
-
-                  Vector xin(fespace.GetTrueVSize());
-                  xin.Randomize();
-                  Vector y_mat(xin);
-                  y_mat = 0.0;
-                  Vector y_assembly(xin);
-                  y_assembly = 0.0;
-                  Vector y_pa(xin);
-                  y_pa = 0.0;
-
-                  paopr->Mult(xin, y_pa);
-                  assemblyform.Mult(xin, y_assembly);
-                  A_explicit.Mult(xin, y_mat);
 
                   y_pa -= y_mat;
                   double pa_error = y_pa.Norml2();
@@ -344,6 +539,9 @@ TEST_CASE("Hcurl/Hdiv pa_coeff")
 
          delete coeff;
          delete coeff2;
+         delete vcoeff;
+         delete mcoeff;
+         delete smcoeff;
       }
 
       delete mesh;
@@ -382,9 +580,9 @@ TEST_CASE("Hcurl/Hdiv mixed pa_coeff")
             vcoeff = new VectorFunctionCoefficient(dimension, &vectorCoeffFunction);
          }
 
-         enum MixedSpaces {HcurlH1, HcurlL2, HdivL2};
+         enum MixedSpaces {HcurlH1, HcurlL2, HdivL2, NumSpaceTypes};
 
-         for (int spaceType = 0; spaceType < 3; ++spaceType)
+         for (int spaceType = 0; spaceType < NumSpaceTypes; ++spaceType)
          {
             if (spaceType == HdivL2 && coeffType == 1)
             {