complex-jacobi minor cleanup

Merging master into feature branch.

.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/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
@@ -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)
@@ -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)
@@ -430,39 +449,58 @@ int main(int argc, char *argv[])
       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;

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
@@ -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)
@@ -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);
@@ -464,11 +483,9 @@ int main(int argc, char *argv[])
       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

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;