Merge remote-tracking branch 'origin/eigensolvers-dev' into arpack-dev

# Conflicts:
#	linalg/eigensolvers.hpp

.gitlab/jobs/dane.yml

@@ -85,8 +85,3 @@ opt_par_gcc_10_pumi:
   extends: .mfem_job_on_dane
   variables:
     SPEC: "%gcc@10.3.1 +pumi"
-
-opt_par_gcc_10_gslib:
-  extends: .mfem_job_on_dane
-  variables:
-    SPEC: "%gcc@10.3.1 +gslib"

.gitlab/jobs/matrix.yml

@@ -63,8 +63,3 @@ opt_mpi_cuda_hypre_cuda_gcc:
   extends: .mfem_job_on_matrix
   variables:
     SPEC: "%gcc@10.3.1 +mpi +cuda cuda_arch=90 ^hypre+cuda"
-
-opt_mpi_cuda_gcc_gslib:
-  extends: .mfem_job_on_matrix
-  variables:
-    SPEC: "%gcc@10.3.1 +mpi +cuda +gslib cuda_arch=90 ^hypre+cuda"

.gitlab/scripts/baseline

@@ -32,9 +32,9 @@ mkdir _${BASELINE_TEST} && cd _${BASELINE_TEST}
 
 # run
 if [[ "${MACHINE_NAME}" == "dane" ]]; then
-  srun --nodes=1 -t 60 --exclusive --reservation=ci ../runtest ../../mfem "${BASELINE_TEST} ${TPLS_DIR}"
+  salloc --nodes=1 -t 60 --exclusive --reservation=ci ../runtest ../../mfem "${BASELINE_TEST} ${TPLS_DIR}"
 elif [[ ${MACHINE_NAME} == "corona" ]]; then
-  srun --nodes=1 -t 60 -p pbatch ../runtest ../../mfem "${BASELINE_TEST} ${TPLS_DIR}"
+  salloc --nodes=1 -t 60 -p pbatch ../runtest ../../mfem "${BASELINE_TEST} ${TPLS_DIR}"
 else
   echo "Unknown machine: MACHINE_NAME=$MACHINE_NAME"
   exit 1

CHANGELOG

@@ -8,10 +8,8 @@
                                https://mfem.org
 
 
-Version 4.9.1 (development)
-===========================
-
-- Policy for AI-assisted contribution added to CONTRIBUTING.md
+Version 4.10 (development)
+==========================
 
 Discretization improvements
 ---------------------------
@@ -25,6 +23,12 @@ Discretization improvements
   Tet rules (d=14-20): Chuluunbaatar et al., Comput. Math. Appl. 124:89-97,
   2022.
 
+
+Version 4.9.1 (development)
+===========================
+
+Discretization improvements
+---------------------------
 - Improved the gridfunction projection routines. Projections work for Scalar,
   Vector and VectorFE, also NURBS versions. Optionally different types of
   projections can be selected, default behaviour has not changed.

CONTRIBUTING.md

@@ -24,14 +24,6 @@ must be made under this license.
 Note also that MFEM has a [Code of Conduct](CODE_OF_CONDUCT.md). By participating
 in the MFEM community, you agree to abide by its rules.
 
-## AI Policy
-- Use of AI code generation in MFEM is allowed but must be disclosed, e.g. by
-  selecting the `AI-assisted` label on the PR.
-- By submitting a PR, the author acknowledges that they have reviewed and
-  understand the changes they are proposing.
-- PR authors are still responsible for correctness, licensing, and attribution
-  of all changes.
-
 If you plan on contributing to MFEM, consider reviewing the
 [issue tracker](https://github.com/mfem/mfem/issues) first to check if a thread
 already exists for your desired feature or the bug you ran into. Use a pull

config/config.hpp.in

@@ -97,6 +97,9 @@
 // Enable MFEM functionality based on the SuiteSparse library.
 // #define MFEM_USE_SUITESPARSE
 
+// Enable MFEM functionality based on the ARPACK library.
+// #define MFEM_USE_ARPACK
+
 // Enable MFEM functionality based on the SuperLU_DIST library.
 // #define MFEM_USE_SUPERLU
 // #define MFEM_USE_SUPERLU5

config/config.mk.in

@@ -32,6 +32,7 @@ MFEM_USE_MEMALLOC      = @MFEM_USE_MEMALLOC@
 MFEM_TIMER_TYPE        = @MFEM_TIMER_TYPE@
 MFEM_USE_SUNDIALS      = @MFEM_USE_SUNDIALS@
 MFEM_USE_SUITESPARSE   = @MFEM_USE_SUITESPARSE@
+MFEM_USE_ARPACK        = @MFEM_USE_ARPACK@
 MFEM_USE_SUPERLU       = @MFEM_USE_SUPERLU@
 MFEM_USE_SUPERLU5      = @MFEM_USE_SUPERLU5@
 MFEM_USE_MUMPS         = @MFEM_USE_MUMPS@

config/defaults.mk

@@ -178,6 +178,7 @@ MFEM_USE_ALGOIM        = NO
 MFEM_USE_UMPIRE        = NO
 MFEM_USE_SIMD          = NO
 MFEM_USE_ADIOS2        = NO
+MFEM_USE_ARPACK        = NO
 MFEM_USE_MKL_CPARDISO  = NO
 MFEM_USE_MKL_PARDISO   = NO
 MFEM_USE_MOONOLITH     = NO
@@ -427,6 +428,14 @@ NETCDF_LIB = $(XLINKER)-rpath,$(NETCDF_DIR)/lib -L$(NETCDF_DIR)/lib\
  $(XLINKER)-rpath,$(HDF5_DIR)/lib -L$(HDF5_DIR)/lib\
  -lnetcdf -lhdf5_hl -lhdf5 $(ZLIB_LIB)
 
+# ARPACK library configuration
+ARPACK_DIR = @MFEM_DIR@/../ARPACK
+ifeq ($(MFEM_USE_MPI),YES)
+ARPACK_LIB = -L$(ARPACK_DIR) -lparpack -larpack
+else
+ARPACK_LIB = -L$(ARPACK_DIR) -larpack
+endif
+
 # PETSc library configuration (version greater or equal to 3.8 or the dev branch)
 PETSC_ARCH := arch-linux2-c-debug
 PETSC_DIR  := $(MFEM_DIR)/../petsc/$(PETSC_ARCH)

examples/CMakeLists.txt

@@ -49,6 +49,13 @@ list(APPEND ALL_EXE_SRCS
   ex41.cpp
   )
 
+if (MFEM_USE_ARPACK)
+list(APPEND ALL_EXE_SRCS
+  ex11.pp
+  ex13.pp
+)
+endif()
+
 if (MFEM_USE_MPI)
   list(APPEND ALL_EXE_SRCS
     ex0p.cpp
@@ -215,7 +222,7 @@ if (MFEM_ENABLE_TESTING)
     add_test(NAME ex1p_ceed_np=${MFEM_MPI_NP}
       COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
       ${MPIEXEC_PREFLAGS}
-      $<TARGET_FILE:ex1p> "-no-vis" "-d" "ceed-cpu" "-pa" "-a"
+      $<TARGET_FILE:ex1p> "-no-vis" "-d ceed-cpu" "-pa" "-a"
       ${MPIEXEC_POSTFLAGS})
   endif()
 endif()

examples/amgx/makefile

@@ -64,7 +64,7 @@ PARALLEL_NAME := Parallel AMGX example
 $(MFEM_LIB_FILE):
 	$(error The MFEM library is not build)
 
-clean: clean-build clean-exec
+clean: clean-build
 
 clean-build:
 	rm -f *.o *~ $(SEQ_EXAMPLES) $(PAR_EXAMPLES)

examples/caliper/makefile

@@ -64,12 +64,12 @@ ex1p-test-par: ex1p
 $(MFEM_LIB_FILE):
 	$(error The MFEM library is not built)
 
-clean: clean-build clean-exec
+clean: clean-build clean-exec $(SUBDIRS_CLEAN)
 
 clean-build:
 	rm -f *.o *~ $(SEQ_EXAMPLES) $(PAR_EXAMPLES)
 	rm -rf *.dSYM *.TVD.*breakpoints
 
 clean-exec:
-	@rm -f refined.mesh mesh.*
-	@rm -f sol.*
+	@rm -f refined.mesh displaced.mesh mesh.* ex5.mesh
+	@rm -f sphere_refined.* sol.* sol_u.* sol_p.* sol_r.* sol_i.*

examples/ex11.cpp

@@ -0,0 +1,298 @@
+//                       MFEM Example 11 - Serial Version
+//
+// Compile with: make ex11
+//
+// Sample runs:  ex11 -m ../data/square-disc.mesh
+//               ex11 -m ../data/star.mesh
+//               ex11 -m ../data/star-mixed.mesh
+//               ex11 -m ../data/periodic-annulus-sector.msh
+//               ex11 -m ../data/square-disc-p2.vtk -o 2
+//               ex11 -m ../data/square-disc-p3.mesh -o 3
+//               ex11 -m ../data/square-disc-nurbs.mesh -o -1
+//               ex11 -m ../data/disc-nurbs.mesh -o -1 -n 20
+//               ex11 -m ../data/star-surf.mesh
+//               ex11 -m ../data/square-disc-surf.mesh
+//               ex11 -m ../data/inline-segment.mesh
+//               ex11 -m ../data/inline-quad.mesh
+//               ex11 -m ../data/inline-tri.mesh
+//               ex11 -m ../data/amr-quad.mesh
+//               ex11 -m ../data/amr-hex.mesh
+//               ex11 -m ../data/mobius-strip.mesh -n 8
+//
+// Description:  This example code demonstrates the use of MFEM to solve the
+//               eigenvalue problem -Delta u = lambda u with homogeneous
+//               Dirichlet boundary conditions.
+//
+//               We compute a number of the lowest eigenmodes by discretizing
+//               the Laplacian and Mass operators using a FE space of the
+//               specified order, or an isoparametric/isogeometric space if
+//               order < 1 (quadratic for quadratic curvilinear mesh, NURBS for
+//               NURBS mesh, etc.)
+//
+//               The example highlights the use of the ARPACK eigenvalue solver
+//               (regular inverse mode). Reusing a single GLVis visualization
+//               window for multiple eigenfunctions is also illustrated.
+//
+//               We recommend viewing Example 1 before viewing this example.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+#ifdef MFEM_USE_ARPACK
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   const char *mesh_file = "../data/star.mesh";
+   int ser_ref_levels = 3;
+   int order = 1;
+   int nev = 5;
+   double dbc_eig = 1e3;
+   bool visualization = 1;
+   bool arp_solver = true;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
+                  "Number of times to refine the mesh uniformly in serial.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   args.AddOption(&nev, "-n", "--num-eigs",
+                  "Number of desired eigenmodes.");
+   args.AddOption(&dbc_eig, "-d", "--dbc-eig",
+                  "Eigenvalues associated with Dirichlet BC "
+                  "(should be larger than the maximum desired eigenvalue).");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   // 2. 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;
+   ifstream imesh(mesh_file);
+   if (!imesh)
+   {
+      cerr << "\nCan not open mesh file: " << mesh_file << '\n' << endl;
+      return 2;
+   }
+   mesh = new Mesh(imesh, 1, 1);
+   imesh.close();
+   int dim = mesh->Dimension();
+
+   // 3. Refine the serial mesh on all processors to increase the resolution. In
+   //    this example we do 'ref_levels' of uniform refinement (2 by default, or
+   //    specified on the command line with -rs).
+   for (int lev = 0; lev < ser_ref_levels; lev++)
+   {
+      mesh->UniformRefinement();
+   }
+
+   // 4. Define a finite element space on the mesh. Here we
+   //    use continuous Lagrange finite elements of the specified order. If
+   //    order < 1, we instead use an isoparametric/isogeometric space.
+   FiniteElementCollection *fec;
+   if (order > 0)
+   {
+      fec = new H1_FECollection(order, dim);
+   }
+   else if (mesh->GetNodes())
+   {
+      fec = mesh->GetNodes()->OwnFEC();
+   }
+   else
+   {
+      fec = new H1_FECollection(order = 1, dim);
+   }
+   FiniteElementSpace *fespace = new FiniteElementSpace(mesh, fec);
+   int size = fespace->GetVSize();
+
+   cout << "Number of unknowns: " << size << endl;
+
+   // 5. Set up the parallel bilinear forms a(.,.) and m(.,.) on the finite
+   //    element space. The first corresponds to the Laplacian operator -Delta,
+   //    while the second is a simple mass matrix needed on the right hand side
+   //    of the generalized eigenvalue problem below. The boundary conditions
+   //    are implemented by elimination with special values on the diagonal to
+   //    shift the Dirichlet eigenvalues out of the computational range. After
+   //    serial and parallel assembly we extract the corresponding parallel
+   //    matrices A and M.
+   ConstantCoefficient one(1.0);
+   Array<int> ess_bdr;
+   if (mesh->bdr_attributes.Size())
+   {
+      ess_bdr.SetSize(mesh->bdr_attributes.Max());
+      ess_bdr = 1;
+   }
+
+   BilinearForm *a = new BilinearForm(fespace);
+   a->AddDomainIntegrator(new DiffusionIntegrator(one));
+   if (mesh->bdr_attributes.Size() == 0)
+   {
+      // Add a mass term if the mesh has no boundary, e.g. periodic mesh or
+      // closed surface.
+      a->AddDomainIntegrator(new MassIntegrator(one));
+   }
+   a->Assemble();
+   if (mesh->bdr_attributes.Size() != 0)
+   {
+      a->EliminateEssentialBCDiag(ess_bdr, dbc_eig);
+   }
+   a->Finalize();
+
+   BilinearForm *m = new BilinearForm(fespace);
+   m->AddDomainIntegrator(new MassIntegrator(one));
+   m->Assemble();
+   if (mesh->bdr_attributes.Size() != 0)
+   {
+      // shift the eigenvalue corresponding to eliminated dofs to a large value
+      m->EliminateEssentialBCDiag(ess_bdr, 1.0);
+   }
+   m->Finalize();
+
+   Solver       * solver = NULL;
+
+#ifndef MFEM_USE_SUITESPARSE
+   // 6. Define a simple symmetric Gauss-Seidel preconditioner and use it to
+   //    solve the system A X = B with PCG.
+   cout << "Building CGSolver" << endl;
+   GSSmoother M(m->SpMat());
+   CGSolver * cg_solver = new CGSolver;
+   cg_solver->SetPreconditioner(M);
+   cg_solver->SetRelTol(1.0e-12);
+   solver = cg_solver;
+#else
+   // 7. If MFEM was compiled with SuiteSparse, use UMFPACK to solve the system.
+   cout << "Building UMFPackSolver" << endl;
+   UMFPackSolver * umf_solver = new UMFPackSolver;
+   umf_solver->Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
+   solver = umf_solver;
+#endif
+   solver->SetOperator(m->SpMat());
+
+   // 7. Define and configure the ARPACK eigensolver
+   SymGenEigensolver * eig_solver = NULL;
+
+   if (arp_solver)
+   {
+      // ArPackSymGen * arpack = new ArPackSymGen();
+      ArPackSAUPD * arpack = new ArPackSAUPD();
+      arpack->SetMode(2);
+      arpack->SetNumModes(nev);
+      arpack->SetMaxIter(400);
+      arpack->SetTol(1e-8);
+      arpack->SetPrintLevel(2);
+      arpack->SetSolver(*solver);
+
+      eig_solver = arpack;
+   }
+
+   eig_solver->SetOperators(*a, *m);
+
+   // 8. Compute the eigenmodes and extract the array of eigenvalues. Define a
+   //    parallel grid function to represent each of the eigenmodes returned by
+   //    the solver.
+   Array<double> eigenvalues;
+   eig_solver->Solve();
+   eig_solver->GetEigenvalues(eigenvalues);
+
+   cout << endl;
+   std::ios::fmtflags old_fmt = cout.flags();
+   cout.setf(std::ios::scientific);
+   std::streamsize old_prec = cout.precision(14);
+   for (int i=0; i<nev; i++)
+   {
+      cout << "Eigenvalue lambda   " << eigenvalues[i] << endl;
+   }
+   cout.precision(old_prec);
+   cout.flags(old_fmt);
+   cout << endl;
+
+   GridFunction x(fespace);
+
+   // 9. Save the refined mesh and the modes in parallel. This output can be
+   //    viewed later using GLVis: "glvis -np <np> -m mesh -g mode".
+   {
+      ostringstream mesh_name, mode_name;
+      mesh_name << "ex11.mesh";
+
+      ofstream mesh_ofs(mesh_name.str().c_str());
+      mesh_ofs.precision(8);
+      mesh->Print(mesh_ofs);
+
+      for (int i=0; i<nev; i++)
+      {
+         // convert eigenvector from Vector to GridFunction
+         x = eig_solver->GetEigenvector(i);
+
+         mode_name << "mode_" << setfill('0') << setw(2) << i;
+
+         ofstream mode_ofs(mode_name.str().c_str());
+         mode_ofs.precision(8);
+         x.Save(mode_ofs);
+         mode_name.str("");
+      }
+   }
+
+   // 10. Send the solution by socket to a GLVis server.
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream mode_sock(vishost, visport);
+      mode_sock.precision(8);
+
+      for (int i=0; i<nev; i++)
+      {
+         cout << "Eigenmode " << i+1 << '/' << nev
+              << ", Lambda = " << eigenvalues[i] << endl;
+
+         // convert eigenvector from Vector to GridFunction
+         x = eig_solver->GetEigenvector(i);
+
+         mode_sock << "solution\n" << *mesh << x << flush
+                   << "window_title 'Eigenmode " << i+1 << '/' << nev
+                   << ", Lambda = " << eigenvalues[i] << "'" << endl;
+
+         char c;
+         cout << "press (q)uit or (c)ontinue --> " << flush;
+         cin >> c;
+
+         if (c != 'c')
+         {
+            break;
+         }
+      }
+      mode_sock.close();
+   }
+
+   // 11. Free the used memory.
+   delete eig_solver;
+   delete solver;
+   delete m;
+   delete a;
+
+   delete fespace;
+   if (order > 0)
+   {
+      delete fec;
+   }
+   delete mesh;
+
+   return 0;
+}
+
+#endif // MFEM_USE_ARPACK

examples/ex11p.cpp

@@ -72,6 +72,8 @@ int main(int argc, char *argv[])
    int seed = 75;
    bool slu_solver  = false;
    bool sp_solver = false;
+   bool lob_solver = true;
+   bool arp_solver = false;
    bool cpardiso_solver = false;
    bool visualization = 1;
 
@@ -97,6 +99,10 @@ int main(int argc, char *argv[])
    args.AddOption(&sp_solver, "-sp", "--strumpack", "-no-sp",
                   "--no-strumpack", "Use the STRUMPACK Solver.");
 #endif
+#ifdef MFEM_USE_ARPACK
+   args.AddOption(&arp_solver, "-arp", "--arpack", "-no-arp",
+                  "--no-arpack", "Use the Parallel ARPACK Solver.");
+#endif
 #ifdef MFEM_USE_MKL_CPARDISO
    args.AddOption(&cpardiso_solver, "-cpardiso", "--cpardiso", "-no-cpardiso",
                   "--no-cpardiso", "Use the MKL CPardiso Solver.");
@@ -113,6 +119,11 @@ int main(int argc, char *argv[])
               << "         Defaulting to SuperLU." << endl;
       sp_solver = false;
    }
+   if (arp_solver)
+   {
+      lob_solver = false;
+   }
+
    // The command line options are also passed to the STRUMPACK
    // solver. So do not exit if some options are not recognized.
    if (!sp_solver)
@@ -243,70 +254,119 @@ int main(int argc, char *argv[])
    // 8. Define and configure the LOBPCG eigensolver and the BoomerAMG
    //    preconditioner for A to be used within the solver. Set the matrices
    //    which define the generalized eigenproblem A x = lambda M x.
+   Solver * solver = NULL;
    Solver * precond = NULL;
    if (!slu_solver && !sp_solver && !cpardiso_solver)
    {
       HypreBoomerAMG * amg = new HypreBoomerAMG(*A);
       amg->SetPrintLevel(0);
       precond = amg;
-   }
-   else
-   {
-#ifdef MFEM_USE_SUPERLU
-      if (slu_solver)
+
+      if (arp_solver)
+      {
+         HyprePCG * pcg = new HyprePCG(*A);
+         pcg->SetTol(1e-12);
+         pcg->SetPreconditioner(*amg);
+         solver = pcg;
+      }
+   }
+#ifdef MFEM_USE_SUPERLU
+   else if (slu_solver)
+   {
+      SuperLUSolver * superlu = new SuperLUSolver(MPI_COMM_WORLD);
+      superlu->SetPrintStatistics(false);
+      superlu->SetSymmetricPattern(true);
+      superlu->SetColumnPermutation(superlu::PARMETIS);
+      superlu->SetOperator(*Arow);
+
+      if (arp_solver)
+      {
+         solver = superlu;
+      }
+      else
       {
-         SuperLUSolver * superlu = new SuperLUSolver(MPI_COMM_WORLD);
-         superlu->SetPrintStatistics(false);
-         superlu->SetSymmetricPattern(true);
-         superlu->SetColumnPermutation(superlu::PARMETIS);
-         superlu->SetOperator(*Arow);
          precond = superlu;
       }
+   }
 #endif
 #ifdef MFEM_USE_STRUMPACK
-      if (sp_solver)
+   else if (sp_solver)
+   {
+      STRUMPACKSolver * strumpack = new STRUMPACKSolver(argc, argv,
+                                                        MPI_COMM_WORLD);
+      strumpack->SetPrintFactorStatistics(true);
+      strumpack->SetPrintSolveStatistics(false);
+      strumpack->SetKrylovSolver(strumpack::KrylovSolver::DIRECT);
+      strumpack->SetReorderingStrategy(strumpack::ReorderingStrategy::METIS);
+      strumpack->SetMatching(strumpack::MatchingJob::NONE);
+      strumpack->SetCompression(strumpack::CompressionType::NONE);
+      strumpack->SetOperator(*Arow);
+      strumpack->SetFromCommandLine();
+      if (arp_solver)
+      {
+         solver = strumpack;
+      }
+      else
       {
-         STRUMPACKSolver * strumpack = new STRUMPACKSolver(MPI_COMM_WORLD, argc, argv);
-         strumpack->SetPrintFactorStatistics(true);
-         strumpack->SetPrintSolveStatistics(false);
-         strumpack->SetKrylovSolver(strumpack::KrylovSolver::DIRECT);
-         strumpack->SetReorderingStrategy(strumpack::ReorderingStrategy::METIS);
-         strumpack->SetMatching(strumpack::MatchingJob::NONE);
-         strumpack->SetCompression(strumpack::CompressionType::NONE);
-         strumpack->SetOperator(*Arow);
-         strumpack->SetFromCommandLine();
          precond = strumpack;
       }
+   }
 #endif
 #ifdef MFEM_USE_MKL_CPARDISO
-      if (cpardiso_solver)
+   else if (cpardiso_solver)
+   {
+      auto cpardiso = new CPardisoSolver(A->GetComm());
+      cpardiso->SetMatrixType(CPardisoSolver::MatType::REAL_STRUCTURE_SYMMETRIC);
+      cpardiso->SetPrintLevel(1);
+      cpardiso->SetOperator(*A);
+      if (arp_solver)
+      {
+         solver = cpardiso;
+      }
+      else
       {
-         auto cpardiso = new CPardisoSolver(A->GetComm());
-         cpardiso->SetMatrixType(CPardisoSolver::MatType::REAL_STRUCTURE_SYMMETRIC);
-         cpardiso->SetPrintLevel(1);
-         cpardiso->SetOperator(*A);
          precond = cpardiso;
       }
-#endif
    }
+#endif
 
-   HypreLOBPCG * lobpcg = new HypreLOBPCG(MPI_COMM_WORLD);
-   lobpcg->SetNumModes(nev);
-   lobpcg->SetRandomSeed(seed);
-   lobpcg->SetPreconditioner(*precond);
-   lobpcg->SetMaxIter(200);
-   lobpcg->SetTol(1e-8);
-   lobpcg->SetPrecondUsageMode(1);
-   lobpcg->SetPrintLevel(1);
-   lobpcg->SetMassMatrix(*M);
-   lobpcg->SetOperator(*A);
+   SymGenEigensolver * eig_solver = NULL;
+
+   if (lob_solver)
+   {
+      HypreLOBPCG * lobpcg = new HypreLOBPCG(MPI_COMM_WORLD);
+      lobpcg->SetNumModes(nev);
+      lobpcg->SetRandomSeed(seed);
+      lobpcg->SetPreconditioner(*precond);
+      lobpcg->SetMaxIter(200);
+      lobpcg->SetTol(1e-8);
+      lobpcg->SetPrecondUsageMode(1);
+      lobpcg->SetPrintLevel(1);
+
+      eig_solver = lobpcg;
+   }
+#ifdef MFEM_USE_ARPACK
+   else if (arp_solver)
+   {
+      ArPackPSAUPD * arpack = new ArPackPSAUPD(MPI_COMM_WORLD);
+      arpack->SetNumModes(nev);
+      arpack->SetMaxIter(400);
+      arpack->SetTol(1e-8);
+      arpack->SetMode(3);
+      arpack->SetPrintLevel(2);
+      arpack->SetSolver(*solver);
+
+      eig_solver = arpack;
+   }
+#endif
+   eig_solver->SetOperators(*A, *M);
 
    // 9. Compute the eigenmodes and extract the array of eigenvalues. Define a
    //    parallel grid function to represent each of the eigenmodes returned by
    //    the solver.
    Array<real_t> eigenvalues;
-   lobpcg->Solve();
-   lobpcg->GetEigenvalues(eigenvalues);
+   eig_solver->Solve();
+   eig_solver->GetEigenvalues(eigenvalues);
    ParGridFunction x(fespace);
 
    // 10. Save the refined mesh and the modes in parallel. This output can be
@@ -321,8 +381,8 @@ int main(int argc, char *argv[])
 
       for (int i=0; i<nev; i++)
       {
-         // convert eigenvector from HypreParVector to ParGridFunction
-         x = lobpcg->GetEigenvector(i);
+         // convert eigenvector from Vector to ParGridFunction
+         x.Distribute(eig_solver->GetEigenvector(i));
 
          mode_name << "mode_" << setfill('0') << setw(2) << i << "."
                    << setfill('0') << setw(6) << myid;
@@ -350,8 +410,8 @@ int main(int argc, char *argv[])
                  << ", Lambda = " << eigenvalues[i] << endl;
          }
 
-         // convert eigenvector from HypreParVector to ParGridFunction
-         x = lobpcg->GetEigenvector(i);
+         // convert eigenvector from Vector to ParGridFunction
+         x.Distribute(eig_solver->GetEigenvector(i));
 
          mode_sock << "parallel " << num_procs << " " << myid << "\n"
                    << "solution\n" << *pmesh << x << flush
@@ -375,7 +435,8 @@ int main(int argc, char *argv[])
    }
 
    // 12. Free the used memory.
-   delete lobpcg;
+   delete eig_solver;
+   delete solver;
    delete precond;
    delete M;
    delete A;

examples/ex11p_arpack.cpp

@@ -0,0 +1,381 @@
+//                       MFEM Example 11 - Parallel Version
+//
+// Compile with: make ex11p
+//
+// Sample runs:  mpirun -np 4 ex11p -m ../data/square-disc.mesh
+//               mpirun -np 4 ex11p -m ../data/star.mesh
+//               mpirun -np 4 ex11p -m ../data/escher.mesh
+//               mpirun -np 4 ex11p -m ../data/fichera.mesh
+//               mpirun -np 4 ex11p -m ../data/square-disc-p2.vtk -o 2
+//               mpirun -np 4 ex11p -m ../data/square-disc-p3.mesh -o 3
+//               mpirun -np 4 ex11p -m ../data/square-disc-nurbs.mesh -o -1
+//               mpirun -np 4 ex11p -m ../data/disc-nurbs.mesh -o -1 -n 20
+//               mpirun -np 4 ex11p -m ../data/pipe-nurbs.mesh -o -1
+//               mpirun -np 4 ex11p -m ../data/ball-nurbs.mesh -o 2
+//               mpirun -np 4 ex11p -m ../data/star-surf.mesh
+//               mpirun -np 4 ex11p -m ../data/square-disc-surf.mesh
+//               mpirun -np 4 ex11p -m ../data/inline-segment.mesh
+//               mpirun -np 4 ex11p -m ../data/amr-quad.mesh
+//               mpirun -np 4 ex11p -m ../data/amr-hex.mesh
+//               mpirun -np 4 ex11p -m ../data/mobius-strip.mesh -n 8
+//               mpirun -np 4 ex11p -m ../data/klein-bottle.mesh -n 10
+//
+// Description:  This example code demonstrates the use of MFEM to solve the
+//               eigenvalue problem -Delta u = lambda u with homogeneous
+//               Dirichlet boundary conditions.
+//
+//               We compute a number of the lowest eigenmodes by discretizing
+//               the Laplacian and Mass operators using a FE space of the
+//               specified order, or an isoparametric/isogeometric space if
+//               order < 1 (quadratic for quadratic curvilinear mesh, NURBS for
+//               NURBS mesh, etc.)
+//
+//               The example highlights the use of the LOBPCG and ARPACK
+//               eigenvalue solvers together with the BoomerAMG preconditioner
+//               in HYPRE, as well as optionally the SuperLU parallel direct
+//               solver. Reusing a single GLVis visualization window for
+//               multiple eigenfunctions is also illustrated.
+//
+//               We recommend viewing Example 1 before viewing this example.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+int main(int argc, char *argv[])
+{
+   // 1. Initialize MPI.
+   int num_procs, myid;
+   MPI_Init(&argc, &argv);
+   MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
+   MPI_Comm_rank(MPI_COMM_WORLD, &myid);
+
+   // 2. Parse command-line options.
+   const char *mesh_file = "../data/star.mesh";
+   int ser_ref_levels = 2;
+   int par_ref_levels = 1;
+   int order = 1;
+   int nev = 5;
+   bool slu_solver = false;
+   bool use_arpack = false;
+   bool visualization = 1;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
+                  "Number of times to refine the mesh uniformly in serial.");
+   args.AddOption(&par_ref_levels, "-rp", "--refine-parallel",
+                  "Number of times to refine the mesh uniformly in parallel.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   args.AddOption(&nev, "-n", "--num-eigs",
+                  "Number of desired eigenmodes.");
+#ifdef MFEM_USE_SUPERLU
+   args.AddOption(&slu_solver, "-slu", "--superlu", "-no-slu",
+                  "--no-superlu", "Use the SuperLU Solver.");
+#endif
+#ifdef MFEM_USE_ARPACK
+   args.AddOption(&use_arpack, "-arpack", "--use-arpack", "-no-arpack",
+                  "--no-arpack",
+                  "Enable or disable the use of ARPACK.");
+#endif
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      MPI_Finalize();
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+
+   // 3. 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;
+   ifstream imesh(mesh_file);
+   if (!imesh)
+   {
+      if (myid == 0)
+      {
+         cerr << "\nCan not open mesh file: " << mesh_file << '\n' << endl;
+      }
+      MPI_Finalize();
+      return 2;
+   }
+   mesh = new Mesh(imesh, 1, 1);
+   imesh.close();
+   int dim = mesh->Dimension();
+
+   // 4. Refine the serial mesh on all processors to increase the resolution. In
+   //    this example we do 'ref_levels' of uniform refinement (2 by default, or
+   //    specified on the command line with -rs).
+   for (int lev = 0; lev < ser_ref_levels; lev++)
+   {
+      mesh->UniformRefinement();
+   }
+
+   // 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
+   //    this mesh further in parallel to increase the resolution (1 time by
+   //    default, or specified on the command line with -rp). Once the parallel
+   //    mesh is defined, the serial mesh can be deleted.
+   ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
+   delete mesh;
+   for (int lev = 0; lev < par_ref_levels; lev++)
+   {
+      pmesh->UniformRefinement();
+   }
+
+   // 6. Define a parallel finite element space on the parallel mesh. Here we
+   //    use continuous Lagrange finite elements of the specified order. If
+   //    order < 1, we instead use an isoparametric/isogeometric space.
+   FiniteElementCollection *fec;
+   if (order > 0)
+   {
+      fec = new H1_FECollection(order, dim);
+   }
+   else if (pmesh->GetNodes())
+   {
+      fec = pmesh->GetNodes()->OwnFEC();
+   }
+   else
+   {
+      fec = new H1_FECollection(order = 1, dim);
+   }
+   ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec);
+   HYPRE_Int size = fespace->GlobalTrueVSize();
+   if (myid == 0)
+   {
+      cout << "Number of unknowns: " << size << endl;
+   }
+
+   // 7. Set up the parallel bilinear forms a(.,.) and m(.,.) on the finite
+   //    element space. The first corresponds to the Laplacian operator -Delta,
+   //    while the second is a simple mass matrix needed on the right hand side
+   //    of the generalized eigenvalue problem below. The boundary conditions
+   //    are implemented by elimination with special values on the diagonal to
+   //    shift the Dirichlet eigenvalues out of the computational range. After
+   //    serial and parallel assembly we extract the corresponding parallel
+   //    matrices A and M.
+   ConstantCoefficient one(1.0);
+   Array<int> ess_bdr;
+   if (pmesh->bdr_attributes.Size())
+   {
+      ess_bdr.SetSize(pmesh->bdr_attributes.Max());
+      ess_bdr = 1;
+   }
+
+   ParBilinearForm *a = new ParBilinearForm(fespace);
+   a->AddDomainIntegrator(new DiffusionIntegrator(one));
+   if (pmesh->bdr_attributes.Size() == 0)
+   {
+      // Add a mass term if the mesh has no boundary, e.g. periodic mesh or
+      // closed surface.
+      a->AddDomainIntegrator(new MassIntegrator(one));
+   }
+   a->Assemble();
+   a->EliminateEssentialBCDiag(ess_bdr, 1.0);
+   a->Finalize();
+
+   ParBilinearForm *m = new ParBilinearForm(fespace);
+   m->AddDomainIntegrator(new MassIntegrator(one));
+   m->Assemble();
+   // shift the eigenvalue corresponding to eliminated dofs to a large value
+   m->EliminateEssentialBCDiag(ess_bdr, numeric_limits<double>::min());
+   m->Finalize();
+
+   HypreParMatrix *A = a->ParallelAssemble();
+   HypreParMatrix *M = m->ParallelAssemble();
+
+#ifdef MFEM_USE_SUPERLU
+   Operator * Arow = NULL;
+   if (slu_solver)
+   {
+      Arow = new SuperLURowLocMatrix(*A);
+   }
+#endif
+
+   delete a;
+   delete m;
+
+   // 8. Define and configure the LOBPCG eigensolver and the BoomerAMG
+   //    preconditioner for A to be used within the solver. Set the matrices
+   //    which define the generalized eigenproblem A x = lambda M x.
+   Eigensolver * esolver = NULL;
+   Solver      *  solver = NULL;
+   Solver      * precond = NULL;
+
+   if (!slu_solver)
+   {
+      HypreBoomerAMG * amg = new HypreBoomerAMG(*A);
+      amg->SetPrintLevel(0);
+      precond = amg;
+
+#ifdef MFEM_USE_ARPACK
+      if ( use_arpack )
+      {
+         HyprePCG * pcg = new HyprePCG(*A);
+         pcg->SetTol(1e-12);
+         pcg->SetMaxIter(200);
+         pcg->SetPreconditioner(*amg);
+         pcg->SetPrintLevel(0);
+
+         solver = pcg;
+      }
+#endif
+   }
+#ifdef MFEM_USE_SUPERLU
+   else
+   {
+      SuperLUSolver * superlu = new SuperLUSolver(MPI_COMM_WORLD);
+      superlu->SetPrintStatistics(false);
+      superlu->SetSymmetricPattern(true);
+      superlu->SetColumnPermutation(superlu::PARMETIS);
+      superlu->SetOperator(*Arow);
+
+      solver  = use_arpack?superlu:NULL;
+      precond = use_arpack?NULL:superlu;
+   }
+#endif
+
+   if ( use_arpack )
+   {
+      ParArPackSym * arpack = new ParArPackSym(MPI_COMM_WORLD);
+      arpack->SetMode(3);
+      arpack->SetPrintLevel(2);
+      arpack->SetSolver(*solver);
+
+      esolver = arpack;
+   }
+   else
+   {
+      HypreLOBPCG * lobpcg = new HypreLOBPCG(MPI_COMM_WORLD);
+      lobpcg->SetPreconditioner(*precond);
+      lobpcg->SetPrecondUsageMode(1);
+      lobpcg->SetPrintLevel(1);
+
+      esolver = lobpcg;
+   }
+
+   esolver->SetNumModes(nev);
+   esolver->SetMaxIter(100);
+   esolver->SetTol(1e-8);
+
+   esolver->SetMassMatrix(*M);
+   esolver->SetOperator(*A);
+
+   // 9. Compute the eigenmodes and extract the array of eigenvalues. Define a
+   //    parallel grid function to represent each of the eigenmodes returned by
+   //    the solver.
+   Array<double> eigenvalues;
+   esolver->Solve();
+   esolver->GetEigenvalues(eigenvalues);
+
+   if ( myid == 0 && use_arpack )
+   {
+      cout << endl;
+      for (int i=0; i<eigenvalues.Size(); i++)
+      {
+         cout << "Eigenvalue lambda " << eigenvalues[i] << endl;
+      }
+      cout << endl;
+   }
+
+   ParGridFunction x(fespace);
+
+   // 10. Save the refined mesh and the modes in parallel. This output can be
+   //     viewed later using GLVis: "glvis -np <np> -m mesh -g mode".
+   {
+      ostringstream mesh_name, mode_name;
+      mesh_name << "mesh." << setfill('0') << setw(6) << myid;
+
+      ofstream mesh_ofs(mesh_name.str().c_str());
+      mesh_ofs.precision(8);
+      pmesh->Print(mesh_ofs);
+
+      for (int i=0; i<nev; i++)
+      {
+         // convert eigenvector from HypreParVector to ParGridFunction
+         x.Distribute(esolver->GetEigenvector(i));
+
+         mode_name << "mode_" << setfill('0') << setw(2) << i << "."
+                   << setfill('0') << setw(6) << myid;
+
+         ofstream mode_ofs(mode_name.str().c_str());
+         mode_ofs.precision(8);
+         x.Save(mode_ofs);
+         mode_name.str("");
+      }
+   }
+
+   // 11. Send the solution by socket to a GLVis server.
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream mode_sock(vishost, visport);
+      mode_sock.precision(8);
+
+      for (int i=0; i<nev; i++)
+      {
+         if ( myid == 0 )
+         {
+            cout << "Eigenmode " << i+1 << '/' << nev
+                 << ", Lambda = " << eigenvalues[i] << endl;
+         }
+
+         // convert eigenvector from HypreParVector to ParGridFunction
+         x.Distribute(esolver->GetEigenvector(i));
+
+         mode_sock << "parallel " << num_procs << " " << myid << "\n"
+                   << "solution\n" << *pmesh << x << flush
+                   << "window_title 'Eigenmode " << i+1 << '/' << nev
+                   << ", Lambda = " << eigenvalues[i] << "'" << endl;
+
+         char c;
+         if (myid == 0)
+         {
+            cout << "press (q)uit or (c)ontinue --> " << flush;
+            cin >> c;
+         }
+         MPI_Bcast(&c, 1, MPI_CHAR, 0, MPI_COMM_WORLD);
+
+         if (c != 'c')
+         {
+            break;
+         }
+      }
+      mode_sock.close();
+   }
+
+   // 12. Free the used memory.
+   delete esolver;
+   delete solver;
+   delete precond;
+   delete M;
+   delete A;
+
+   delete fespace;
+   if (order > 0)
+   {
+      delete fec;
+   }
+   delete pmesh;
+
+   MPI_Finalize();
+
+   return 0;
+}

examples/ex12p.cpp

@@ -276,8 +276,8 @@ int main(int argc, char *argv[])
 
       for (int i=0; i<nev; i++)
       {
-         // convert eigenvector from HypreParVector to ParGridFunction
-         x = lobpcg->GetEigenvector(i);
+         // convert eigenvector from Vector to ParGridFunction
+         x.Distribute(lobpcg->GetEigenvector(i));
 
          mode_name << "mode_" << setfill('0') << setw(2) << i << "."
                    << setfill('0') << setw(6) << myid;
@@ -303,7 +303,7 @@ int main(int argc, char *argv[])
       pmesh->Print(adios2output);
       for (int i=0; i<nev; i++)
       {
-         x = lobpcg->GetEigenvector(i);
+         x.Distribute(lobpcg->GetEigenvector(i));
          // x is a temporary that must be saved immediately
          x.Save(adios2output, "mode_" + std::to_string(i));
       }
@@ -326,8 +326,8 @@ int main(int argc, char *argv[])
                  << ", Lambda = " << eigenvalues[i] << endl;
          }
 
-         // convert eigenvector from HypreParVector to ParGridFunction
-         x = lobpcg->GetEigenvector(i);
+         // convert eigenvector from Vector to ParGridFunction
+         x.Distribute(lobpcg->GetEigenvector(i));
 
          mode_sock << "parallel " << num_procs << " " << myid << "\n"
                    << "solution\n" << *pmesh << x << flush

examples/ex13.cpp

@@ -0,0 +1,282 @@
+//                       MFEM Example 13
+//
+// Compile with: make ex3p
+//
+// Sample runs:  ex13 -m ../data/star.mesh -s 5
+//               ex13 -m ../data/square-disc.mesh -o 2 -n 4 // minres fails to conv.
+//               ex13 -m ../data/beam-hex.mesh
+//               ex13 -m ../data/square-disc.mesh -rs 1 -s 26
+//               ex13 -m ../data/square-disc-nurbs.mesh -rs 3 -s 26
+//               ex13 -m ../data/amr-quad.mesh -o 2 // minres fails to conv.
+//               ex13 -m ../data/mobius-strip.mesh -n 8
+//
+// Description:  This example code solves a simple 3D electromagnetic
+//               eigenmode problem corresponding to the second order
+//               Maxwell equation curl curl E = lambda E with boundary
+//               condition E x n = 0. We discretize with Nedelec finite
+//               elements.
+//
+//               The example demonstrates the use of H(curl) finite element
+//               spaces with the curl-curl and the (vector finite element) mass
+//               bilinear form, as well as the use of the ARPACK eigenmode
+//               solver for symmetric matrices using the shift-invert mode.
+//
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+#ifdef MFEM_USE_ARPACK
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   const char *mesh_file = "../data/beam-tet.mesh";
+   int order = 1;
+   int nev = 5;
+   int sr = 2;
+   double sigma = 11.0;
+   bool visualization = 1;
+   bool arp_solver = true;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree).");
+   args.AddOption(&nev, "-n", "--num-eigs",
+                  "Number of desired eigenmodes.");
+   args.AddOption(&sr, "-rs", "--refine-serial",
+                  "Number of times to refine the mesh uniformly in serial.");
+   args.AddOption(&sigma, "-s", "--shift",
+                  "Average of the desired eigenvalue range.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   // 2. 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;
+   ifstream imesh(mesh_file);
+   if (!imesh)
+   {
+      cerr << "\nCan not open mesh file: " << mesh_file << '\n' << endl;
+      return 2;
+   }
+   mesh = new Mesh(imesh, 1, 1);
+   imesh.close();
+   int dim = mesh->Dimension();
+
+   // 3. Refine the mesh to increase the resolution. In this example we do
+   //    'ref_levels' of uniform refinement.
+   {
+      int ref_levels = sr;
+      for (int l = 0; l < ref_levels; l++)
+      {
+         mesh->UniformRefinement();
+      }
+   }
+
+   // 4. Define a finite element space on the mesh. Here we use the lowest
+   //    order Nedelec finite elements, but we can easily switch
+   //    to higher-order spaces by changing the value of p.
+   FiniteElementCollection *fec = new ND_FECollection(order, dim);
+   FiniteElementSpace *fespace = new FiniteElementSpace(mesh, fec);
+   int size = fespace->GetVSize();
+
+   cout << "Number of unknowns: " << size << endl;
+   cout << "Number of boundary attributes: " << mesh->bdr_attributes.Max()
+        << endl;
+
+   // 5. Set up the parallel bilinear form corresponding to the EM diffusion
+   //    operator curl muinv curl - sigma I, by adding the curl-curl and the
+   //    mass domain integrators and finally imposing homogeneous Dirichlet
+   //    boundary conditions. The boundary conditions are implemented by
+   //    marking all the boundary attributes from the mesh as essential
+   //    (Dirichlet). After serial and parallel assembly we extract the
+   //    parallel matrices A and M.
+   Coefficient *muinv    = new ConstantCoefficient(1.0);
+   Coefficient *negSigma = new ConstantCoefficient(-sigma);
+
+   BilinearForm *a = new BilinearForm(fespace);
+   a->AddDomainIntegrator(new CurlCurlIntegrator(*muinv));
+   a->AddDomainIntegrator(new VectorFEMassIntegrator(*negSigma));
+   a->Assemble();
+   Array<int> ess_bdr(mesh->bdr_attributes.Max());
+   ess_bdr = 1;
+   a->EliminateEssentialBC(ess_bdr);
+   a->Finalize();
+
+   BilinearForm *m = new BilinearForm(fespace);
+   m->AddDomainIntegrator(new VectorFEMassIntegrator());
+   m->Assemble();
+   m->EliminateEssentialBCDiag(ess_bdr, sqrt(numeric_limits<double>::min()));
+   m->Finalize();
+
+   // 6. Define a parallel grid function to approximate each of the
+   // eigenmodes returned by the solver.  Use this as a template to
+   // create a special multi-vector object needed by the eigensolver
+   // which is then initialized with random values.
+   GridFunction x(fespace);
+   x = 0.0;
+
+   // 7. Define and configure the GMRES
+   //    solver to be used within the eigensolver.
+   Solver    * solver = NULL;
+   if ( false )
+   {
+      GMRESSolver *  gmres = new GMRESSolver();
+
+      gmres->SetOperator(*a);
+      gmres->SetRelTol(1e-8);
+      gmres->SetMaxIter(1000);
+      gmres->SetPrintLevel(0);
+      solver = gmres;
+   }
+   else
+   {
+#ifndef MFEM_USE_SUITESPARSE
+      cout << "Building MINRESSolver" << endl;
+      MINRESSolver * minres = new MINRESSolver();
+
+      minres->SetRelTol(1e-12);
+      minres->SetMaxIter(1000);
+      minres->SetPrintLevel(0);
+      solver = minres;
+#else
+      // 7. If MFEM was compiled with SuiteSparse, use UMFPACK to solve the system.
+      cout << "Building UMFPackSolver" << endl;
+      UMFPackSolver * umf_solver = new UMFPackSolver;
+      umf_solver->Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
+      solver = umf_solver;
+#endif
+   }
+   solver->SetOperator(a->SpMat());
+
+   // 7. Define and configure the ARPACK eigensolver
+   SymGenEigensolver * eig_solver = NULL;
+   if (arp_solver)
+   {
+      ArPackSAUPD * arpack = new ArPackSAUPD();
+
+      arpack->SetNumModes(nev);
+      arpack->SetMaxIter(400);
+      arpack->SetTol(1e-8);
+      arpack->SetShift(sigma);
+      arpack->SetMode(3);
+      arpack->SetPrintLevel(2);
+      arpack->SetSolver(*solver);
+
+      eig_solver = arpack;
+   }
+
+   eig_solver->SetOperators(*a, *m);
+
+   // Obtain the eigenvalues and eigenvectors
+   Array<double> eigenvalues(nev);
+   eigenvalues = -1.0;
+
+   // arpack->Solve(eigenvalues, *eigenvectors);
+   eig_solver->Solve();
+
+   eig_solver->GetEigenvalues(eigenvalues);
+
+   cout << endl;
+   std::ios::fmtflags old_fmt = cout.flags();
+   cout.setf(std::ios::scientific);
+   std::streamsize old_prec = cout.precision(14);
+   for (int i=0; i<min(nev,eigenvalues.Size()); i++)
+   {
+      cout << "Eigenvalue lambda   " << eigenvalues[i] << endl;
+   }
+   cout.precision(old_prec);
+   cout.flags(old_fmt);
+   cout << endl;
+
+   VisItDataCollection visit_dc("Example13", mesh);
+   GridFunction ** mode = new GridFunction*[min(nev,eigenvalues.Size())];
+   for (int i=0; i<min(nev,eigenvalues.Size()); i++)
+   {
+      mode[i] = new GridFunction(fespace);
+      *mode[i] = eig_solver->GetEigenvector(i);
+
+      ostringstream modeName;
+      modeName << "mode_" << setfill('0') << setw(2) << i;
+      visit_dc.RegisterField(modeName.str().c_str(),mode[i]);
+   }
+   visit_dc.Save();
+
+   // 8. Save the refined mesh and the modes. This output can
+   //    be viewed later using GLVis: "glvis -m mesh -g mode".
+   {
+      ofstream mesh_ofs("refined.mesh");
+      mesh_ofs.precision(8);
+      mesh->Print(mesh_ofs);
+
+      for (int i=0; i<min(nev,eigenvalues.Size()); i++)
+      {
+         x = eig_solver->GetEigenvector(i);
+
+         ostringstream modeName;
+         modeName << "mode_" << setfill('0') << setw(2) << i;
+
+         ofstream mode_ofs(modeName.str().c_str());
+         mode_ofs.precision(8);
+         x.Save(mode_ofs);
+         modeName.str("");
+      }
+   }
+
+   // 9. Send the solution by socket to a GLVis server.
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream mode_sock(vishost, visport);
+      mode_sock.precision(8);
+
+      for (int i=0; i<min(nev,eigenvalues.Size()); i++)
+      {
+         x = eig_solver->GetEigenvector(i);
+
+         mode_sock << "solution\n" << *mesh << x << flush;
+
+         char c;
+         cout << "press (q)uit or (c)ontinue --> " << flush;
+         cin >> c;
+
+         if (c != 'c')
+         {
+            break;
+         }
+      }
+      mode_sock.close();
+   }
+
+   // 10. Free the used memory.
+   delete a;
+   delete m;
+   delete negSigma;
+   delete muinv;
+   delete eig_solver;
+   delete solver;
+   // delete X;
+   delete fespace;
+   delete fec;
+   delete mesh;
+
+   return 0;
+}
+
+#endif // MFEM_USE_ARPACK

examples/ex13p.cpp

@@ -215,8 +215,8 @@ int main(int argc, char *argv[])
 
       for (int i=0; i<nev; i++)
       {
-         // convert eigenvector from HypreParVector to ParGridFunction
-         x = ame->GetEigenvector(i);
+         // convert eigenvector from Vector to ParGridFunction
+         x.Distribute(ame->GetEigenvector(i));
 
          mode_name << "mode_" << setfill('0') << setw(2) << i << "."
                    << setfill('0') << setw(6) << myid;
@@ -244,8 +244,8 @@ int main(int argc, char *argv[])
                  << ", Lambda = " << eigenvalues[i] << endl;
          }
 
-         // convert eigenvector from HypreParVector to ParGridFunction
-         x = ame->GetEigenvector(i);
+         // convert eigenvector from Vector to ParGridFunction
+         x.Distribute(ame->GetEigenvector(i));
 
          mode_sock << "parallel " << num_procs << " " << myid << "\n"
                    << "solution\n" << *pmesh << x << flush

examples/ex32p.cpp

@@ -228,7 +228,7 @@ int main(int argc, char *argv[])
       for (int i=0; i<nev; i++)
       {
          // convert eigenvector from HypreParVector to ParGridFunction
-         x = ame->GetEigenvector(i);
+         x.Distribute(ame->GetEigenvector(i));
          curl.Mult(x, dx);
 
          mode_name << "mode_" << setfill('0') << setw(2) << i << "."
@@ -295,7 +295,7 @@ int main(int argc, char *argv[])
             }
 
             // convert eigenvector from HypreParVector to ParGridFunction
-            x = ame->GetEigenvector(i);
+            x.Distribute(ame->GetEigenvector(i));
             curl.Mult(x, dx);
 
             {
@@ -469,7 +469,7 @@ int main(int argc, char *argv[])
             }
 
             // convert eigenvector from HypreParVector to ParGridFunction
-            x = ame->GetEigenvector(i);
+            x.Distribute(ame->GetEigenvector(i));
             curl.Mult(x, dx);
 
             {
@@ -599,7 +599,7 @@ int main(int argc, char *argv[])
             }
 
             // convert eigenvector from HypreParVector to ParGridFunction
-            x = ame->GetEigenvector(i);
+            x.Distribute(ame->GetEigenvector(i));
             curl.Mult(x, dx);
 
             mode_sock << "parallel " << num_procs << " " << myid << "\n"

examples/ex35p.cpp

@@ -658,7 +658,7 @@ void ScalarWaveGuide(int mode, ParGridFunction &x)
    lobpcg.SetOperator(*A);
    lobpcg.Solve();
 
-   x = lobpcg.GetEigenvector(mode);
+   x.Distribute(lobpcg.GetEigenvector(mode));
 
    delete A;
    delete M;
@@ -714,7 +714,7 @@ void VectorWaveGuide(int mode, ParGridFunction &x)
    ame.SetOperator(*A);
    ame.Solve();
 
-   x = ame.GetEigenvector(mode);
+   x.Distribute(ame.GetEigenvector(mode));
 
    delete A;
    delete M;
@@ -780,7 +780,7 @@ void PseudoScalarWaveGuide(int mode, ParGridFunction &x_l2)
    lobpcg.SetOperator(*A);
    lobpcg.Solve();
 
-   x = lobpcg.GetEigenvector(mode);
+   x.Distribute(lobpcg.GetEigenvector(mode));
 
    x_l2.ProjectCoefficient(xCoef);
 

examples/ginkgo/makefile

@@ -76,4 +76,4 @@ clean-build:
 	rm -rf *.dSYM *.TVD.*breakpoints
 
 clean-exec:
-	@rm -f refined.mesh sol.gf mesh.* sol.*
+	@rm -f refined.mesh sol.gf

examples/makefile

@@ -31,6 +31,9 @@ SEQ_DEVICE_EXAMPLES = ex1 ex3 ex4 ex5 ex6 ex9 ex14 ex22 ex24 ex25 ex26 ex34
 PAR_DEVICE_EXAMPLES = ex1p ex2p ex3p ex4p ex5p ex6p ex7p ex9p ex13p ex14p \
    ex22p ex24p ex25p ex26p ex34p ex35p
 
+ifeq ($(MFEM_USE_ARPACK),YES)
+   SEQ_EXAMPLES += ex11 ex13
+endif
 ifeq ($(MFEM_USE_LAPACK),YES)
    SEQ_EXAMPLES += ex38
 endif
@@ -71,7 +74,6 @@ endif
 
 SUBDIRS_ALL = $(addsuffix /all,$(SUBDIRS))
 SUBDIRS_TEST = $(addsuffix /test,$(SUBDIRS))
-SUBDIRS_TEST_NOCLEAN = $(addsuffix /test-noclean,$(SUBDIRS))
 SUBDIRS_CLEAN = $(addsuffix /clean,$(SUBDIRS))
 SUBDIRS_TPRINT = $(addsuffix /test-print,$(SUBDIRS))
 
@@ -88,9 +90,8 @@ SUBDIRS_TPRINT = $(addsuffix /test-print,$(SUBDIRS))
 
 all: $(EXAMPLES) $(SUBDIRS_ALL)
 
-.PHONY: $(SUBDIRS_ALL) $(SUBDIRS_TEST) $(SUBDIRS_TEST_NOCLEAN) \
-   $(SUBDIRS_CLEAN) $(SUBDIRS_TPRINT)
-$(SUBDIRS_ALL) $(SUBDIRS_TEST) $(SUBDIRS_TEST_NOCLEAN) $(SUBDIRS_CLEAN):
+.PHONY: $(SUBDIRS_ALL) $(SUBDIRS_TEST) $(SUBDIRS_CLEAN) $(SUBDIRS_TPRINT)
+$(SUBDIRS_ALL) $(SUBDIRS_TEST) $(SUBDIRS_CLEAN):
 	$(MAKE) -C $(@D) $(@F)
 $(SUBDIRS_TPRINT):
 	@$(MAKE) -C $(@D) $(@F)
@@ -109,7 +110,6 @@ endif
 MFEM_TESTS = EXAMPLES
 include $(MFEM_TEST_MK)
 test: $(SUBDIRS_TEST)
-test-noclean: $(SUBDIRS_TEST_NOCLEAN)
 test-print: $(SUBDIRS_TPRINT)
 
 # Testing: Parallel vs. serial runs

fem/bounds.cpp

@@ -41,14 +41,9 @@ void PLBound::Setup(const int nb_i, const int ncp_i,
    tol = tol_i;
    lbound.SetSize(ncp, nb);
    ubound.SetSize(ncp, nb);
-   lbound_t.SetSize(nb, ncp);
-   ubound_t.SetSize(nb, ncp);
    nodes.SetSize(nb);
    weights.SetSize(nb);
    control_points.SetSize(ncp);
-   xhat.SetSize(nb);
-   what.SetSize(nb);
-   cphat.SetSize(ncp);
 
    auto scalenodes = [](const Vector &in, const real_t a, const real_t b) -> Vector
    {
@@ -95,10 +90,6 @@ void PLBound::Setup(const int nb_i, const int ncp_i,
       MFEM_ABORT("Unsupported interval points. Use [0,1].\n");
    }
    control_points = scalenodes(control_points, 0.0, 1.0); // rescale to [0,1]
-   for (int i = 0; i < ncp; i++)
-   {
-      cphat(i) = 2.0*control_points(i) - 1.0;
-   }
 
    Poly_1D::Basis &basis1d(poly1d.GetBasis(nb-1, b_type));
 
@@ -154,8 +145,6 @@ void PLBound::Setup(const int nb_i, const int ncp_i,
                lbound(j,i) = std::max(lbound(j,i),0_r);
             }
          }
-         lbound_t(i,j) = lbound(j,i);
-         ubound_t(i,j) = ubound(j,i);
       }
    }
 
@@ -187,11 +176,6 @@ void PLBound::Setup(const int nb_i, const int ncp_i,
          nodes(i) = irule.IntPoint(i).x;
       }
    }
-   for (int i = 0; i < nb; i++)
-   {
-      xhat(i) = 2.0*nodes(i) - 1.0;
-      what(i) = 2.0*weights(i);
-   }
 
    if (b_type == 2)
    {
@@ -771,4 +755,4 @@ void PLBound::Print(std::ostream &outp) const
    ubound.Print(outp);
 }
 
-}
+}

fem/bounds.hpp

@@ -13,7 +13,6 @@
 #define MFEM_BOUNDS
 
 #include "../config/config.hpp"
-#include "../general/forall.hpp"
 #include "fespace.hpp"
 
 namespace mfem
@@ -61,9 +60,7 @@ private:
    bool proj = true; // Use linear projection to compute bounds.
    real_t tol = 0.0; // offset bounds to avoid round-off errors
    Vector nodes, weights, control_points;
-   Vector xhat, what, cphat;
    DenseMatrix lbound, ubound; // ncp x nb matrices with bounds of all bases
-   DenseMatrix lbound_t, ubound_t; // nb x ncp transposes for device kernel
    // Some auxillary storage for computing the bounds with Bernstein
    DenseMatrix basisMatNodes; // Bernstein bases at equispaced nodes
    DenseMatrix basisMatInt;   // Bernstein bases at GLL nodes
@@ -116,10 +113,7 @@ public:
     *  @details This projection increases the computational cost but results in
     *  tighter bounds.
     */
-   void SetProjectionFlagForBounding(bool proj_)
-   {
-      proj = proj_;
-   }
+   void SetProjectionFlagForBounding(bool proj_) { proj = proj_; }
 
    /** @brief Compute piecewise linear bounds for the lexicographically-ordered
     *  nodal coefficients in @a coeff in 1D/2D/3D.
@@ -143,23 +137,9 @@ public:
    /// Get number of control points used to compute the bounds.
    int GetNControlPoints() const { return ncp; }
 
-   /// Get the underlying 1D basis type.
-   int GetBasisType() const { return b_type; }
-
    /// Get 1D control point locations (lexicographic order) in [0,1].
    const Vector &GetControlPoints() const { return control_points; }
 
-   /** @brief Compute element-wise bounds from a lexicographic E-vector.
-    *
-    *  @details The expected layout of @a e_vec is `ND x VDIM x NE`, where
-    *  `ND = nb^rdim`, `VDIM = fes_vdim`, and `NE` is the number of elements.
-    *  The output layout matches GridFunction::GetElementBounds:
-    *  `NE x active_vdim`, with the element index varying fastest.
-    */
-   void GetElementBoundsKernel(const int rdim, const int fes_vdim,
-                               const Vector &e_vec, Vector &lower,
-                               Vector &upper, const int vdim = 0) const;
-
    /** @brief Get lower and upper bounding matrix (ncp^dim x nb^dim)
     *
     *  @details The matrices can be used to compute the bounds at control points
@@ -203,600 +183,6 @@ private:
               const int cp_type_i, const real_t tol_i);
 };
 
-namespace internal
-{
-
-struct PLBoundDeviceData
-{
-   int nb;
-   int ncp;
-   const real_t *xhat;
-   const real_t *what;
-   const real_t *cphat;
-   const real_t *lbound;
-   const real_t *ubound;
-};
-
-template<int T_NB = 0, bool T_PROJ = true>
-inline void GetElementBoundsKernel1D(const PLBoundDeviceData &data,
-                                     const int fes_vdim,
-                                     const int ne,
-                                     const Vector &e_vec,
-                                     Vector &lower,
-                                     Vector &upper,
-                                     const int comp0,
-                                     const int ncomp)
-{
-   constexpr int GENERIC_MAX_ND = 32;
-   constexpr int MAX_ND = T_NB ? T_NB : GENERIC_MAX_ND;
-   constexpr int BLOCK_X = 2*MAX_ND;
-
-   const int nd = T_NB ? T_NB : data.nb;
-   MFEM_VERIFY(nd <= MAX_ND,
-               "Device element bounds kernel supports up to 32 "
-               "1D degrees of freedom.");
-
-   const auto E = Reshape(e_vec.Read(), nd, fes_vdim, ne);
-   auto L = Reshape(lower.Write(), ne, ncomp);
-   auto U = Reshape(upper.Write(), ne, ncomp);
-
-   mfem::forall_2D<BLOCK_X>(ne*ncomp, BLOCK_X, 1,
-                            [=] MFEM_HOST_DEVICE (int ec)
-   {
-      const int e = ec % ne;
-      const int c = ec / ne;
-      const int vc = comp0 + c;
-      const real_t *coeff = &E(0, vc, e);
-      const int tid = MFEM_THREAD_ID(x);
-
-      MFEM_SHARED real_t sproj[MAX_ND];
-      MFEM_SHARED real_t ssum0[MAX_ND];
-      MFEM_SHARED real_t ssum1[MAX_ND];
-      MFEM_SHARED real_t smin[BLOCK_X];
-      MFEM_SHARED real_t smax[BLOCK_X];
-      MFEM_SHARED real_t sa0;
-      MFEM_SHARED real_t sa1;
-
-      MFEM_FOREACH_THREAD(i, x, nd)
-      {
-         if constexpr (T_PROJ)
-         {
-            const real_t x = data.xhat[i];
-            const real_t w = data.what[i];
-            ssum0[i] = 0.5*coeff[i]*w;
-            ssum1[i] = 1.5*coeff[i]*w*x;
-         }
-         else
-         {
-            ssum0[i] = 0.0;
-            ssum1[i] = 0.0;
-         }
-      }
-      MFEM_SYNC_THREAD;
-
-      MFEM_FOREACH_THREAD(ii, x, 1)
-      {
-         sa0 = 0.0;
-         sa1 = 0.0;
-         for (int i = 0; i < nd; i++)
-         {
-            sa0 += ssum0[i];
-            sa1 += ssum1[i];
-         }
-      }
-      MFEM_SYNC_THREAD;
-
-      MFEM_FOREACH_THREAD(i, x, nd)
-      {
-         if constexpr (T_PROJ)
-         {
-            const real_t x = data.xhat[i];
-            sproj[i] = coeff[i] - sa0 - sa1*x;
-         }
-         else
-         {
-            sproj[i] = coeff[i];
-         }
-      }
-      MFEM_SYNC_THREAD;
-
-      real_t lower_local = HUGE_VAL;
-      real_t upper_local = -HUGE_VAL;
-      MFEM_FOREACH_THREAD(j, x, data.ncp)
-      {
-         real_t lo = 0.0;
-         real_t hi = 0.0;
-         if constexpr (T_PROJ)
-         {
-            const real_t xcp = data.cphat[j];
-            lo = sa0 + sa1*xcp;
-            hi = lo;
-         }
-
-         for (int i = 0; i < nd; i++)
-         {
-            const real_t val = sproj[i];
-            const real_t lv = data.lbound[j + i*data.ncp]*val;
-            const real_t uv = data.ubound[j + i*data.ncp]*val;
-            lo += lv < uv ? lv : uv;
-            hi += lv > uv ? lv : uv;
-         }
-         lower_local = lower_local < lo ? lower_local : lo;
-         upper_local = upper_local > hi ? upper_local : hi;
-      }
-
-      smin[tid] = lower_local;
-      smax[tid] = upper_local;
-      MFEM_SYNC_THREAD;
-
-      MFEM_FOREACH_THREAD(ii, x, 1)
-      {
-         real_t lower_ec = smin[0];
-         real_t upper_ec = smax[0];
-         const int nthreads = MFEM_THREAD_SIZE(x);
-         const int nactive = data.ncp < nthreads ? data.ncp : nthreads;
-         for (int t = 1; t < nactive; t++)
-         {
-            lower_ec = lower_ec < smin[t] ? lower_ec : smin[t];
-            upper_ec = upper_ec > smax[t] ? upper_ec : smax[t];
-         }
-         L(e, c) = lower_ec;
-         U(e, c) = upper_ec;
-      }
-   });
-}
-
-template<int T_NB = 0, int T_NCP = 0, bool T_PROJ = true>
-inline void GetElementBoundsKernel2D(const PLBoundDeviceData &data,
-                                     const int fes_vdim,
-                                     const int ne,
-                                     const Vector &e_vec,
-                                     Vector &lower,
-                                     Vector &upper,
-                                     const int comp0,
-                                     const int ncomp)
-{
-   constexpr int DEFAULT_MAX_NB = 8;
-   constexpr int DEFAULT_MAX_CP = 3*DEFAULT_MAX_NB;
-   constexpr int MAX_NB = T_NB ? T_NB : DEFAULT_MAX_NB;
-   constexpr int MAX_CP = T_NCP ? T_NCP : DEFAULT_MAX_CP;
-   constexpr int MAX_THREADS = MAX_CP*MAX_CP;
-
-   const int nb = data.nb;
-   const int ncp = data.ncp;
-   const int nd = nb*nb;
-   MFEM_VERIFY(nb <= MAX_NB,
-               "Device 2D element bounds kernel exceeds its compile-time "
-               "1D degree bound.");
-   MFEM_VERIFY(ncp <= MAX_CP,
-               "Device 2D element bounds kernel exceeds its compile-time "
-               "control-point bound.");
-   MFEM_VERIFY(ncp*ncp <= MAX_THREADS,
-               "Device 2D element bounds kernel exceeds its compile-time "
-               "thread-block bound.");
-
-   const auto E = Reshape(e_vec.Read(), nd, fes_vdim, ne);
-   auto L = Reshape(lower.Write(), ne, ncomp);
-   auto U = Reshape(upper.Write(), ne, ncomp);
-
-   mfem::forall_2D<MAX_THREADS>(ne*ncomp, ncp, ncp,
-                                [=] MFEM_HOST_DEVICE (int ec)
-   {
-      const int e = ec % ne;
-      const int c = ec / ne;
-      const int vc = comp0 + c;
-      const real_t *coeff = &E(0, vc, e);
-      const int tx = MFEM_THREAD_ID(x);
-      const int ty = MFEM_THREAD_ID(y);
-
-      MFEM_SHARED real_t sproj[MAX_NB*MAX_NB];
-      MFEM_SHARED real_t srow_min[MAX_NB*MAX_CP];
-      MFEM_SHARED real_t srow_max[MAX_NB*MAX_CP];
-      MFEM_SHARED real_t srow_a0[MAX_NB];
-      MFEM_SHARED real_t srow_a1[MAX_NB];
-      MFEM_SHARED real_t sa0[MAX_CP];
-      MFEM_SHARED real_t sa1[MAX_CP];
-      MFEM_SHARED real_t smin[MAX_THREADS];
-      MFEM_SHARED real_t smax[MAX_THREADS];
-
-      // Stage 1a: for each nodal row, form the per-node contributions to the
-      // row-wise linear fit used by the first 1D bounding solve.
-      MFEM_FOREACH_THREAD(jrow, y, nb)
-      {
-         const real_t *row_coeff = coeff + jrow*nb;
-         const int row_ncp_off = jrow*MAX_CP;
-         MFEM_FOREACH_THREAD(i, x, nb)
-         {
-            if constexpr (T_PROJ)
-            {
-               const real_t x = data.xhat[i];
-               const real_t w = data.what[i];
-               srow_min[row_ncp_off + i] = 0.5*row_coeff[i]*w;
-               srow_max[row_ncp_off + i] = 1.5*row_coeff[i]*w*x;
-            }
-            else
-            {
-               srow_min[row_ncp_off + i] = 0.0;
-               srow_max[row_ncp_off + i] = 0.0;
-            }
-         }
-      }
-      MFEM_SYNC_THREAD;
-
-      // Stage 1b: reduce the row-wise projection coefficients a0/a1.
-      if constexpr (T_PROJ)
-      {
-         MFEM_FOREACH_THREAD(jrow, y, nb)
-         {
-            const int row_ncp_off = jrow*MAX_CP;
-            real_t a0 = 0.0;
-            real_t a1 = 0.0;
-            MFEM_FOREACH_THREAD(ii, x, 1)
-            {
-               for (int i = 0; i < nb; i++)
-               {
-                  a0 += srow_min[row_ncp_off + i];
-                  a1 += srow_max[row_ncp_off + i];
-               }
-               srow_a0[jrow] = a0;
-               srow_a1[jrow] = a1;
-            }
-         }
-         MFEM_SYNC_THREAD;
-      }
-
-      // Stage 1c: subtract the row-wise linear fit once and cache the
-      // projected row coefficients for reuse across all x-control points.
-      MFEM_FOREACH_THREAD(jrow, y, nb)
-      {
-         const real_t *row_coeff = coeff + jrow*nb;
-         MFEM_FOREACH_THREAD(i, x, nb)
-         {
-            if constexpr (T_PROJ)
-            {
-               const real_t x = data.xhat[i];
-               sproj[jrow*MAX_NB + i] = row_coeff[i]
-                                        - srow_a0[jrow] - srow_a1[jrow]*x;
-            }
-            else
-            {
-               sproj[jrow*MAX_NB + i] = row_coeff[i];
-            }
-         }
-      }
-      MFEM_SYNC_THREAD;
-
-      // Stage 1d: solve the first 1D bounding problem along each nodal row and
-      // store bounds at every x-direction control point.
-      MFEM_FOREACH_THREAD(icp, x, ncp)
-      {
-         MFEM_FOREACH_THREAD(jrow, y, nb)
-         {
-            const int row_cp_off = jrow*ncp;
-            real_t lo = 0.0;
-            real_t hi = 0.0;
-            if constexpr (T_PROJ)
-            {
-               const real_t xcp = data.cphat[icp];
-               lo = srow_a0[jrow] + srow_a1[jrow]*xcp;
-               hi = lo;
-            }
-            for (int i = 0; i < nb; i++)
-            {
-               const real_t val = sproj[jrow*MAX_NB + i];
-               const real_t lv = data.lbound[icp + i*data.ncp]*val;
-               const real_t uv = data.ubound[icp + i*data.ncp]*val;
-               lo += lv < uv ? lv : uv;
-               hi += lv > uv ? lv : uv;
-            }
-            srow_min[row_cp_off + icp] = lo;
-            srow_max[row_cp_off + icp] = hi;
-         }
-      }
-      MFEM_SYNC_THREAD;
-
-      // Stage 2a: from the row bounds, form the per-row contributions to the
-      // second 1D projection solve in the y-direction.
-      MFEM_FOREACH_THREAD(icp, x, ncp)
-      {
-         MFEM_FOREACH_THREAD(jrow, y, nb)
-         {
-            const int row_cp_off = jrow*ncp;
-            if constexpr (T_PROJ)
-            {
-               const real_t x = data.xhat[jrow];
-               const real_t w = data.what[jrow];
-               const real_t t = 0.5*(srow_min[row_cp_off + icp] +
-                                     srow_max[row_cp_off + icp]);
-               smin[row_cp_off + icp] = 0.5*t*w;
-               smax[row_cp_off + icp] = 1.5*t*w*x;
-            }
-            else
-            {
-               smin[row_cp_off + icp] = 0.0;
-               smax[row_cp_off + icp] = 0.0;
-            }
-         }
-      }
-      MFEM_SYNC_THREAD;
-
-      // Stage 2b: reduce the y-direction projection coefficients for each
-      // x-control-point column.
-      MFEM_FOREACH_THREAD(jj, y, 1)
-      {
-         MFEM_FOREACH_THREAD(icp, x, ncp)
-         {
-            real_t a0 = 0.0;
-            real_t a1 = 0.0;
-            for (int jrow = 0; jrow < nb; jrow++)
-            {
-               a0 += smin[jrow*ncp + icp];
-               a1 += smax[jrow*ncp + icp];
-            }
-            sa0[icp] = a0;
-            sa1[icp] = a1;
-         }
-      }
-      MFEM_SYNC_THREAD;
-
-      // Stage 2c: subtract the y-direction linear fit from the intermediate
-      // row bounds so the final tensor-product bound uses the perturbation.
-      if constexpr (T_PROJ)
-      {
-         MFEM_FOREACH_THREAD(icp, x, ncp)
-         {
-            MFEM_FOREACH_THREAD(jrow, y, nb)
-            {
-               const int row_cp_off = jrow*ncp;
-               const real_t x = data.xhat[jrow];
-               const real_t t = sa0[icp] + sa1[icp]*x;
-               srow_min[row_cp_off + icp] -= t;
-               srow_max[row_cp_off + icp] -= t;
-            }
-         }
-      }
-      MFEM_SYNC_THREAD;
-
-      // Stage 3: each thread now owns one 2D control point (icp, kcp) and
-      // accumulates its final lower/upper bound from the row-bound data.
-      MFEM_FOREACH_THREAD(icp, x, ncp)
-      {
-         MFEM_FOREACH_THREAD(kcp, y, ncp)
-         {
-            real_t lo = 0.0;
-            real_t hi = 0.0;
-            if constexpr (T_PROJ)
-            {
-               const real_t xcp = data.cphat[kcp];
-               lo = sa0[icp] + sa1[icp]*xcp;
-               hi = lo;
-            }
-            for (int jrow = 0; jrow < nb; jrow++)
-            {
-               const real_t w0 = srow_min[jrow*ncp + icp];
-               const real_t w1 = srow_max[jrow*ncp + icp];
-               const real_t lb = data.lbound[kcp + jrow*data.ncp];
-               const real_t ub = data.ubound[kcp + jrow*data.ncp];
-               const real_t v0 = lb*w0;
-               const real_t v1 = ub*w0;
-               const real_t v2 = lb*w1;
-               const real_t v3 = ub*w1;
-               real_t vlo = v0 < v1 ? v0 : v1;
-               real_t vhi = v0 > v1 ? v0 : v1;
-               vlo = vlo < v2 ? vlo : v2;
-               vlo = vlo < v3 ? vlo : v3;
-               vhi = vhi > v2 ? vhi : v2;
-               vhi = vhi > v3 ? vhi : v3;
-               lo += vlo;
-               hi += vhi;
-            }
-            const int slot = kcp*ncp + icp;
-            smin[slot] = lo;
-            smax[slot] = hi;
-         }
-      }
-      MFEM_SYNC_THREAD;
-
-      const int lane = ty*ncp + tx;
-      const int nactive = ncp*ncp;
-      const int nthreads = MFEM_THREAD_SIZE(x)*MFEM_THREAD_SIZE(y);
-
-      // Reduce all 2D control-point bounds to one lower/upper pair per
-      // (element, component).
-      if (nthreads == 1)
-      {
-         if (tx == 0 && ty == 0)
-         {
-            real_t lower_ec = smin[0];
-            real_t upper_ec = smax[0];
-            for (int t = 1; t < nactive; t++)
-            {
-               lower_ec = lower_ec < smin[t] ? lower_ec : smin[t];
-               upper_ec = upper_ec > smax[t] ? upper_ec : smax[t];
-            }
-            L(e, c) = lower_ec;
-            U(e, c) = upper_ec;
-         }
-      }
-      else
-      {
-         for (int stride = (nactive + 1)/2; stride > 0;
-              stride = (stride + 1)/2)
-         {
-            if (lane < stride && lane + stride < nactive)
-            {
-               smin[lane] = smin[lane] < smin[lane + stride] ?
-                            smin[lane] : smin[lane + stride];
-               smax[lane] = smax[lane] > smax[lane + stride] ?
-                            smax[lane] : smax[lane + stride];
-            }
-            MFEM_SYNC_THREAD;
-            if (stride == 1) { break; }
-         }
-
-         if (lane == 0)
-         {
-            L(e, c) = smin[0];
-            U(e, c) = smax[0];
-         }
-      }
-   });
-}
-
-} // namespace internal
-
-inline void PLBound::GetElementBoundsKernel(const int rdim, const int fes_vdim,
-                                            const Vector &e_vec,
-                                            Vector &lower, Vector &upper,
-                                            const int vdim) const
-{
-   MFEM_VERIFY(b_type != BasisType::Positive,
-               "Bernstein device bounds are not implemented.");
-   if (rdim == 3)
-   {
-      MFEM_ABORT("Device element bounds kernel currently only supports 1D/2D.");
-   }
-   MFEM_VERIFY(rdim == 1 || rdim == 2, "Invalid element dimension.");
-   MFEM_VERIFY(vdim >= -1 && vdim <= fes_vdim, "Invalid vector component.");
-   const int nd = static_cast<int>(std::pow(nb, rdim));
-   const int ne = e_vec.Size()/(nd*fes_vdim);
-   const int ncomp = (vdim > 0) ? 1 : fes_vdim;
-
-   lower.SetSize(ne*ncomp, e_vec);
-   upper.SetSize(ne*ncomp, e_vec);
-   lower.UseDevice(true);
-   upper.UseDevice(true);
-
-   if (!proj)
-   {
-      MFEM_ABORT("Device element bounds kernel currently requires projection "
-                 "enabled.");
-   }
-
-   const real_t *dxhat = xhat.Read();
-   const real_t *dwhat = what.Read();
-   const real_t *dcphat = cphat.Read();
-   const real_t *dlbound = lbound.Read();
-   const real_t *dubound = ubound.Read();
-
-   internal::PLBoundDeviceData data
-   {
-      nb,
-      ncp,
-      dxhat,
-      dwhat,
-      dcphat,
-      dlbound,
-      dubound
-   };
-
-   const int comp0 = (vdim > 0) ? (vdim - 1) : 0;
-
-   if (rdim == 1)
-   {
-      switch (nb)
-      {
-         case 2: return internal::GetElementBoundsKernel1D<2, true>(data, fes_vdim, ne,
-                                                                     e_vec, lower, upper,
-                                                                     comp0, ncomp);
-         case 3: return internal::GetElementBoundsKernel1D<3, true>(data, fes_vdim, ne,
-                                                                     e_vec, lower, upper,
-                                                                     comp0, ncomp);
-         case 4: return internal::GetElementBoundsKernel1D<4, true>(data, fes_vdim, ne,
-                                                                     e_vec, lower, upper,
-                                                                     comp0, ncomp);
-         case 5: return internal::GetElementBoundsKernel1D<5, true>(data, fes_vdim, ne,
-                                                                     e_vec, lower, upper,
-                                                                     comp0, ncomp);
-         case 6: return internal::GetElementBoundsKernel1D<6, true>(data, fes_vdim, ne,
-                                                                     e_vec, lower, upper,
-                                                                     comp0, ncomp);
-         case 7: return internal::GetElementBoundsKernel1D<7, true>(data, fes_vdim, ne,
-                                                                     e_vec, lower, upper,
-                                                                     comp0, ncomp);
-         case 8: return internal::GetElementBoundsKernel1D<8, true>(data, fes_vdim, ne,
-                                                                     e_vec, lower, upper,
-                                                                     comp0, ncomp);
-         case 9: return internal::GetElementBoundsKernel1D<9, true>(data, fes_vdim, ne,
-                                                                     e_vec, lower, upper,
-                                                                     comp0, ncomp);
-         case 10: return internal::GetElementBoundsKernel1D<10, true>(data, fes_vdim, ne,
-                                                                       e_vec, lower, upper,
-                                                                       comp0, ncomp);
-         default: return internal::GetElementBoundsKernel1D<0, true>(data, fes_vdim, ne,
-                                                                      e_vec, lower, upper,
-                                                                      comp0, ncomp);
-      }
-   }
-#define MFEM_PLBOUND_2D_DISPATCH(NB, NCP) \
-   return internal::GetElementBoundsKernel2D<NB, NCP, true>(data, fes_vdim, ne, \
-                                                            e_vec, lower, upper, \
-                                                            comp0, ncomp)
-   switch (nb)
-   {
-      case 2:
-         switch (ncp)
-         {
-            case 4: MFEM_PLBOUND_2D_DISPATCH(2, 4);
-            case 6: MFEM_PLBOUND_2D_DISPATCH(2, 6);
-            case 8: MFEM_PLBOUND_2D_DISPATCH(2, 8);
-         }
-         break;
-      case 3:
-         switch (ncp)
-         {
-            case 6: MFEM_PLBOUND_2D_DISPATCH(3, 6);
-            case 9: MFEM_PLBOUND_2D_DISPATCH(3, 9);
-            case 12: MFEM_PLBOUND_2D_DISPATCH(3, 12);
-         }
-         break;
-      case 4:
-         switch (ncp)
-         {
-            case 8: MFEM_PLBOUND_2D_DISPATCH(4, 8);
-            case 12: MFEM_PLBOUND_2D_DISPATCH(4, 12);
-            case 16: MFEM_PLBOUND_2D_DISPATCH(4, 16);
-         }
-         break;
-      case 5:
-         switch (ncp)
-         {
-            case 10: MFEM_PLBOUND_2D_DISPATCH(5, 10);
-            case 15: MFEM_PLBOUND_2D_DISPATCH(5, 15);
-            case 20: MFEM_PLBOUND_2D_DISPATCH(5, 20);
-         }
-         break;
-      case 6:
-         switch (ncp)
-         {
-            case 12: MFEM_PLBOUND_2D_DISPATCH(6, 12);
-            case 18: MFEM_PLBOUND_2D_DISPATCH(6, 18);
-            case 24: MFEM_PLBOUND_2D_DISPATCH(6, 24);
-         }
-         break;
-      case 7:
-         switch (ncp)
-         {
-            case 14: MFEM_PLBOUND_2D_DISPATCH(7, 14);
-            case 21: MFEM_PLBOUND_2D_DISPATCH(7, 21);
-            case 28: MFEM_PLBOUND_2D_DISPATCH(7, 28);
-         }
-         break;
-      case 8:
-         switch (ncp)
-         {
-            case 16: MFEM_PLBOUND_2D_DISPATCH(8, 16);
-            case 24: MFEM_PLBOUND_2D_DISPATCH(8, 24);
-            case 32: MFEM_PLBOUND_2D_DISPATCH(8, 32);
-         }
-         break;
-   }
-#undef MFEM_PLBOUND_2D_DISPATCH
-   return internal::GetElementBoundsKernel2D<0, 0, true>(data, fes_vdim, ne,
-                                                         e_vec, lower, upper,
-                                                         comp0, ncomp);
-}
-
 } // namespace mfem
 
 #endif // MFEM_BOUNDS

fem/dgmassinv_kernels.hpp

@@ -57,7 +57,7 @@ void DGMassApply(const int e,
       }
       else if (DIM == 3)
       {
-         SmemPAMassApply3D_Element<TD1D,TQ1D,NBZ,ACCUM>(e, NE, B, pa_data, x, y);
+         SmemPAMassApply3D_Element<TD1D,TQ1D,ACCUM>(e, NE, B, pa_data, x, y);
       }
       else
       {

fem/gridfunc.cpp

@@ -5252,30 +5252,6 @@ void GridFunction::GetElementBounds(const PLBound &plb,
                                     Vector &lower, Vector &upper,
                                     const int vdim) const
 {
-   if (UseDevice() && Device::Allows(Backend::DEVICE_MASK) &&
-       plb.GetBasisType() != BasisType::Positive &&
-       UsesTensorBasis(*fes))
-   {
-      const FiniteElement &fe = *fes->GetTypicalFE();
-      const int rdim = fe.GetDim();
-      const int fes_dim = fes->GetVDim();
-      const int nel = fes->GetNE();
-      const int nd = fe.GetDof();
-
-      Vector e_vec(nd*fes_dim*nel, Device::GetDeviceMemoryType());
-      e_vec.UseDevice(true);
-      const ElementRestrictionOperator *elem_restr =
-         fes->GetElementRestriction(ElementDofOrdering::LEXICOGRAPHIC);
-      MFEM_VERIFY(elem_restr != nullptr,
-                  "Element restriction is required for device bounds.");
-      elem_restr->Mult(*this, e_vec);
-
-      plb.GetElementBoundsKernel(rdim, fes_dim, e_vec, lower, upper, vdim);
-      lower.HostRead();
-      upper.HostRead();
-      return;
-   }
-
    int nel  = fes->GetNE();
    int fes_dim  = fes->GetVDim();
    lower.SetSize(nel*(vdim > 0 ? 1 :fes_dim));

fem/gslib.cpp

@@ -106,7 +106,9 @@ FindPointsGSLIB::FindPointsGSLIB()
    gsl_comm = new gslib::comm;
    cr       = new gslib::crystal;
 #ifdef MFEM_USE_MPI
-   if (!Mpi::IsInitialized()) { Mpi::Init(); }
+   int initialized = 0;
+   MPI_Initialized(&initialized);
+   if (!initialized) { MPI_Init(NULL, NULL); }
    MPI_Comm comm = MPI_COMM_WORLD;
    comm_init(gsl_comm, comm);
 #else
@@ -2622,7 +2624,9 @@ GSOPGSLIB::GSOPGSLIB(Array<long long> &ids)
    gsl_comm = new gslib::comm;
    cr       = new gslib::crystal;
 #ifdef MFEM_USE_MPI
-   if (!Mpi::IsInitialized()) { Mpi::Init(); }
+   int initialized;
+   MPI_Initialized(&initialized);
+   if (!initialized) { MPI_Init(NULL, NULL); }
    MPI_Comm comm = MPI_COMM_WORLD;
    comm_init(gsl_comm, comm);
 #else

fem/integ/bilininteg_mass_kernels.hpp

@@ -181,12 +181,6 @@ constexpr int NBZ(int D1D)
 {
    return ipow(2, D(D1D) >= 0 ? D(D1D) : 0);
 }
-constexpr int NBZ3D(int MDQ)
-{
-   return MDQ > 0 ? std::min<int>(
-             (128 + MDQ * MDQ * MDQ - 1) / (MDQ * MDQ * MDQ), 64)
-          : 1;
-}
 }
 
 // Shared memory PA Mass Diagonal 2D kernel
@@ -810,23 +804,19 @@ void PAMassApply3D_Element(const int e,
    }
 }
 
-template <int T_D1D, int T_Q1D, int TBATCH, bool ACCUMULATE = true>
-MFEM_HOST_DEVICE inline void
-SmemPAMassApply3D_Element(const int e, const int NE, const real_t *b_,
-                          const real_t *d_, const real_t *x_, real_t *y_,
-                          int d1d = 0, int q1d = 0)
+template<int T_D1D, int T_Q1D, bool ACCUMULATE = true>
+MFEM_HOST_DEVICE inline
+void SmemPAMassApply3D_Element(const int e,
+                               const int NE,
+                               const real_t *b_,
+                               const real_t *d_,
+                               const real_t *x_,
+                               real_t *y_,
+                               const int d1d = 0,
+                               const int q1d = 0)
 {
-   static_assert(TBATCH > 0, "TBATCH must be positive");
-#if defined(__CUDA_ARCH__) || defined(__HIP_DEVICE_COMPILE__)
-   constexpr int tbatch = TBATCH;
-   const int tidz = MFEM_THREAD_ID(z);
-#else
-   // host always batch size 1
-   constexpr int tbatch = 1;
-   constexpr int tidz = 0;
-#endif
-   const int D1D = T_D1D ? T_D1D : d1d;
-   const int Q1D = T_Q1D ? T_Q1D : q1d;
+   constexpr int D1D = T_D1D ? T_D1D : d1d;
+   constexpr int Q1D = T_Q1D ? T_Q1D : q1d;
    constexpr int MQ1 = T_Q1D ? T_Q1D : DofQuadLimits::MAX_Q1D;
    constexpr int MD1 = T_D1D ? T_D1D : DofQuadLimits::MAX_D1D;
    constexpr int MDQ = (MQ1 > MD1) ? MQ1 : MD1;
@@ -839,37 +829,33 @@ SmemPAMassApply3D_Element(const int e, const int NE, const real_t *b_,
    MFEM_SHARED real_t sDQ[MQ1*MD1];
    real_t (*B)[MD1] = (real_t (*)[MD1]) sDQ;
    real_t (*Bt)[MQ1] = (real_t (*)[MQ1]) sDQ;
-   MFEM_SHARED real_t sm0[tbatch][MDQ*MDQ*MDQ];
-   MFEM_SHARED real_t sm1[tbatch][MDQ*MDQ*MDQ];
-   real_t (*X)[MD1][MD1]   = (real_t (*)[MD1][MD1]) (sm0+tidz);
-   real_t (*DDQ)[MD1][MQ1] = (real_t (*)[MD1][MQ1]) (sm1+tidz);
-   real_t (*DQQ)[MQ1][MQ1] = (real_t (*)[MQ1][MQ1]) (sm0+tidz);
-   real_t (*QQQ)[MQ1][MQ1] = (real_t (*)[MQ1][MQ1]) (sm1+tidz);
-   real_t (*QQD)[MQ1][MD1] = (real_t (*)[MQ1][MD1]) (sm0+tidz);
-   real_t (*QDD)[MD1][MD1] = (real_t (*)[MD1][MD1]) (sm1+tidz);
-   MFEM_FOREACH_THREAD(dy, y, D1D)
+   MFEM_SHARED real_t sm0[MDQ*MDQ*MDQ];
+   MFEM_SHARED real_t sm1[MDQ*MDQ*MDQ];
+   real_t (*X)[MD1][MD1]   = (real_t (*)[MD1][MD1]) sm0;
+   real_t (*DDQ)[MD1][MQ1] = (real_t (*)[MD1][MQ1]) sm1;
+   real_t (*DQQ)[MQ1][MQ1] = (real_t (*)[MQ1][MQ1]) sm0;
+   real_t (*QQQ)[MQ1][MQ1] = (real_t (*)[MQ1][MQ1]) sm1;
+   real_t (*QQD)[MQ1][MD1] = (real_t (*)[MQ1][MD1]) sm0;
+   real_t (*QDD)[MD1][MD1] = (real_t (*)[MD1][MD1]) sm1;
+   MFEM_FOREACH_THREAD(dy,y,D1D)
    {
-      MFEM_FOREACH_THREAD(dx, x, D1D)
+      MFEM_FOREACH_THREAD(dx,x,D1D)
       {
          MFEM_UNROLL(MD1)
          for (int dz = 0; dz < D1D; ++dz)
          {
-            X[dz][dy][dx] = x(dx, dy, dz, e);
+            X[dz][dy][dx] = x(dx,dy,dz,e);
          }
       }
-      MFEM_FOREACH_THREAD(dx, x, Q1D) { B[dx][dy] = b(dx, dy); }
-   }
-   if (tidz == 0)
-   {
-      MFEM_FOREACH_THREAD(dy, y, D1D)
+      MFEM_FOREACH_THREAD(dx,x,Q1D)
       {
-         MFEM_FOREACH_THREAD(dx, x, Q1D) { B[dx][dy] = b(dx, dy); }
+         B[dx][dy] = b(dx,dy);
       }
    }
    MFEM_SYNC_THREAD;
-   MFEM_FOREACH_THREAD(dy, y, D1D)
+   MFEM_FOREACH_THREAD(dy,y,D1D)
    {
-      MFEM_FOREACH_THREAD(qx, x, Q1D)
+      MFEM_FOREACH_THREAD(qx,x,Q1D)
       {
          real_t u[D1D];
          MFEM_UNROLL(MD1)
@@ -894,9 +880,9 @@ SmemPAMassApply3D_Element(const int e, const int NE, const real_t *b_,
       }
    }
    MFEM_SYNC_THREAD;
-   MFEM_FOREACH_THREAD(qy, y, Q1D)
+   MFEM_FOREACH_THREAD(qy,y,Q1D)
    {
-      MFEM_FOREACH_THREAD(qx, x, Q1D)
+      MFEM_FOREACH_THREAD(qx,x,Q1D)
       {
          real_t u[D1D];
          MFEM_UNROLL(MD1)
@@ -921,9 +907,9 @@ SmemPAMassApply3D_Element(const int e, const int NE, const real_t *b_,
       }
    }
    MFEM_SYNC_THREAD;
-   MFEM_FOREACH_THREAD(qy, y, Q1D)
+   MFEM_FOREACH_THREAD(qy,y,Q1D)
    {
-      MFEM_FOREACH_THREAD(qx, x, Q1D)
+      MFEM_FOREACH_THREAD(qx,x,Q1D)
       {
          real_t u[Q1D];
          MFEM_UNROLL(MQ1)
@@ -943,22 +929,22 @@ SmemPAMassApply3D_Element(const int e, const int NE, const real_t *b_,
          MFEM_UNROLL(MQ1)
          for (int qz = 0; qz < Q1D; qz++)
          {
-            QQQ[qz][qy][qx] = u[qz] * d(qx, qy, qz, e);
+            QQQ[qz][qy][qx] = u[qz] * d(qx,qy,qz,e);
          }
       }
    }
    MFEM_SYNC_THREAD;
-   if (tidz == 0)
+   MFEM_FOREACH_THREAD(di,y,D1D)
    {
-      MFEM_FOREACH_THREAD(di, y, D1D)
+      MFEM_FOREACH_THREAD(q,x,Q1D)
       {
-         MFEM_FOREACH_THREAD(q, x, Q1D) { Bt[di][q] = b(q, di); }
+         Bt[di][q] = b(q,di);
       }
    }
    MFEM_SYNC_THREAD;
-   MFEM_FOREACH_THREAD(qy, y, Q1D)
+   MFEM_FOREACH_THREAD(qy,y,Q1D)
    {
-      MFEM_FOREACH_THREAD(dx, x, D1D)
+      MFEM_FOREACH_THREAD(dx,x,D1D)
       {
          real_t u[Q1D];
          MFEM_UNROLL(MQ1)
@@ -983,9 +969,9 @@ SmemPAMassApply3D_Element(const int e, const int NE, const real_t *b_,
       }
    }
    MFEM_SYNC_THREAD;
-   MFEM_FOREACH_THREAD(dy, y, D1D)
+   MFEM_FOREACH_THREAD(dy,y,D1D)
    {
-      MFEM_FOREACH_THREAD(dx, x, D1D)
+      MFEM_FOREACH_THREAD(dx,x,D1D)
       {
          real_t u[Q1D];
          MFEM_UNROLL(MQ1)
@@ -1010,9 +996,9 @@ SmemPAMassApply3D_Element(const int e, const int NE, const real_t *b_,
       }
    }
    MFEM_SYNC_THREAD;
-   MFEM_FOREACH_THREAD(dy, y, D1D)
+   MFEM_FOREACH_THREAD(dy,y,D1D)
    {
-      MFEM_FOREACH_THREAD(dx, x, D1D)
+      MFEM_FOREACH_THREAD(dx,x,D1D)
       {
          real_t u[D1D];
          MFEM_UNROLL(MD1)
@@ -1034,11 +1020,11 @@ SmemPAMassApply3D_Element(const int e, const int NE, const real_t *b_,
          {
             if (ACCUMULATE)
             {
-               y(dx, dy, dz, e) += u[dz];
+               y(dx,dy,dz,e) += u[dz];
             }
             else
             {
-               y(dx, dy, dz, e) = u[dz];
+               y(dx,dy,dz,e) = u[dz];
             }
          }
       }
@@ -1129,8 +1115,8 @@ inline void PAMassApply3D(const int NE,
    });
 }
 
-// Shared memory PA Mass Apply 3D kernel
-template<int T_D1D = 0, int T_Q1D = 0, int TBATCH=1>
+// Shared memory PA Mass Apply 2D kernel
+template<int T_D1D = 0, int T_Q1D = 0>
 inline void SmemPAMassApply3D(const int NE,
                               const Array<real_t> &b_,
                               const Array<real_t> &bt_,
@@ -1140,9 +1126,6 @@ inline void SmemPAMassApply3D(const int NE,
                               const int d1d = 0,
                               const int q1d = 0)
 {
-   static_assert(T_D1D > 0, "T_D1D must be positive");
-   static_assert(T_Q1D > 0, "T_Q1D must be positive");
-   static_assert(TBATCH > 0, "TBATCH must be positive");
    MFEM_CONTRACT_VAR(bt_);
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -1154,11 +1137,9 @@ inline void SmemPAMassApply3D(const int NE,
    const auto d = d_.Read();
    const auto x = x_.Read();
    auto y = y_.ReadWrite();
-   mfem::forall_2D_batch<T_Q1D * T_Q1D * TBATCH>(NE, Q1D, Q1D, TBATCH,
-                                                 [=] MFEM_HOST_DEVICE(int e)
+   mfem::forall_2D<T_Q1D*T_Q1D>(NE, Q1D, Q1D, [=] MFEM_HOST_DEVICE (int e)
    {
-      internal::SmemPAMassApply3D_Element<T_D1D, T_Q1D, TBATCH>(e, NE, b, d, x,
-                                                                y, d1d, q1d);
+      internal::SmemPAMassApply3D_Element<T_D1D,T_Q1D>(e, NE, b, d, x, y, d1d, q1d);
    });
 }
 
@@ -1413,16 +1394,7 @@ ApplyKernelType MassIntegrator::ApplyPAKernels::Kernel()
 {
    if constexpr (DIM == 1) { return internal::PAMassApply1D; }
    else if constexpr (DIM == 2) { return internal::SmemPAMassApply2D<T_D1D,T_Q1D>; }
-   else if constexpr (DIM == 3)
-   {
-      constexpr int MDQ = T_D1D >= T_Q1D ? T_D1D : T_Q1D;
-      // max 64 threads in z limit in cuda and hip
-      if constexpr (MDQ > 0)
-      {
-         return internal::SmemPAMassApply3D<T_D1D, T_Q1D,
-                internal::mass::NBZ3D(MDQ)>;
-      }
-   }
+   else if constexpr (DIM == 3) { return internal::SmemPAMassApply3D<T_D1D, T_Q1D>; }
    MFEM_ABORT("");
 }
 

general/forall.hpp

@@ -1090,12 +1090,6 @@ inline void forall_2D_batch(int N, int X, int Y, int BZ, lambda &&body)
    ForallWrap<2>(true, N, body, X, Y, BZ);
 }
 
-template<int MAX_THREADS_PER_BLOCK, typename lambda>
-inline void forall_2D_batch(int N, int X, int Y, int BZ, lambda &&body)
-{
-   ForallWrap<2, MAX_THREADS_PER_BLOCK>(true, N, body, X, Y, BZ);
-}
-
 template<typename lambda>
 inline void forall_3D(int N, int X, int Y, int Z, lambda &&body)
 {

linalg/CMakeLists.txt

@@ -23,6 +23,7 @@ list(APPEND SRCS
   complex_operator.cpp
   constraints.cpp
   densemat.cpp
+  eigensolvers.cpp
   filteredsolver.cpp
   handle.cpp
   matrix.cpp
@@ -55,6 +56,7 @@ list(APPEND HDRS
   dinvariants.hpp
   dtensor.hpp
   dual.hpp
+  eigensolvers.hpp
   filteredsolver.hpp
   handle.hpp
   invariants.hpp
@@ -101,6 +103,11 @@ if (MFEM_USE_MPI)
   endif()
 endif()
 
+if (MFEM_USE_ARPACK)
+  list(APPEND SRCS arpack.cpp)
+  list(APPEND HDRS arpack.hpp)
+endif()
+
 if (MFEM_USE_SUNDIALS)
   list(APPEND SRCS sundials.cpp)
   list(APPEND HDRS sundials.hpp)

linalg/arpack.hpp

@@ -0,0 +1,271 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef MFEM_ARPACK
+#define MFEM_ARPACK
+
+#include "../config/config.hpp"
+
+#ifdef MFEM_USE_ARPACK
+
+#include <string>
+
+#ifdef MFEM_USE_MPI
+#include <mpi.h>
+#include "hypre.hpp"
+#endif
+
+#include "operator.hpp"
+
+#define SSAUPD  ssaupd_
+#define SSEUPD  sseupd_
+#define DSAUPD  dsaupd_
+#define DSEUPD  dseupd_
+
+#ifdef MFEM_USE_MPI
+#define PSSAUPD pssaupd_
+#define PSSEUPD psseupd_
+#define PDSAUPD pdsaupd_
+#define PDSEUPD pdseupd_
+#endif
+
+extern "C" void SSAUPD(int *ido, char *bmat, int *n,
+                       char *which, int *nev, float *tol, float *resid,
+                       int *ncv, float *v, int *ldv,
+                       int *iparam, int *ipntr,
+                       float *workd, float *workl, int *lworkl, int *info);
+
+extern "C" void SSEUPD(int *, char *, int *, float *,
+                       float *, int *, float *, char *, int *, char *,
+                       int *, float *, float *, int *, float *,
+                       int *, int  *, int *, float *,
+                       float *, int *, int *);
+
+extern "C" void DSAUPD(int *ido, char *bmat, int *n,
+                       char *which, int *nev, double *tol, double *resid,
+                       int *ncv, double *v, int *ldv,
+                       int *iparam, int *ipntr,
+                       double *workd, double *workl, int *lworkl, int *info);
+
+extern "C" void DSEUPD(int *, char *, int *, double *,
+                       double *, int *, double *, char *, int *, char *,
+                       int *, double *, double *, int *, double *,
+                       int *, int  *, int *, double *,
+                       double *, int *, int *);
+
+#ifdef MFEM_USE_MPI
+
+extern "C" void PSSAUPD(int *comm, int *ido, char *bmat, int *n,
+                        char *which, int *nev, float *tol, float *resid,
+                        int *ncv, float *v, int *ldv,
+                        int *iparam, int *ipntr,
+                        float *workd, float *workl, int *lworkl, int *info);
+
+extern "C" void PSSEUPD(int *comm, int *, char *, int *, float *,
+                        float *, int *, float *, char *, int *, char *,
+                        int *, float *, float *, int *, float *,
+                        int *, int  *, int *, float *,
+                        float *, int *, int *);
+
+extern "C" void PDSAUPD(int *comm, int *ido, char *bmat, int *n,
+                        char *which, int *nev, double *tol, double *resid,
+                        int *ncv, double *v, int *ldv,
+                        int *iparam, int *ipntr,
+                        double *workd, double *workl, int *lworkl, int *info);
+
+extern "C" void PDSEUPD(int *comm, int *, char *, int *, double *,
+                        double *, int *, double *, char *, int *, char *,
+                        int *, double *, double *, int *, double *,
+                        int *, int  *, int *, double *,
+                        double *, int *, int *);
+
+#endif
+
+extern "C" {
+   void arpackgetcommdbg_(int *,int *,int *);
+   void arpacksetcommdbg_(int *,int *,int *);
+   void arpacksymdbg_(int *,int *,int *,int *,int *,int *,int *);
+   void arpacknonsymdbg_(int *,int *,int *,int *,int *,int *,int *);
+   void arpackcmplxdbg_(int *,int *,int *,int *,int *,int *,int *);
+}
+
+namespace mfem
+{
+
+/// Wrapper for the ARPACK routine SSAUPD or DSAUPD
+class ArPackSAUPD : public SymEigensolver, public SymGenEigensolver
+{
+public:
+
+   ArPackSAUPD();
+   virtual ~ArPackSAUPD();
+
+   /** ARPACK modes are described in section 3.5 of the ARPACK manual.
+       Mode 1: regular mode to solve A x = lambda x
+               No solver and no mass matrix are needed.
+       Mode 2: regular inverse mode to solve A x = lambda M x
+               Both A and M are needed and the solver should compute M^{-1}.
+       Mode 3: shift-invert mode to solve either A x = lambda x
+               or A x = lambda M x
+               Mass matrix is optional.  The solver should compute
+               (A-sigma I)^{-1} or (A-sigma M)^{-1}.  The shift parameter,
+               sigma, also needs to be set with SetShift().
+       Mode 4: Buckling mode to solve K x = lambda K_G x
+               K is set using SetMassMatrix(), K_G is set using SetOperator(),
+          and the solver should compute (K-sigma K_G)^{-1}.  The shift
+               parameter, sigma, also needs to be set with SetShift().
+       Mode 5: Cayley mode to solve A x = lambda M x
+               Both A and M are needed and the solver should compute
+               (A - sigma M)^{-1}.  The shift parameter, sigma, also needs
+               to be set with SetShift().
+    */
+   void SetMode(int mode);
+
+   inline void SetTol(real_t tol) override         {      tol_ = tol;      }
+   inline void SetMaxIter(int max_iter) override   { max_iter_ = max_iter; }
+   inline void SetPrintLevel(int logging) override {  logging_ = logging;  }
+   inline void SetShift(real_t sigma)              {    sigma_ = sigma;    }
+   inline void SetNumModes(int num_eigs) override  {      nev_ = num_eigs; }
+
+   virtual void SetSolver(Solver & solver);
+   virtual void SetOperator(const Operator & A) override;
+   virtual void SetMassMatrix(const Operator & M);
+   virtual void SetOperators(const Operator & A, const Operator & B) override
+   { SetOperator(A); SetMassMatrix(B); }
+
+   void Solve() override;
+
+   virtual int GetNumConverged() const override { return iparam_[4]; }
+
+   /// Collect the converged eigenvalues
+   virtual void GetEigenvalues(Array<real_t> & eigenvalues) const override;
+
+   /// Extract a single eigenvector
+   virtual const Vector & GetEigenvector(unsigned int i) const override;
+
+   /// Transfer ownership of the converged eigenvectors
+   Vector ** StealEigenvectors() override;
+
+protected:
+
+   int myid_;       // Index of this processor
+   int max_iter_;
+   int logging_;
+
+   // The following variables are for ARPACK
+   int nloc_;       // number of items stored locally
+   int nev_;        // number of requested eigenvalues
+   int ncv_;        // number of ritz vectors
+   int rvec_;       // boolean to return eigenvectors as well
+   int mode_;       // 1 = standard, 2 = generalized, 3 = shift invert,
+   // 4 = buckling, 5 = Cayley
+   int lworkl_;     // length of lworkl_ work array
+   int iparam_[12]; // arpack parameters
+   int ipntr_[12];  // arpack pointers
+
+   char bmat_;      // I for standard problem, G for generalized
+   char which_[3];  // spectrum portion: LA, SA, LM, SM, BE
+   char hwmny_;     // DSEUPD: A for all eigenvalues, S for some
+
+   real_t tol_;     // relative accuracy bound for Ritz values
+   real_t sigma_;   // eigenvalue shift parameter
+
+   int    * select_;// workspace used during eigenvalue computation
+   real_t * dv_;    // Ritz values
+   real_t * v_;     // ncv Lanczos basis vectors
+   real_t * resid_; // residual vector
+   real_t * workd_; // work array for 3 vectors used in Arnoldi iteration
+   real_t * workl_; // work array
+
+   // Operators and Vectors needed outside of ARPACK
+   Solver   * solver_;
+   const Operator * A_;
+   const Operator * B_;
+
+   Vector * w_;
+   Vector * x_;
+   Vector * y_;
+   Vector * z_;
+
+   mutable Vector ** eigenvectors_;
+
+   std::string solverName_;
+
+   void reverseComm();
+
+   int reverseCommMode1();
+   int reverseCommMode2();
+   int reverseCommMode3();
+   int reverseCommMode4();
+   int reverseCommMode5();
+
+   virtual void prepareEigenvectors() const;
+
+   void printErrors(const int & info, const int iparam[],
+                    const char & bmat, const int & n,
+                    const char which[],
+                    const int & nev, const int & ncv,
+                    const int & lworkl );
+
+private:
+
+   virtual int computeNlocf() { return nloc_; }
+   virtual int computeIter(int & ido);
+   virtual int computeEigs();
+
+};
+
+#ifdef MFEM_USE_MPI
+
+class ArPackPSAUPD : public ArPackSAUPD
+{
+public:
+   ArPackPSAUPD(MPI_Comm comm);
+   virtual ~ArPackPSAUPD() {}
+
+   void SetOperator(const Operator & A);
+   void SetMassMatrix(const Operator & M);
+
+   /// Collect the converged eigenvalues
+   void GetEigenvalues(Array<real_t> & eigenvalues) const;
+
+   /// Extract a single eigenvector
+   const Vector & GetEigenvector(unsigned int i) const;
+
+   /// Transfer ownership of the converged eigenvectors
+   // HypreParVector ** StealEigenvectors();
+   Vector ** StealEigenvectors();
+
+protected:
+
+   void prepareEigenvectors() const;
+
+private:
+
+   MPI_Comm comm_;
+   MPI_Fint commf_; // Fortran style MPI communicator
+   int numProcs_;   // Number of processors
+
+   mutable HYPRE_Int * part_; // parallel partitioning for eigenvectors
+
+   int computeNlocf();
+   int computeIter(int & ido);
+   int computeEigs();
+
+};
+
+#endif // MFEM_USE_MPI
+
+};
+
+#endif // MFEM_USE_ARPACK
+
+#endif // MFEM_ARPACK

linalg/eigensolvers.cpp

@@ -0,0 +1,20 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "linalg.hpp"
+#include "eigensolvers.hpp"
+
+using namespace std;
+
+namespace mfem
+{
+
+};

linalg/eigensolvers.hpp

@@ -0,0 +1,396 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef MFEM_EIGENSOLVERS
+#define MFEM_EIGENSOLVERS
+
+#include "vector.hpp"
+
+namespace mfem
+{
+
+/// Abstract Eigenequation
+/// Defines the operator of the linear eigenvalue equation
+///    A x_i = lambda_i x_i
+/// Where A is a real-valued operator, the lambda_i are the eigenvalues,
+/// and x_i are the eigenvectors.
+class Eigenequation
+{
+protected:
+   Eigenequation() = default;
+
+public:
+   virtual ~Eigenequation() = default;
+
+   /// @brief Set the operator A of the eigenvalue equation
+   virtual void SetOperator(const Operator & A) = 0;
+};
+
+/// Abstract Complex-valued Eigenequation
+/// Defines the operator of the linear eigenvalue equation
+///    A x_i = lambda_i x_i
+/// Where A is a complex-valued operator, the lambda_i are the eigenvalues,
+/// and x_i are the eigenvectors.
+class ComplexEigenequation
+{
+protected:
+   ComplexEigenequation() = default;
+
+public:
+   virtual ~ComplexEigenequation() = default;
+
+   /// @brief Set the real and imaginary parts of the operator A
+   virtual void SetOperator(const Operator & Ar, const Operator & Ai) = 0;
+};
+
+/// Abstract Generalized Eigenequation
+/// Defines the operator of the linear eigenvalue equation
+///    A x_i = lambda_i B x_i
+/// Where A and B are real-valued operators, the lambda_i are the eigenvalues,
+/// and x_i are the eigenvectors.
+class GenEigenequation
+{
+protected:
+   GenEigenequation() = default;
+
+public:
+   virtual ~GenEigenequation() = default;
+
+   /// @brief Set the operators A and B of the generalized eigenvalue equation
+   virtual void SetOperators(const Operator & A, const Operator & B) = 0;
+};
+
+/// Abstract Complex-valued Generalized Eigenequation
+/// Defines the operator of the linear eigenvalue equation
+///    A x_i = lambda_i B x_i
+/// Where A and B are complex-valued operators, the lambda_i are the
+/// eigenvalues, and x_i are the eigenvectors.
+class ComplexGenEigenequation
+{
+protected:
+   ComplexGenEigenequation() = default;
+
+public:
+   virtual ~ComplexGenEigenequation() = default;
+
+   /// @brief Set the real and imaginary parts of the operators A and B
+   virtual void SetOperators(const Operator & Ar, const Operator & Ai,
+                             const Operator & Br, const Operator & Bi) = 0;
+};
+
+/// Abstract Eigensolver
+/// Computes eigenvalue/eigenvector pairs for the linear system
+///    A x_i = lambda_i x_i
+/// Where the lambda_i are the eigenvalues and x_i are the eigenvectors.
+class EigensolverBase
+{
+protected:
+   EigensolverBase() = default;
+
+public:
+   virtual ~EigensolverBase() = default;
+
+   /// @brief Stopping criteria based on numerical tolerance
+   ///
+   /// @note This may be defined differently by different solvers.
+   virtual void SetTol(real_t tol) = 0;
+
+   /// @brief Stopping criteria based on number of iterations required to
+   /// reach convergence.
+   ///
+   /// @note This may also be defined differently in different solvers.
+   virtual void SetMaxIter(int max_iter) = 0;
+
+   /// @brief Controls the type and amount of information printed to
+   ///        standard output.
+   virtual void SetPrintLevel(int logging) = 0;
+
+   /// @brief Set the number of desired eigenmodes to compute
+   virtual void SetNumModes(int num_eigs) = 0;
+
+   /// @brief Get the number of converged eigenmodes
+   virtual int GetNumConverged() const = 0;
+
+   /// @brief Perform the eigenvalue solve
+   virtual void Solve() = 0;
+};
+
+/// Symmetric Eigensolver
+/// If A^T = A the linear system must have real-valued eigenvalues
+/// and eigenvectors.
+class SymEigensolver : public EigensolverBase, public Eigenequation
+{
+protected:
+   SymEigensolver() = default;
+
+public:
+   virtual ~SymEigensolver() = default;
+
+   /// @brief Collect the converged eigenvalues
+   ///
+   /// The length of the array should equal the number of converged eigenvalues.
+   virtual void GetEigenvalues(Array<real_t> & eigenvalues) const = 0;
+
+   /// @brief Extract a single eigenvector
+   ///
+   /// The index i should be in the range [0, numConverged). The
+   virtual const Vector & GetEigenvector(unsigned int i) const = 0;
+
+   /// @brief Transfer ownership of the converged eigenvectors
+   ///
+   /// The array should contain numConverged vectors.
+   virtual Vector ** StealEigenvectors() = 0;
+};
+
+/// Symmetric Generalized Eigensolver
+/// If A^T = A and M^T = M the linear system must have real-valued eigenvalues
+/// and eigenvectors.
+class SymGenEigensolver : public EigensolverBase, public GenEigenequation
+{
+protected:
+   SymGenEigensolver() = default;
+
+public:
+   virtual ~SymGenEigensolver() = default;
+
+   /// @brief Collect the converged eigenvalues
+   ///
+   /// The length of the array should equal the number of converged eigenvalues.
+   virtual void GetEigenvalues(Array<real_t> & eigenvalues) const = 0;
+
+   /// @brief Extract a single eigenvector
+   ///
+   /// The index i should be in the range [0, numConverged). The
+   virtual const Vector & GetEigenvector(unsigned int i) const = 0;
+
+   /// @brief Transfer ownership of the converged eigenvectors
+   ///
+   /// The array should contain numConverged vectors.
+   virtual Vector ** StealEigenvectors() = 0;
+};
+
+/// Hermetian Eigensolver
+/// If A^H = A the linear system must have real-valued eigenvalues
+/// but may have complex-valued eigenvectors.
+class HermEigensolver : public EigensolverBase, public ComplexEigenequation
+{
+protected:
+   HermEigensolver() = default;
+
+public:
+   virtual ~HermEigensolver() = default;
+
+   /// @brief Collect the converged eigenvalues
+   ///
+   /// The length of the array should equal the number of converged eigenvalues.
+   virtual void GetEigenvalues(Array<real_t> & eigenvalues) const = 0;
+
+   /// @brief Extract a single eigenvector
+   ///
+   /// The index i should be in the range [0, 2*numConverged). The
+   /// vectors corresponding to even indices are the real parts of the
+   /// converged eigenvectors and the odd indices correspond to the
+   /// imaginary parts.
+   virtual const Vector & GetEigenvector(unsigned int i) const = 0;
+
+   /// @brief Transfer ownership of the converged eigenvectors
+   ///
+   /// The array should contain 2*numConverged vectors with the even
+   /// indices corresponding to the real parts of the converged
+   /// eigenvectors and the odd indices corresponding to the imaginary
+   /// parts.
+   virtual Vector ** StealEigenvectors() = 0;
+};
+
+/// Hermetian Generalized Eigensolver
+/// If A^H = A and M^H = M the linear system must have real-valued eigenvalues
+/// but may have complex-valued eigenvectors.
+class HermGenEigensolver :
+   public EigensolverBase, public ComplexGenEigenequation
+{
+protected:
+   HermGenEigensolver() = default;
+
+public:
+   virtual ~HermGenEigensolver() = default;
+
+   /// @brief Collect the converged eigenvalues
+   ///
+   /// The length of the array should equal the number of converged eigenvalues.
+   virtual void GetEigenvalues(Array<real_t> & eigenvalues) const = 0;
+
+   /// @brief Extract a single eigenvector
+   ///
+   /// The index i should be in the range [0, 2*numConverged). The
+   /// vectors corresponding to even indices are the real parts of the
+   /// converged eigenvectors and the odd indices correspond to the
+   /// imaginary parts.
+   virtual const Vector & GetEigenvector(unsigned int i) const = 0;
+
+   /// @brief Transfer ownership of the converged eigenvectors
+   ///
+   /// The array should contain 2*numConverged vectors with the even
+   /// indices corresponding to the real parts of the converged
+   /// eigenvectors and the odd indices corresponding to the imaginary
+   /// parts.
+   virtual Vector ** StealEigenvectors() = 0;
+};
+
+/// Non-Symmetric Eigensolver
+/// For general real-valued operators A the linear system must have
+/// eigenvalues and eigenvectors which form complex conjugate pairs.
+class NonSymEigensolver : public EigensolverBase, public Eigenequation
+{
+protected:
+   NonSymEigensolver() = default;
+
+public:
+   virtual ~NonSymEigensolver() = default;
+
+   /// @brief Collect the converged eigenvalues
+   ///
+   /// The length of the array should be the number of converged
+   /// eigenvalues. The complex-valued eigenvalues can be constructed
+   /// as: lambda_{2*j} = eig[2*j]+i*eig[2*j+1] and
+   /// lambda_{2*j+1} = eig[2*j]-i*eig[2*j+1]
+   /// With j in the range [0, numConverged/2)
+   virtual void GetEigenvalues(Array<real_t> & eig) const = 0;
+
+   /// @brief Extract a single eigenvector
+   ///
+   /// The index i should be in the range [0, numConverged). The
+   /// vectors corresponding to even indices are the real parts of the
+   /// converged eigenvectors and the odd indices correspond to the
+   /// imaginary parts. If needed, the complex conjugate pairs of
+   /// eigenvectors can be constructed in the same manner described
+   /// for the eigenvalues.
+   virtual const Vector & GetEigenvector(unsigned int i) const = 0;
+
+   /// @brief Transfer ownership of the converged eigenvectors
+   ///
+   /// The array should contain numConverged vectors with the even
+   /// indices corresponding to the real parts of the converged
+   /// eigenvectors and the odd indices corresponding to the imaginary
+   /// parts.
+   virtual Vector ** StealEigenvectors() = 0;
+};
+
+/// Non-Symmetric Eigensolver
+/// For general real-valued operators A and M the linear system must have
+/// eigenvalues and eigenvectors which form complex conjugate pairs.
+class NonSymGenEigensolver : public EigensolverBase, public GenEigenequation
+{
+protected:
+   NonSymGenEigensolver() = default;
+
+public:
+   virtual ~NonSymGenEigensolver() = default;
+
+   /// @brief Collect the converged eigenvalues
+   ///
+   /// The length of the array should be the number of converged
+   /// eigenvalues. The complex-valued eigenvalues can be constructed
+   /// as: lambda_{2*j} = eig[2*j]+i*eig[2*j+1] and
+   /// lambda_{2*j+1} = eig[2*j]-i*eig[2*j+1]
+   /// With j in the range [0, numConverged/2)
+   virtual void GetEigenvalues(Array<real_t> & eig) const = 0;
+
+   /// @brief Extract a single eigenvector
+   ///
+   /// The index i should be in the range [0, numConverged). The
+   /// vectors corresponding to even indices are the real parts of the
+   /// converged eigenvectors and the odd indices correspond to the
+   /// imaginary parts. If needed, the complex conjugate pairs of
+   /// eigenvectors can be constructed in the same manner described
+   /// for the eigenvalues.
+   virtual const Vector & GetEigenvector(unsigned int i) const = 0;
+
+   /// @brief Transfer ownership of the converged eigenvectors
+   ///
+   /// The array should contain numConverged vectors with the even
+   /// indices corresponding to the real parts of the converged
+   /// eigenvectors and the odd indices corresponding to the imaginary
+   /// parts.
+   virtual Vector ** StealEigenvectors() = 0;
+};
+
+/// Complex Eigensolver
+/// Can have arbitrary complex-valued eigenvalues and eigenvectors
+class ComplexEigensolver : public EigensolverBase, public ComplexEigenequation
+{
+protected:
+   ComplexEigensolver() = default;
+
+public:
+   virtual ~ComplexEigensolver() = default;
+
+   /// @brief Collect the converged eigenvalues
+   ///
+   /// The length of the array should be twice the number of converged
+   /// eigenvalues. The complex-valued eigenvalues can be constructed
+   /// as: lambda_j = eig[2*j]+i*eig[2*j+1]
+   virtual void GetEigenvalues(Array<real_t> & eig) const = 0;
+
+   /// @brief Extract a single eigenvector
+   ///
+   /// The index i should be in the range [0, 2*numConverged). The
+   /// vectors corresponding to even indices are the real parts of the
+   /// converged eigenvectors and the odd indices correspond to the
+   /// imaginary parts.
+   virtual const Vector & GetEigenvector(unsigned int i) const = 0;
+
+   /// @brief Transfer ownership of the converged eigenvectors
+   ///
+   /// The array should contain 2*numConverged vectors with the even
+   /// indices corresponding to the real parts of the converged
+   /// eigenvectors and the odd indices corresponding to the imaginary
+   /// parts.
+   virtual Vector ** StealEigenvectors() = 0;
+};
+
+/// Complex Generalized Eigensolver
+/// Can have arbitrary complex-valued eigenvalues and eigenvectors
+class ComplexGenEigensolver :
+   public EigensolverBase, public ComplexGenEigenequation
+{
+protected:
+   ComplexGenEigensolver() = default;
+
+public:
+   virtual ~ComplexGenEigensolver() = default;
+
+   /// @brief Collect the converged eigenvalues
+   ///
+   /// The length of the array should be twice the number of converged
+   /// eigenvalues. The complex-valued eigenvalues can be constructed
+   /// as: lambda_j = eig[2*j]+i*eig[2*j+1]
+   virtual void GetEigenvalues(Array<real_t> & eig) const = 0;
+
+   /// @brief Extract a single eigenvector
+   ///
+   /// The index i should be in the range [0, 2*numConverged). The
+   /// vectors corresponding to even indices are the real parts of the
+   /// converged eigenvectors and the odd indices correspond to the
+   /// imaginary parts.
+   virtual const Vector & GetEigenvector(unsigned int i) const = 0;
+
+   /// @brief Transfer ownership of the converged eigenvectors
+   ///
+   /// The array should contain 2*numConverged vectors with the even
+   /// indices corresponding to the real parts of the converged
+   /// eigenvectors and the odd indices corresponding to the imaginary
+   /// parts.
+   virtual Vector ** StealEigenvectors() = 0;
+};
+
+}
+
+#endif

linalg/hypre.cpp

@@ -6556,7 +6556,7 @@ HypreLOBPCG::SetPreconditioner(Solver & precond)
 }
 
 void
-HypreLOBPCG::SetOperator(Operator & A)
+HypreLOBPCG::SetOperator(const Operator & A)
 {
    HYPRE_BigInt locSize = A.Width();
 
@@ -6603,7 +6603,7 @@ HypreLOBPCG::SetOperator(Operator & A)
 }
 
 void
-HypreLOBPCG::SetMassMatrix(Operator & M)
+HypreLOBPCG::SetMassMatrix(const Operator & M)
 {
    matvec_fn.MatvecCreate  = this->OperatorMatvecCreate;
    matvec_fn.Matvec        = this->OperatorMatvec;
@@ -6624,7 +6624,7 @@ HypreLOBPCG::GetEigenvalues(Array<real_t> & eigs) const
    }
 }
 
-const HypreParVector &
+const Vector &
 HypreLOBPCG::GetEigenvector(unsigned int i) const
 {
    return multi_vec->GetVector(i);
@@ -6866,6 +6866,24 @@ HypreAME::SetPreconditioner(HypreSolver & precond)
    ams_precond = &precond;
 }
 
+void
+HypreAME::SetOperators(const Operator & opA, const Operator & opB)
+{
+   const HypreParMatrix * A = dynamic_cast<const HypreParMatrix *>(&opA);
+   if (A == NULL)
+   {
+      mfem_error("HypreAME::SetOperator : first operator not HypreParMatrix!");
+   }
+   SetOperator(*A);
+
+   const HypreParMatrix * B = dynamic_cast<const HypreParMatrix *>(&opB);
+   if (B == NULL)
+   {
+      mfem_error("HypreAME::SetOperator : second operator not HypreParMatrix!");
+   }
+   SetMassMatrix(*B);
+}
+
 void
 HypreAME::SetOperator(const HypreParMatrix & A)
 {
@@ -6924,7 +6942,7 @@ HypreAME::createDummyVectors() const
    }
 }
 
-const HypreParVector &
+const Vector &
 HypreAME::GetEigenvector(unsigned int i) const
 {
    if ( eigenvectors == NULL )
@@ -6935,7 +6953,7 @@ HypreAME::GetEigenvector(unsigned int i) const
    return *eigenvectors[i];
 }
 
-HypreParVector **
+Vector **
 HypreAME::StealEigenvectors()
 {
    if ( eigenvectors == NULL )
@@ -6944,11 +6962,11 @@ HypreAME::StealEigenvectors()
    }
 
    // Set the local pointers to NULL so that they won't be deleted later
-   HypreParVector ** vecs = eigenvectors;
+   Vector ** vecs = (Vector**)eigenvectors;
    eigenvectors = NULL;
    multi_vec = NULL;
 
-   return vecs;
+   return (Vector**)vecs;
 }
 
 }

linalg/hypre.hpp

@@ -18,7 +18,9 @@
 
 #include "../general/globals.hpp"
 #include "sparsemat.hpp"
+#include "eigensolvers.hpp"
 #include "hypre_parcsr.hpp"
+#include "eigensolvers.hpp"
 #include <mpi.h>
 
 // Enable internal hypre timing routines
@@ -2146,7 +2148,7 @@ public:
     A. Knyazev, M. Argentati, I. Lashuk, and E. Ovtchinnikov, SISC, 29(5),
     2224-2239, 2007.
 */
-class HypreLOBPCG
+class HypreLOBPCG : public SymGenEigensolver
 {
 private:
    MPI_Comm comm;
@@ -2236,38 +2238,43 @@ public:
    HypreLOBPCG(MPI_Comm comm);
    ~HypreLOBPCG();
 
-   void SetTol(real_t tol);
+   void SetTol(real_t tol) override;
    // not implemented in HYPRE
    // real_t GetTol() const;
    void SetRelTol(real_t rel_tol);
    // not implemented in HYPRE
    // real_t GetRelTol() const;
-   void SetMaxIter(int max_iter);
+   void SetMaxIter(int max_iter) override;
    // not implemented in HYPRE
    // int GetMaxIter() const;
-   void SetPrintLevel(int logging);
-   void SetNumModes(int num_eigs) { nev = num_eigs; }
+   void SetPrintLevel(int logging) override;
+   void SetNumModes(int num_eigs) override { nev = num_eigs; }
    void SetPrecondUsageMode(int pcg_mode);
    void SetRandomSeed(int s) { seed = s; }
    void SetInitialVectors(int num_vecs, HypreParVector ** vecs);
 
    // The following four methods support general operators
    void SetPreconditioner(Solver & precond);
-   void SetOperator(Operator & A);
-   void SetMassMatrix(Operator & M);
+   void SetOperators(const Operator & A, const Operator & B) override
+   { SetOperator(A); SetMassMatrix(B); }
+   void SetOperator(const Operator & A);
+   void SetMassMatrix(const Operator & M);
    void SetSubSpaceProjector(Operator & proj) { subSpaceProj = &proj; }
 
    /// Solve the eigenproblem
-   void Solve();
+   void Solve() override;
+
+   int GetNumConverged() const override { return nev; }
 
    /// Collect the converged eigenvalues
-   void GetEigenvalues(Array<real_t> & eigenvalues) const;
+   void GetEigenvalues(Array<real_t> & eigenvalues) const override;
 
    /// Extract a single eigenvector
-   const HypreParVector & GetEigenvector(unsigned int i) const;
+   const Vector & GetEigenvector(unsigned int i) const override;
 
    /// Transfer ownership of the converged eigenvectors
-   HypreParVector ** StealEigenvectors() { return multi_vec->StealVectors(); }
+   Vector ** StealEigenvectors() override
+   { return (Vector**)multi_vec->StealVectors(); }
 };
 
 /** AME eigenvalue solver in hypre
@@ -2292,7 +2299,7 @@ public:
     mass matrix but it seems unlikely that this would be useful so it is not the
     default behavior.
 */
-class HypreAME
+class HypreAME : public SymGenEigensolver
 {
 private:
    int myid;
@@ -2321,28 +2328,31 @@ public:
    HypreAME(MPI_Comm comm);
    ~HypreAME();
 
-   void SetTol(real_t tol);
+   void SetTol(real_t tol) override;
    void SetRelTol(real_t rel_tol);
-   void SetMaxIter(int max_iter);
-   void SetPrintLevel(int logging);
-   void SetNumModes(int num_eigs);
+   void SetMaxIter(int max_iter) override;
+   void SetPrintLevel(int logging) override;
+   void SetNumModes(int num_eigs) override;
 
    // The following four methods support operators of type HypreParMatrix.
    void SetPreconditioner(HypreSolver & precond);
+   void SetOperators(const Operator & opA, const Operator & opB) override;
    void SetOperator(const HypreParMatrix & A);
    void SetMassMatrix(const HypreParMatrix & M);
 
    /// Solve the eigenproblem
-   void Solve();
+   void Solve() override;
+
+   int GetNumConverged() const override { return nev; }
 
    /// Collect the converged eigenvalues
-   void GetEigenvalues(Array<real_t> & eigenvalues) const;
+   void GetEigenvalues(Array<real_t> & eigenvalues) const override;
 
    /// Extract a single eigenvector
-   const HypreParVector & GetEigenvector(unsigned int i) const;
+   const Vector & GetEigenvector(unsigned int i) const override;
 
    /// Transfer ownership of the converged eigenvectors
-   HypreParVector ** StealEigenvectors();
+   Vector ** StealEigenvectors() override;
 };
 
 }

linalg/linalg.hpp

@@ -28,6 +28,7 @@
 #include "symmat.hpp"
 #include "ode.hpp"
 #include "solvers.hpp"
+#include "eigensolvers.hpp"
 #include "handle.hpp"
 #include "invariants.hpp"
 #include "constraints.hpp"
@@ -57,6 +58,10 @@
 #include "ginkgo.hpp"
 #endif
 
+#ifdef MFEM_USE_ARPACK
+#include "arpack.hpp"
+#endif
+
 #ifdef MFEM_USE_MKL_PARDISO
 #include "pardiso.hpp"
 #endif

linalg/operator.hpp

@@ -844,6 +844,22 @@ public:
 };
 
 
+/// Zero Operator N: x -> 0.
+class ZeroOperator : public Operator
+{
+public:
+   /// Create an zero operator of size @a n.
+   explicit ZeroOperator(int n) : Operator(n) { }
+
+   /// Operator application
+   void Mult(const Vector &x, Vector &y) const override
+   { y.SetSize(width); y = 0_r; }
+
+   /// Application of the transpose
+   void MultTranspose(const Vector &x, Vector &y) const override
+   { y.SetSize(width); y = 0_r; }
+};
+
 /// Identity Operator I: x -> x.
 class IdentityOperator : public Operator
 {

makefile

@@ -123,20 +123,15 @@ EXAMPLE_SUBDIRS = amgx caliper ginkgo hiop petsc pumi sundials superlu moonolith
 EXAMPLE_DIRS := examples $(addprefix examples/,$(EXAMPLE_SUBDIRS))
 EXAMPLE_TEST_DIRS := examples
 
-MINIAPP_ALL_SUBDIRS = common electromagnetics meshing performance tools \
+MINIAPP_SUBDIRS = common electromagnetics meshing performance tools \
  toys nurbs gslib adjoint solvers shifted mtop parelag tribol autodiff dfem \
  hooke multidomain dpg hdiv-linear-solver spde diag-smoothers contact \
  fluids/navier fluids/schrodinger-flow plasma plasma/pic
-MINIAPP_RECURSIVE_SUBDIRS = plasma/pic
-MINIAPP_SUBDIRS := $(filter-out \
-   $(MINIAPP_RECURSIVE_SUBDIRS),$(MINIAPP_ALL_SUBDIRS))
-MINIAPP_ALL_DIRS := $(addprefix miniapps/,$(MINIAPP_ALL_SUBDIRS))
 MINIAPP_DIRS := $(addprefix miniapps/,$(MINIAPP_SUBDIRS))
 MINIAPP_TEST_DIRS := $(filter-out %/common,$(MINIAPP_DIRS))
 MINIAPP_USE_COMMON := $(addprefix miniapps/,electromagnetics meshing tools \
- toys gslib shifted dpg diag-smoothers fluids/navier plasma plasma/pic)
+ toys shifted dpg diag-smoothers fluids/navier plasma plasma/pic)
 
-EM_ALL_DIRS = $(EXAMPLE_DIRS) $(MINIAPP_ALL_DIRS)
 EM_DIRS = $(EXAMPLE_DIRS) $(MINIAPP_DIRS)
 
 TEST_SUBDIRS = unit
@@ -151,7 +146,7 @@ MFEM_BUILD_DIR ?= .
 BUILD_DIR := $(MFEM_BUILD_DIR)
 BUILD_REAL_DIR := $(abspath $(BUILD_DIR))
 ifneq ($(BUILD_REAL_DIR),$(MFEM_REAL_DIR))
-   BUILD_SUBDIRS = $(DIRS) config $(EM_ALL_DIRS) doc $(TEST_DIRS)
+   BUILD_SUBDIRS = $(DIRS) config $(EM_DIRS) doc $(TEST_DIRS)
    CONFIG_FILE_DEF = -DMFEM_CONFIG_FILE='"$(BUILD_REAL_DIR)/config/_config.hpp"'
    BLD := $(if $(BUILD_REAL_DIR:$(CURDIR)=),$(BUILD_DIR)/,)
    $(if $(word 2,$(BLD)),$(error Spaces in BLD = "$(BLD)" are not supported))
@@ -307,7 +302,7 @@ endif
 MFEM_REQ_LIB_DEPS = SUPERLU MUMPS METIS FMS CONDUIT SIDRE LAPACK SUNDIALS\
  SUITESPARSE STRUMPACK GINKGO GNUTLS HDF5 NETCDF SLEPC PETSC MPFR PUMI HIOP\
  GSLIB OCCA CEED RAJA UMPIRE MKL_CPARDISO MKL_PARDISO AMGX MAGMA CALIPER PARELAG\
- TRIBOL BENCHMARK MOONOLITH ALGOIM
+ TRIBOL BENCHMARK MOONOLITH ALGOIM ARPACK
 
 
 PETSC_ERROR_MSG = $(if $(PETSC_FOUND),,. PETSC config not found: $(PETSC_VARS))
@@ -376,7 +371,8 @@ MFEM_DEFINES = MFEM_VERSION MFEM_VERSION_STRING MFEM_GIT_STRING MFEM_USE_MPI\
  MFEM_USE_SIMD MFEM_USE_ADIOS2 MFEM_USE_MKL_CPARDISO MFEM_USE_MKL_PARDISO MFEM_USE_AMGX\
  MFEM_USE_MAGMA MFEM_USE_MUMPS MFEM_USE_ADFORWARD MFEM_USE_CODIPACK MFEM_USE_CALIPER\
  MFEM_USE_BENCHMARK MFEM_USE_PARELAG MFEM_USE_TRIBOL MFEM_USE_ALGOIM MFEM_USE_ENZYME\
- MFEM_SOURCE_DIR MFEM_INSTALL_DIR MFEM_SHARED_BUILD MFEM_USE_DOUBLE MFEM_USE_SINGLE
+ MFEM_SOURCE_DIR MFEM_INSTALL_DIR MFEM_SHARED_BUILD MFEM_USE_DOUBLE MFEM_USE_SINGLE\
+MFEM_USE_ARPACK
 
 # List of makefile variables that will be written to config.mk:
 MFEM_CONFIG_VARS = MFEM_CXX MFEM_HOST_CXX MFEM_CPPFLAGS MFEM_CXXFLAGS\
@@ -488,10 +484,10 @@ $(OBJECT_FILES): $(BLD)%.o: $(SRC)%.cpp $(CONFIG_MK)
 
 all: examples miniapps $(TEST_DIRS)
 
-.PHONY: miniapps $(EM_ALL_DIRS) $(TEST_DIRS)
+.PHONY: miniapps $(EM_DIRS) $(TEST_DIRS)
 miniapps: $(MINIAPP_DIRS)
 $(MINIAPP_USE_COMMON): miniapps/common
-$(EM_ALL_DIRS) $(TEST_DIRS): lib
+$(EM_DIRS) $(TEST_DIRS): lib
 	$(MAKE) -C $(BLD)$(@)
 
 .PHONY: doc
@@ -699,7 +695,7 @@ local-config:
 .PHONY: build-config
 build-config:
 	for d in $(BUILD_SUBDIRS); do mkdir -p $(BLD)$${d}; done
-	for dir in "" $(addsuffix /,config $(EM_ALL_DIRS) doc $(TEST_DIRS)); do\
+	for dir in "" $(addsuffix /,config $(EM_DIRS) doc $(TEST_DIRS)); do \
 	   printf "# Auto-generated file.\n%s\n%s\n" \
 	      "MFEM_DIR = $(MFEM_REAL_DIR)" \
 	      "include \$$(MFEM_DIR)/$${dir}makefile" \
@@ -738,6 +734,7 @@ status info:
 	$(info MFEM_TIMER_TYPE        = $(MFEM_TIMER_TYPE))
 	$(info MFEM_USE_SUNDIALS      = $(MFEM_USE_SUNDIALS))
 	$(info MFEM_USE_SUITESPARSE   = $(MFEM_USE_SUITESPARSE))
+	$(info MFEM_USE_ARPACK        = $(MFEM_USE_ARPACK))
 	$(info MFEM_USE_SUPERLU       = $(MFEM_USE_SUPERLU))
 	$(info MFEM_USE_SUPERLU5      = $(MFEM_USE_SUPERLU5))
 	$(info MFEM_USE_MUMPS         = $(MFEM_USE_MUMPS))
@@ -801,15 +798,13 @@ status info:
 
 ASTYLE = $(ASTYLE_BIN) --options=$(SRC)config/mfem.astylerc
 ASTYLE_VER = "Artistic Style Version 3.1"
-FORMAT_FILES = $(foreach dir,$(DIRS) $(EM_ALL_DIRS) config,$(dir)/*.?pp)
+FORMAT_FILES = $(foreach dir,$(DIRS) $(EM_DIRS) config,$(dir)/*.?pp)
 TESTS_SUBDIRS = unit benchmarks convergence mem_manager par-mesh-format
-UNIT_TESTS_SUBDIRS = general linalg mesh fem miniapps ceed enzyme dfem
-MINIAPPS_SUBDIRS = dpg/util hooke/operators hooke/preconditioners \
-   hooke/materials hooke/kernels
+UNIT_TESTS_SUBDIRS = general linalg mesh fem miniapps ceed enzyme
+MINIAPPS_SUBDIRS = dpg/util hooke/operators hooke/preconditioners hooke/materials hooke/kernels
 FORMAT_FILES += $(foreach dir,$(TESTS_SUBDIRS),tests/$(dir)/*.?pp)
 FORMAT_FILES += $(foreach dir,$(UNIT_TESTS_SUBDIRS),tests/unit/$(dir)/*.?pp)
 FORMAT_FILES += $(foreach dir,$(MINIAPPS_SUBDIRS),miniapps/$(dir)/*.?pp)
-FORMAT_FILES += config/cmake/config.hpp.in config/config.hpp.in mfem*.hpp
 FORMAT_EXCLUDE = general/tinyxml2.cpp tests/unit/catch.hpp
 FORMAT_LIST = $(filter-out $(FORMAT_EXCLUDE),$(wildcard $(FORMAT_FILES)))
 
@@ -840,29 +835,14 @@ mfem_check_command = \
 # Verify the C++ code styling in MFEM and check that std::cout and std::cerr are
 # not used in the library (use mfem::out and mfem::err instead).
 style:
+	@echo "Applying C++ code style..."
 	@astyle_version="$$($(ASTYLE_BIN) --version)";\
 	 if [ "$$astyle_version" != $(ASTYLE_VER) ]; then\
 	    printf "%s\n" "Invalid astyle version: '$$astyle_version'"\
 	           "Please use: '"$(ASTYLE_VER)"'";\
 	    exit 1;\
 	 fi
-	@err_code=0; \
-	if command -v git 2>&1 > /dev/null && [ -d $(MFEM_DIR)/.git ]; then \
-	   echo "Checking if all git files are selected for formatting ..."; \
-	   ls -1 $(FORMAT_FILES) | sort > format-files-make.txt; \
-	   git -C $(MFEM_DIR) ls-files '*.[ch]pp*' | sort \
-	      > format-files-git.txt; \
-	   cat format-files-make.txt format-files-git.txt | sort | uniq \
-	      > format-files-make-plus-git.txt; \
-	   rm -f format-files-git.txt; \
-	   $(call mfem_check_command,\
-	      diff format-files-make.txt format-files-make-plus-git.txt | \
-	      grep "^> ",\
-	      "All git files are selected for formatting",\
-	      "The above git files are NOT selected for formatting"); \
-	   rm -f format-files-make.txt format-files-make-plus-git.txt; \
-	fi; \
-	echo "Applying C++ code style...";\
+	@err_code=0;\
 	$(call mfem_check_command,\
 	    $(ASTYLE) $(FORMAT_LIST) | grep Formatted,\
 	    "No source files were changed",\

miniapps/gslib/CMakeLists.txt

@@ -46,16 +46,8 @@ if (MFEM_USE_GSLIB)
 
   if (MFEM_ENABLE_TESTING)
     foreach (test "schwarz_ex1" "field-diff" "findpts" "field-interp")
-      if (MFEM_USE_MPI)
-        add_test(NAME ${test}
-          COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 1
-          ${MPIEXEC_PREFLAGS}
-          $<TARGET_FILE:${test}> -no-vis
-          ${MPIEXEC_POSTFLAGS})
-      else()
-        add_test(NAME ${test}
-          COMMAND $<TARGET_FILE:${test}> -no-vis)
-      endif()
+      add_test(NAME ${test}
+        COMMAND $<TARGET_FILE:${test}> -no-vis)
     endforeach()
   endif()
 

miniapps/gslib/makefile

@@ -83,16 +83,11 @@ include $(MFEM_TEST_MK)
 
 # Testing: Parallel vs. serial runs
 RUN_MPI = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) $(MFEM_MPI_NP)
-ifeq ($(MFEM_USE_MPI),YES)
-   RUN_MPI_SER = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) 1
-else
-   RUN_MPI_SER =
-endif
 TEST_NAME := GSLIB miniapp
 %-test-par: %
 	@$(call mfem-test,$<, $(RUN_MPI), $(TEST_NAME))
 %-test-seq: %
-	@$(call mfem-test,$<, $(RUN_MPI_SER), $(TEST_NAME))
+	@$(call mfem-test,$<,, $(TEST_NAME))
 
 # Testing: "test" target and mfem-test* variables are defined in config/test.mk
 

miniapps/nurbs/nurbs_ex11p.cpp

@@ -329,8 +329,8 @@ int main(int argc, char *argv[])
 
       for (int i=0; i<nev; i++)
       {
-         // convert eigenvector from HypreParVector to ParGridFunction
-         x = lobpcg->GetEigenvector(i);
+         // convert eigenvector from Vector to ParGridFunction
+         x.Distribute(lobpcg->GetEigenvector(i));
 
          mode_name << "mode_" << setfill('0') << setw(2) << i << "."
                    << setfill('0') << setw(6) << myid;
@@ -357,8 +357,8 @@ int main(int argc, char *argv[])
                  << ", Lambda = " << eigenvalues[i] << endl;
          }
 
-         // convert eigenvector from HypreParVector to ParGridFunction
-         x = lobpcg->GetEigenvector(i);
+         // convert eigenvector from Vector to ParGridFunction
+         x.Distribute(lobpcg->GetEigenvector(i));
 
          mode_sock << "parallel " << num_procs << " " << myid << "\n"
                    << "solution\n" << *pmesh << x << flush

miniapps/plasma/makefile

@@ -28,12 +28,6 @@ endif
 
 PLASMA_SUBDIRS = pic
 
-SUBDIRS_ALL = $(addsuffix /all,$(PLASMA_SUBDIRS))
-SUBDIRS_TEST = $(addsuffix /test,$(PLASMA_SUBDIRS))
-SUBDIRS_TEST_NOCLEAN = $(addsuffix /test-noclean,$(PLASMA_SUBDIRS))
-SUBDIRS_CLEAN = $(addsuffix /clean,$(PLASMA_SUBDIRS))
-SUBDIRS_TPRINT = $(addsuffix /test-print,$(PLASMA_SUBDIRS))
-
 .SUFFIXES:
 .SUFFIXES: .o .cpp .mk
 .PHONY: all lib-common clean clean-build clean-exec
@@ -43,24 +37,31 @@ COMMON_LIB = -L$(MFEM_BUILD_DIR)/miniapps/common -lmfem-common
 
 # If MFEM_SHARED is set, add the ../common rpath
 COMMON_LIB += $(if $(MFEM_SHARED:YES=),,\
-   $(MFEM_XLINKER)-rpath,$(abspath $(MFEM_BUILD_DIR)/miniapps/common))
+   $(if $(MFEM_USE_CUDA:YES=),$(CXX_XLINKER),$(CUDA_XLINKER))-rpath,$(abspath\
+   $(MFEM_BUILD_DIR)/miniapps/common))
+
+COMMON_O=
 
 # Remove built-in rules
 %: %.cpp
 %.o: %.cpp
 
-all: $(MINIAPPS) $(SUBDIRS_ALL) 
+all: $(MINIAPPS) subdirs
 
-.PHONY: $(SUBDIRS_ALL) $(SUBDIRS_TEST) $(SUBDIRS_TEST_NOCLEAN) \
-   $(SUBDIRS_CLEAN) $(SUBDIRS_TPRINT)
-$(SUBDIRS_ALL) $(SUBDIRS_TEST) $(SUBDIRS_TEST_NOCLEAN) $(SUBDIRS_CLEAN):
-	$(MAKE) -C $(@D) $(@F)
-$(SUBDIRS_TPRINT):
-	@$(MAKE) -C $(@D) $(@F)
+.PHONY: subdirs $(PLASMA_SUBDIRS)
+subdirs: $(PLASMA_SUBDIRS)
+$(PLASMA_SUBDIRS): lib-common
+	$(MAKE) -C $(BLD)$(@)
 
 # Rules for building the miniapps
-%: $(SRC)%.cpp $(MFEM_LIB_FILE) $(CONFIG_MK) | lib-common
-	$(MFEM_CXX) $(MFEM_LINK_FLAGS) $< -o $@ $(COMMON_LIB) $(MFEM_LIBS)
+%: $(SRC)%.cpp $(COMMON_O) $(MFEM_LIB_FILE) $(CONFIG_MK) | lib-common
+	$(MFEM_CXX) $(MFEM_LINK_FLAGS) $< -o $@ $(COMMON_O) $(COMMON_LIB) \
+        $(MFEM_LIBS)
+
+# Rules for compiling miniapp dependencies
+$(COMMON_O) $(addsuffix _solver.o,$(MINIAPPS)): \
+%.o: $(SRC)%.cpp $(SRC)%.hpp $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) -c $(<) -o $(@)
 
 # Rule for building lib-common
 lib-common:
@@ -68,9 +69,6 @@ lib-common:
 
 MFEM_TESTS = MINIAPPS
 include $(MFEM_TEST_MK)
-test: $(SUBDIRS_TEST)
-test-noclean: $(SUBDIRS_TEST_NOCLEAN)
-test-print: $(SUBDIRS_TPRINT)
 
 # Testing: Specific execution options
 RUN_MPI = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) $(MFEM_MPI_NP)
@@ -81,9 +79,14 @@ RUN_MPI = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) $(MFEM_MPI_NP)
 $(MFEM_LIB_FILE):
 	$(error The MFEM library is not built)
 
-clean: clean-build clean-exec $(SUBDIRS_CLEAN)
+ALL_CLEAN_SUBDIRS = $(addsuffix /clean,$(PLASMA_SUBDIRS))
+.PHONY: $(ALL_CLEAN_SUBDIRS)
+$(ALL_CLEAN_SUBDIRS):
+	$(MAKE) -C $(BLD)$(@D) $(@F)
 
-clean-build:
+clean: clean-build clean-exec
+
+clean-build: $(addsuffix /clean,$(PLASMA_SUBDIRS))
 	rm -f *.o *~ $(SEQ_MINIAPPS) $(PAR_MINIAPPS)
 	rm -rf *.dSYM *.TVD.*breakpoints
 

miniapps/plasma/pic/makefile

@@ -20,7 +20,6 @@ CONFIG_MK = $(or $(wildcard $(MFEM_BUILD_DIR)/config/config.mk),\
 MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 
-SEQ_MINIAPPS =
 PAR_MINIAPPS =
 
 ifeq ($(MFEM_USE_GSLIB),YES)
@@ -28,9 +27,9 @@ ifeq ($(MFEM_USE_GSLIB),YES)
 endif
 
 ifeq ($(MFEM_USE_MPI),NO)
-   MINIAPPS = $(SEQ_MINIAPPS)
+   MINIAPPS =
 else
-   MINIAPPS = $(PAR_MINIAPPS) $(SEQ_MINIAPPS)
+   MINIAPPS = $(PAR_MINIAPPS)
 endif
 
 .SUFFIXES:
@@ -51,7 +50,7 @@ COMMON_LIB += $(if $(MFEM_SHARED:YES=),,\
 all: $(MINIAPPS)
 
 # Rules for building the miniapps
-%: $(SRC)%.cpp $(MFEM_LIB_FILE) $(CONFIG_MK) | lib-common
+electrostatic-pic: electrostatic-pic.cpp $(MFEM_LIB_FILE) $(CONFIG_MK) | lib-common
 	$(MFEM_CXX) $(MFEM_FLAGS) -c $<
 	$(MFEM_CXX) $(MFEM_LINK_FLAGS) -o $@ $@.o $(COMMON_LIB) $(MFEM_LIBS)
 

miniapps/tools/CMakeLists.txt

@@ -61,10 +61,6 @@ if (MFEM_USE_MPI)
       LIBRARIES mfem-common)
       add_dependencies(gridfunction-bounds copy_miniapps_tools_data)
 
-    add_mfem_miniapp(random-gridfunction-bounds
-      MAIN random-gridfunction-bounds.cpp
-      LIBRARIES mfem)
-
     add_mfem_miniapp(plor-transfer
       MAIN plor-transfer.cpp LIBRARIES mfem)
 

miniapps/tools/makefile

@@ -23,8 +23,7 @@ MFEM_LIB_FILE = mfem_is_not_built
 SEQ_MINIAPPS = display-basis load-dc convert-dc get-values lor-transfer \
 	tmop-check-metric tmop-metric-magnitude compare-dc
 
-PAR_MINIAPPS = nodal-transfer plor-transfer gridfunction-bounds \
-	random-gridfunction-bounds
+PAR_MINIAPPS = nodal-transfer plor-transfer gridfunction-bounds
 
 ifeq ($(MFEM_USE_MPI),NO)
    MINIAPPS = $(SEQ_MINIAPPS)
@@ -80,7 +79,7 @@ RUN_MPI = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) $(MFEM_MPI_NP)
 # Do not test: display-basis, load-dc, convert-dc, get-values, lor-transfer, plor-transfer
 NO_TEST_APPS = display-basis load-dc convert-dc get-values lor-transfer \
 	plor-transfer tmop-check-metric tmop-metric-magnitude gridfunction-bounds \
-	random-gridfunction-bounds compare-dc
+	compare-dc
 $(foreach app,$(NO_TEST_APPS),$(app)-test-seq $(app)-test-par):
 	@true
 

miniapps/tools/random-gridfunction-bounds.cpp

@@ -1,263 +0,0 @@
-// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
-// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
-// LICENSE and NOTICE for details. LLNL-CODE-806117.
-//
-// This file is part of the MFEM library. For more information and source code
-// availability visit https://mfem.org.
-//
-// MFEM is free software; you can redistribute it and/or modify it under the
-// terms of the BSD-3 license. We welcome feedback and contributions, see file
-// CONTRIBUTING.md for details.
-//
-//    ---------------------------------------------------------------------
-//     Compute bounds of a random grid function on a generated tensor mesh
-//    ---------------------------------------------------------------------
-//
-// This miniapp generates a 1D segment mesh or 2D quad mesh, builds a random
-// discontinuous grid function, computes element-wise piecewise linear bounds,
-// and visualizes the input field together with the lower and upper bounds.
-//
-// Compile with: make random-gridfunction-bounds
-//
-// Sample runs:
-//   mpirun -np 4 random-gridfunction-bounds
-//   mpirun -np 4 random-gridfunction-bounds -nx 64 -o 6 -ref 3 -d hip
-
-#include "mfem.hpp"
-
-#include <algorithm>
-#include <type_traits>
-
-using namespace mfem;
-using namespace std;
-
-void VisualizeField(ParMesh &pmesh, ParGridFunction &input,
-                    char *title, int pos_x, int pos_y);
-
-int main(int argc, char *argv[])
-{
-   Mpi::Init(argc, argv);
-   Hypre::Init();
-
-   int dim = 2;
-   int nx = 16;
-   int order = 4;
-   int num_comp = 2;
-   int ref = 2;
-   int niter = 1000;
-   int seed = 12345;
-   bool kernel_only = true;
-   bool visualization = false;
-   const char *device_config = "cpu";
-
-   OptionsParser args(argc, argv);
-   args.AddOption(&dim, "-dim", "--dimension",
-                  "Dimension of the generated tensor-product mesh (1 or 2).");
-   args.AddOption(&nx, "-nx", "--num-elements",
-                  "Number of elements in each mesh direction.");
-   args.AddOption(&order, "-o", "--order",
-                  "Polynomial degree of the random discontinuous field.");
-   args.AddOption(&num_comp, "-nc", "--num-components",
-                  "Number of vector components in the ParFiniteElementSpace.");
-   args.AddOption(&ref, "-ref", "--piecewise-linear-ref-factor",
-                  "Scaling factor for the resolution of the piecewise linear "
-                  "bounds. If less than 2, the resolution is picked "
-                  "automatically.");
-   args.AddOption(&niter, "-ni", "--num-iters",
-                  "Number of times to evaluate the bounds.");
-   args.AddOption(&seed, "-rs", "--random-seed",
-                  "Random seed used to initialize the field.");
-   args.AddOption(&device_config, "-d", "--device",
-                  "Device configuration string, see Device::Configure().");
-   args.AddOption(&kernel_only, "-ko", "--kernel-only",
-                  "-no-ko", "--no-kernel-only",
-                  "Run only PLBound::GetElementBoundsKernel on a prebuilt "
-                  "element E-vector.");
-   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
-                  "--no-visualization",
-                  "Enable or disable GLVis visualization.");
-   args.ParseCheck();
-
-   MFEM_VERIFY(dim == 1 || dim == 2, "dim must be 1 or 2.");
-   MFEM_VERIFY(nx > 0, "nx must be positive.");
-   MFEM_VERIFY(order >= 0, "order must be non-negative.");
-   MFEM_VERIFY(num_comp > 0, "num_comp must be positive.");
-   MFEM_VERIFY(niter > 0, "niter must be positive.");
-
-   Device device(device_config);
-   if (Mpi::Root()) { device.Print(); }
-
-   Mesh mesh = (dim == 1) ?
-               Mesh::MakeCartesian1D(nx, 1.0) :
-               Mesh::MakeCartesian2D(nx, nx, Element::QUADRILATERAL, true,
-                                     1.0, 1.0);
-   ParMesh pmesh(MPI_COMM_WORLD, mesh);
-
-   const int mesh_dim = pmesh.Dimension();
-   L2_FECollection fec(order, mesh_dim, BasisType::GaussLobatto);
-   ParFiniteElementSpace fes(&pmesh, &fec, num_comp, Ordering::byNODES);
-   ParGridFunction input(&fes);
-   input.Randomize(seed + Mpi::WorldRank());
-   input.UseDevice(true);
-
-   L2_FECollection fec_pc(0, mesh_dim);
-   ParFiniteElementSpace fes_pc(&pmesh, &fec_pc, num_comp, Ordering::byNODES);
-   ParGridFunction lowerb(&fes_pc), upperb(&fes_pc);
-   Vector lower_vec, upper_vec;
-
-   PLBound plb(&fes, ref*(fes.GetMaxElementOrder() + 1));
-   if (kernel_only)
-   {
-      const FiniteElement &fe = *fes.GetTypicalFE();
-      const int rdim = fe.GetDim();
-      const int nd = fe.GetDof();
-      const int fes_dim = fes.GetVDim();
-      Vector e_vec(nd*fes_dim*fes.GetNE(), Device::GetDeviceMemoryType());
-      e_vec.UseDevice(true);
-
-      const ElementRestrictionOperator *elem_restr =
-         fes.GetElementRestriction(ElementDofOrdering::LEXICOGRAPHIC);
-      MFEM_VERIFY(elem_restr != nullptr,
-                  "Element restriction is required for kernel-only mode.");
-      elem_restr->Mult(input, e_vec);
-
-      for (int i = 0; i < niter; i++)
-      {
-         plb.GetElementBoundsKernel(rdim, fes_dim, e_vec, lower_vec, upper_vec);
-      }
-   }
-   else
-   {
-      for (int i = 0; i < niter; i++)
-      {
-         input.GetElementBounds(plb, lower_vec, upper_vec);
-      }
-   }
-
-   const real_t *lower_data = lower_vec.HostRead();
-   const real_t *upper_data = upper_vec.HostRead();
-
-   // Build a host reference from the lexicographic E-vector and the scalar
-   // PLBound::GetNDBounds path to avoid re-entering the device dispatch.
-   const bool use_dev = input.UseDevice();
-   PLBound plb_host(&fes, ref*(fes.GetMaxElementOrder() + 1));
-   Vector lower_ref, upper_ref;
-   const FiniteElement &fe = *fes.GetTypicalFE();
-   const int rdim = fe.GetDim();
-   const int nd = fe.GetDof();
-   const int nel = fes.GetNE();
-   const int fes_dim = fes.GetVDim();
-   Vector e_vec_ref(nd*fes_dim*nel);
-   lower_ref.SetSize(nel*fes_dim);
-   upper_ref.SetSize(nel*fes_dim);
-   const ElementRestrictionOperator *elem_restr =
-      fes.GetElementRestriction(ElementDofOrdering::LEXICOGRAPHIC);
-   MFEM_VERIFY(elem_restr != nullptr,
-               "Element restriction is required for host reference.");
-   input.UseDevice(false);
-   input.HostRead();
-   elem_restr->Mult(input, e_vec_ref);
-   input.UseDevice(use_dev);
-   const real_t *e_ref_data = e_vec_ref.HostRead();
-
-   for (int d = 0; d < fes_dim; d++)
-   {
-      for (int e = 0; e < nel; e++)
-      {
-         Vector coeff(nd);
-         for (int i = 0; i < nd; i++)
-         {
-            coeff(i) = e_ref_data[i + nd*(d + fes_dim*e)];
-         }
-         Vector lower_c, upper_c;
-         plb_host.GetNDBounds(rdim, coeff, lower_c, upper_c);
-         lower_ref(e + d*nel) = lower_c.Min();
-         upper_ref(e + d*nel) = upper_c.Max();
-      }
-   }
-   const real_t *lower_ref_data = lower_ref.HostRead();
-   const real_t *upper_ref_data = upper_ref.HostRead();
-
-   MFEM_VERIFY(lower_vec.Size() == lower_ref.Size() &&
-               upper_vec.Size() == upper_ref.Size(),
-               "Reference element-bound vectors have inconsistent sizes.");
-
-   real_t lower_diff = 0.0;
-   real_t upper_diff = 0.0;
-   for (int i = 0; i < lower_vec.Size(); i++)
-   {
-      lower_diff = std::max(lower_diff,
-                            std::abs(lower_data[i] - lower_ref_data[i]));
-   }
-   for (int i = 0; i < upper_vec.Size(); i++)
-   {
-      upper_diff = std::max(upper_diff,
-                            std::abs(upper_data[i] - upper_ref_data[i]));
-   }
-   MPI_Allreduce(MPI_IN_PLACE, &lower_diff, 1, MPITypeMap<real_t>::mpi_type,
-                 MPI_MAX, pmesh.GetComm());
-   MPI_Allreduce(MPI_IN_PLACE, &upper_diff, 1, MPITypeMap<real_t>::mpi_type,
-                 MPI_MAX, pmesh.GetComm());
-
-   const real_t verify_tol = std::is_same<real_t, float>::value ?
-                             real_t(1.0e-5) : real_t(1.0e-12);
-   MFEM_VERIFY(lower_diff <= verify_tol && upper_diff <= verify_tol,
-               "Device element bounds do not match host reference.");
-
-   lowerb = lower_vec;
-   upperb = upper_vec;
-
-   real_t lower_min = lowerb.Min();
-   real_t upper_max = upperb.Max();
-   MPI_Allreduce(MPI_IN_PLACE, &lower_min, 1, MPITypeMap<real_t>::mpi_type,
-                 MPI_MIN, pmesh.GetComm());
-   MPI_Allreduce(MPI_IN_PLACE, &upper_max, 1, MPITypeMap<real_t>::mpi_type,
-                 MPI_MAX, pmesh.GetComm());
-
-   if (Mpi::Root())
-   {
-      cout << "dim: " << mesh_dim << '\n'
-           << "nx: " << nx << '\n'
-           << "order: " << order << '\n'
-           << "num components: " << num_comp << '\n'
-           << "PL bound control-point factor: " << ref << '\n'
-           << "iterations: " << niter << '\n'
-           << "kernel-only mode: " << (kernel_only ? "yes" : "no") << '\n'
-           << "host/device lower max diff: " << lower_diff << '\n'
-           << "host/device upper max diff: " << upper_diff << '\n'
-           << "global lower bound minimum: " << lower_min << '\n'
-           << "global upper bound maximum: " << upper_max << endl;
-   }
-
-   if (visualization)
-   {
-      char title1[] = "Random input gridfunction";
-      char title2[] = "Element-wise lower bound";
-      char title3[] = "Element-wise upper bound";
-      VisualizeField(pmesh, input, title1, 0, 0);
-      VisualizeField(pmesh, lowerb, title2, 450, 0);
-      VisualizeField(pmesh, upperb, title3, 900, 0);
-   }
-
-   return 0;
-}
-
-void VisualizeField(ParMesh &pmesh, ParGridFunction &input,
-                    char *title, int pos_x, int pos_y)
-{
-   socketstream sock;
-   if (pmesh.GetMyRank() == 0)
-   {
-      sock.open("localhost", 19916);
-      sock << "solution\n";
-   }
-   pmesh.PrintAsOne(sock);
-   input.SaveAsOne(sock);
-   if (pmesh.GetMyRank() == 0)
-   {
-      sock << "window_title '" << title << "'\n"
-           << "window_geometry "
-           << pos_x << " " << pos_y << " " << 400 << " " << 400 << "\n"
-           << "keys jRmclApppppppppppp//]]]]]]]]" << endl;
-   }
-}