Merge branch 'master' into eigensolvers-dev

examples/ex11p.cpp

@@ -321,8 +321,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;
@@ -350,8 +350,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/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/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);
 

linalg/CMakeLists.txt

@@ -55,6 +55,7 @@ list(APPEND HDRS
   dinvariants.hpp
   dtensor.hpp
   dual.hpp
+  eigensolvers.hpp
   filteredsolver.hpp
   handle.hpp
   invariants.hpp

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,7 +6962,7 @@ 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;
 

linalg/hypre.hpp

@@ -18,6 +18,7 @@
 
 #include "../general/globals.hpp"
 #include "sparsemat.hpp"
+#include "eigensolvers.hpp"
 #include "hypre_parcsr.hpp"
 #include <mpi.h>
 
@@ -2146,7 +2147,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 +2237,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 +2298,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 +2327,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"

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
 {

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