From c5f3ca14f0b5be84f624a69c604cd9f7ef394cb4 Mon Sep 17 00:00:00 2001
From: yianzeng <57109024+yianzeng@users.noreply.github.com>
Date: Tue, 13 Jun 2023 16:40:39 +0100
Subject: [PATCH 01/22] Create eigensolver.py

---
 firedrake/eigensolver.py | 117 +++++++++++++++++++++++++++++++++++++++
 1 file changed, 117 insertions(+)
 create mode 100644 firedrake/eigensolver.py

diff --git a/firedrake/eigensolver.py b/firedrake/eigensolver.py
new file mode 100644
index 0000000000..1f19ff97b3
--- /dev/null
+++ b/firedrake/eigensolver.py
@@ -0,0 +1,117 @@
+from firedrake.assemble import assemble
+from firedrake.function import Function
+from firedrake import solving_utils
+from firedrake import utils
+from firedrake.petsc import PETSc, OptionsManager, flatten_parameters
+from firedrake.logging import warning
+from firedrake.exceptions import ConvergenceError
+try:
+    from slepc4py import SLEPc
+except ImportError:
+    import sys
+    warning("Unable to import SLEPc, eigenvalue computation not possible (try firedrake-update --slepc)")
+    sys.exit(0)
+__all__ = ["LinearEigenproblem",
+           "LinearEigensolver"]
+class LinearEigenproblem():
+    r"""Linear eigenvalue problem A(u; v) = lambda * M (u; v)."""
+    def __init__(self, A, M=None, bcs=None):
+        r"""
+        :param A: the bilinear form A(u, v)
+        :param M: the mass form M(u, v) (optional)
+        :param bcs: the boundary conditions (optional)
+        """
+        self.A = A  # LHS 
+        args = A.arguments()
+        v, u = args[0], args[1]
+        if M:
+            self.M = M
+        else:
+            from ufl import inner, dx
+            self.M = inner(u, v) * dx
+        self.output_space = u.function_space()
+        self.bcs = bcs
+
+    def dirichlet_bcs(self):  # cargo cult
+        r"""Return an iterator over the Dirichlet boundary conditions in self.bcs."""
+        for bc in self.bcs:
+            yield from bc.dirichlet_bcs()
+
+    @utils.cached_property
+    def dm(self):  # cargo cult
+        r"""Return the function space's distributed mesh associated with self.output_space (a cached property)."""
+        return self.output_space.dm
+
+class LinearEigensolver(OptionsManager): 
+    r"""Solve a :class:`LinearEigenproblem`."""
+    # DEFAULT_SNES_PARAMETERS = {"snes_type": "newtonls",
+    #                            "snes_linesearch_type": "basic"}
+
+    # # Looser default tolerance for KSP inside SNES.
+    DEFAULT_EPS_PARAMETERS = solving_utils.DEFAULT_KSP_PARAMETERS.copy() 
+
+    DEFAULT_EPS_PARAMETERS = {"eps_gen_non_hermitian": None, 
+                              "st_pc_factor_shift_type": "NONZERO",
+                              "eps_type": "krylovschur",
+                              "eps_largest_imaginary": None,
+                              "eps_tol":1e-10}  # to change
+
+    def __init__(self, problem, n_evals, *, options_prefix=None, solver_parameters=None):
+        r'''
+        :param problem: :class:`LinearEigenproblem` to solve.
+        :param n_evals: number of eigenvalues to compute.
+        :param options_prefix: options prefix to use for the eigensolver.
+        :param solver_parameters: PETSc options for the eigenvalue problem.
+        '''
+
+        self.es = SLEPc.EPS().create(comm=problem.dm.comm)
+        self._problem = problem 
+        self.n_evals = n_evals
+        solver_parameters = flatten_parameters(solver_parameters or {})
+        for key in self.DEFAULT_EPS_PARAMETERS:
+            value = self.DEFAULT_EPS_PARAMETERS[key]
+            solver_parameters.setdefault(key, value)
+        super().__init__(solver_parameters, options_prefix)
+        self.set_from_options(self.es)
+    
+    def check_es_convergence(self):
+        r"""Check the convergence of the eigenvalue problem."""
+        r = self.es.getConvergedReason()
+        try:
+            reason = SLEPc.EPS.ConvergedReasons[r]
+        except KeyError:
+            reason = "unknown reason (petsc4py enum incomplete?), try with -eps_converged_reason"
+        if r < 0:
+            raise ConvergenceError(r"""Eigenproblem failed to converge after %d iterations.
+        Reason:
+        %s""" % (self.es.getIterationNumber(), reason))
+
+    def solve(self):
+        r"""Solve the eigenproblem, return the number of converged eigenvalues."""
+        if self._problem.bcs is None: # Neumann BCs 
+            self.A_mat = assemble(self._problem.A, bcs=self._problem.bcs).M.handle
+            self.M_mat = assemble(self._problem.M, bcs=self._problem.bcs).M.handle 
+        else:  # Dirichlet BCs - in this case we are solving Mu = lambda Au, so the eigenvalue needs to be inverted
+            self.A_mat = assemble(self._problem.M, bcs=self._problem.bcs, weight=0.).M.handle
+            self.M_mat = assemble(self._problem.A, bcs=self._problem.bcs, weight=1.).M.handle
+        self.es.setOperators(self.A_mat, self.M_mat)
+        self.es.setDimensions(nev=self.n_evals,ncv=2*self.n_evals)  # SLEPc recommended params
+        with self.inserted_options():
+            self.es.solve()
+        nconv = self.es.getConverged()
+        if nconv == 0:
+            raise ConvergenceError("Did not converge any eigenvalues.") 
+        return nconv
+
+    def eigenvalue(self, i):
+        r"""Return the `i`th eigenvalue of the solved problem."""
+        return self.es.getEigenvalue(i)
+
+    def eigenfunction(self, i):
+        r"""Return the `i`th eigenfunction of the solved problem."""
+        eigenmodes_real = Function(self._problem.output_space)  # fn of V
+        eigenmodes_imag = Function(self._problem.output_space)
+        with eigenmodes_real.dat.vec_wo as vr:
+            with eigenmodes_imag.dat.vec_wo as vi:
+                self.es.getEigenvector(i, vr, vi)  # gets the i-th eigenvector
+        return eigenmodes_real, eigenmodes_imag  # returns Firedrake fns

From bb838628a12c556a05ebadc84ddb9ea610869e3c Mon Sep 17 00:00:00 2001
From: yianzeng <57109024+yianzeng@users.noreply.github.com>
Date: Tue, 13 Jun 2023 16:42:04 +0100
Subject: [PATCH 02/22] Create test_eigensolver.py

---
 tests/regression/test_eigensolver.py | 103 +++++++++++++++++++++++++++
 1 file changed, 103 insertions(+)
 create mode 100644 tests/regression/test_eigensolver.py

diff --git a/tests/regression/test_eigensolver.py b/tests/regression/test_eigensolver.py
new file mode 100644
index 0000000000..bf50ee26f0
--- /dev/null
+++ b/tests/regression/test_eigensolver.py
@@ -0,0 +1,103 @@
+# tests/regression/test_helmholtz.py
+
+from os.path import abspath, dirname, join
+import numpy as np
+import pytest
+
+from firedrake import *
+
+cwd = abspath(dirname(__file__))
+
+
+def poisson_1d(n, quadrilateral=False, degree=1, mesh=None):
+    # Create mesh and define function space
+    if mesh is None:
+        mesh = IntervalMesh(n, 0, pi)
+    V = FunctionSpace(mesh, "CG", degree)
+    # Define variational problem
+    u = TrialFunction(V)
+    v = TestFunction(V)
+    a = (inner(grad(u), grad(v))) * dx
+
+    # Create eigenproblem with boundary conditions
+    bc = DirichletBC(V, 0.0, "on_boundary")
+    eigenprob = LinearEigenproblem(a, bcs=bc)
+    
+    # Create corresponding eigensolver, looking for n eigenvalues
+    eigensolver = LinearEigensolver(eigenprob, n)
+    ncov = eigensolver.solve()
+    print(ncov)
+
+    # boffi solns 
+    h = pi /n
+    true_values = np.zeros(ncov-2)
+    estimates = np.zeros(ncov-2)
+    for k in range(ncov-2):
+        true_val = 6 / h**2
+        true_val *= (1-cos((k+1)*h))/(2+cos((k+1)*h)) # k+1 because we skip the trivial 0 eigenvalue
+        true_values[k] = true_val
+
+        estimates[k] = 1/eigensolver.eigenvalue(k).real # takes real part 
+    return true_values, estimates
+
+@pytest.mark.parametrize(('n', 'quadrilateral', 'degree', 'tolerance'),
+                         [(5, False, 1, 1e-14),
+                          (10, False, 1, 1e-14),
+                          (20, False, 1, 1e-14),
+                          (30, False, 1, 1e-14)])
+def test_poisson_eigenvalue_convergence(n, quadrilateral, degree, tolerance):
+    true_values, estimates = poisson_1d(n, quadrilateral=quadrilateral, degree=degree)
+    assert np.allclose(true_values, estimates, rtol=tolerance)
+
+
+# tests/regression/test_helmholtz.py
+
+from os.path import abspath, dirname, join
+import numpy as np
+import pytest
+
+from firedrake import *
+
+cwd = abspath(dirname(__file__))
+
+
+def poisson_2d(n, quadrilateral=False, degree=1, mesh=None):
+    # Create mesh and define function space
+    if mesh is None:
+        mesh = UnitSquareMesh(n, 0, pi)
+    V = FunctionSpace(mesh, "CG", degree)
+    # Define variational problem
+    u = TrialFunction(V)
+    v = TestFunction(V)
+    a = (inner(grad(u), grad(v))) * dx
+
+    # Create eigenproblem with boundary conditions
+    bc = DirichletBC(V, 0.0, "on_boundary")
+    eigenprob = LinearEigenproblem(a, bcs=bc)
+    
+    # Create corresponding eigensolver, looking for n eigenvalues
+    eigensolver = LinearEigensolver(eigenprob, n)
+    ncov = eigensolver.solve()
+    print(ncov)
+
+    # boffi solns 
+    h = pi /n
+    true_values = np.zeros(ncov-2)
+    estimates = np.zeros(ncov-2)
+    for k in range(ncov-2):
+        true_val = 6 / h**2
+        true_val *= (1-cos((k+1)*h))/(2+cos((k+1)*h)) # k+1 because we skip the trivial 0 eigenvalue
+        true_values[k] = true_val
+        estimates[k] = 1/eigensolver.eigenvalue(k).real # takes real part 
+    return true_values, estimates
+
+@pytest.mark.parametrize(('n', 'quadrilateral', 'degree', 'tolerance'),
+                         [(5, False, 1, 1e-14),
+                          (10, False, 1, 1e-14),
+                          (20, False, 1, 1e-14),
+                          (30, False, 1, 1e-14)])
+def test_poisson_eigenvalue_convergence(n, quadrilateral, degree, tolerance):
+    true_values, estimates = poisson_1d(n, quadrilateral=quadrilateral, degree=degree)
+    assert np.allclose(true_values, estimates, rtol=tolerance)
+
+

From 6441d55e1b743526c88237a4213c449e38b2ea54 Mon Sep 17 00:00:00 2001
From: yianzeng <57109024+yianzeng@users.noreply.github.com>
Date: Tue, 13 Jun 2023 16:43:11 +0100
Subject: [PATCH 03/22] Update test_eigensolver.py

---
 tests/regression/test_eigensolver.py | 2 --
 1 file changed, 2 deletions(-)

diff --git a/tests/regression/test_eigensolver.py b/tests/regression/test_eigensolver.py
index bf50ee26f0..24aef7f992 100644
--- a/tests/regression/test_eigensolver.py
+++ b/tests/regression/test_eigensolver.py
@@ -99,5 +99,3 @@ def poisson_2d(n, quadrilateral=False, degree=1, mesh=None):
 def test_poisson_eigenvalue_convergence(n, quadrilateral, degree, tolerance):
     true_values, estimates = poisson_1d(n, quadrilateral=quadrilateral, degree=degree)
     assert np.allclose(true_values, estimates, rtol=tolerance)
-
-

From ffb849d278a1428dd8dcadf872144fc063a5eae3 Mon Sep 17 00:00:00 2001
From: yianzeng <57109024+yianzeng@users.noreply.github.com>
Date: Wed, 14 Jun 2023 14:09:08 +0100
Subject: [PATCH 04/22] Update test_eigensolver.py

Added docstrings, added a line to make sure compared values are sorted
---
 tests/regression/test_eigensolver.py | 72 +++++-----------------------
 1 file changed, 11 insertions(+), 61 deletions(-)

diff --git a/tests/regression/test_eigensolver.py b/tests/regression/test_eigensolver.py
index 24aef7f992..d3e3527d04 100644
--- a/tests/regression/test_eigensolver.py
+++ b/tests/regression/test_eigensolver.py
@@ -8,8 +8,10 @@
 
 cwd = abspath(dirname(__file__))
 
-
-def poisson_1d(n, quadrilateral=False, degree=1, mesh=None):
+def evals(n, quadrilateral=False, degree=1, mesh=None):
+    '''We base this test on the 1D Poisson problem with Dirichlet boundary conditions,
+    outlined in part 1 of Daniele Boffi's 'Finite element approximation
+    of eigenvalue problems' Acta Numerica 2010'''
     # Create mesh and define function space
     if mesh is None:
         mesh = IntervalMesh(n, 0, pi)
@@ -26,7 +28,6 @@ def poisson_1d(n, quadrilateral=False, degree=1, mesh=None):
     # Create corresponding eigensolver, looking for n eigenvalues
     eigensolver = LinearEigensolver(eigenprob, n)
     ncov = eigensolver.solve()
-    print(ncov)
 
     # boffi solns 
     h = pi /n
@@ -38,64 +39,13 @@ def poisson_1d(n, quadrilateral=False, degree=1, mesh=None):
         true_values[k] = true_val
 
         estimates[k] = 1/eigensolver.eigenvalue(k).real # takes real part 
-    return true_values, estimates
-
-@pytest.mark.parametrize(('n', 'quadrilateral', 'degree', 'tolerance'),
-                         [(5, False, 1, 1e-14),
-                          (10, False, 1, 1e-14),
-                          (20, False, 1, 1e-14),
-                          (30, False, 1, 1e-14)])
-def test_poisson_eigenvalue_convergence(n, quadrilateral, degree, tolerance):
-    true_values, estimates = poisson_1d(n, quadrilateral=quadrilateral, degree=degree)
-    assert np.allclose(true_values, estimates, rtol=tolerance)
-
-
-# tests/regression/test_helmholtz.py
-
-from os.path import abspath, dirname, join
-import numpy as np
-import pytest
-
-from firedrake import *
-
-cwd = abspath(dirname(__file__))
-
-
-def poisson_2d(n, quadrilateral=False, degree=1, mesh=None):
-    # Create mesh and define function space
-    if mesh is None:
-        mesh = UnitSquareMesh(n, 0, pi)
-    V = FunctionSpace(mesh, "CG", degree)
-    # Define variational problem
-    u = TrialFunction(V)
-    v = TestFunction(V)
-    a = (inner(grad(u), grad(v))) * dx
-
-    # Create eigenproblem with boundary conditions
-    bc = DirichletBC(V, 0.0, "on_boundary")
-    eigenprob = LinearEigenproblem(a, bcs=bc)
-    
-    # Create corresponding eigensolver, looking for n eigenvalues
-    eigensolver = LinearEigensolver(eigenprob, n)
-    ncov = eigensolver.solve()
-    print(ncov)
-
-    # boffi solns 
-    h = pi /n
-    true_values = np.zeros(ncov-2)
-    estimates = np.zeros(ncov-2)
-    for k in range(ncov-2):
-        true_val = 6 / h**2
-        true_val *= (1-cos((k+1)*h))/(2+cos((k+1)*h)) # k+1 because we skip the trivial 0 eigenvalue
-        true_values[k] = true_val
-        estimates[k] = 1/eigensolver.eigenvalue(k).real # takes real part 
-    return true_values, estimates
+    return sorted(true_values), sorted(estimates) # sort in case order of evals returned differently
 
 @pytest.mark.parametrize(('n', 'quadrilateral', 'degree', 'tolerance'),
-                         [(5, False, 1, 1e-14),
-                          (10, False, 1, 1e-14),
-                          (20, False, 1, 1e-14),
-                          (30, False, 1, 1e-14)])
-def test_poisson_eigenvalue_convergence(n, quadrilateral, degree, tolerance):
-    true_values, estimates = poisson_1d(n, quadrilateral=quadrilateral, degree=degree)
+                         [(5, False, 1, 1e-13),
+                          (10, False, 1, 1e-13),
+                          (20, False, 1, 1e-13),
+                          (30, False, 1, 1e-13)])
+def test_evals_1d(n, quadrilateral, degree, tolerance):
+    true_values, estimates = evals(n, quadrilateral=quadrilateral, degree=degree)
     assert np.allclose(true_values, estimates, rtol=tolerance)

From 984e5f62e94fb8743eb9d1ae2b395fa6a4ac9c27 Mon Sep 17 00:00:00 2001
From: David Ham <David.Ham@imperial.ac.uk>
Date: Thu, 15 Jun 2023 18:15:03 +0100
Subject: [PATCH 05/22] lint and other cleanup

---
 firedrake/eigensolver.py             | 49 ++++++++++++++--------------
 tests/regression/test_eigensolver.py | 27 +++++++--------
 2 files changed, 39 insertions(+), 37 deletions(-)

diff --git a/firedrake/eigensolver.py b/firedrake/eigensolver.py
index 1f19ff97b3..495c044cad 100644
--- a/firedrake/eigensolver.py
+++ b/firedrake/eigensolver.py
@@ -1,18 +1,19 @@
 from firedrake.assemble import assemble
 from firedrake.function import Function
-from firedrake import solving_utils
 from firedrake import utils
-from firedrake.petsc import PETSc, OptionsManager, flatten_parameters
-from firedrake.logging import warning
+from firedrake.petsc import OptionsManager, flatten_parameters
 from firedrake.exceptions import ConvergenceError
 try:
     from slepc4py import SLEPc
 except ImportError:
-    import sys
-    warning("Unable to import SLEPc, eigenvalue computation not possible (try firedrake-update --slepc)")
-    sys.exit(0)
+    raise ImportError(
+        "Unable to import SLEPc, eigenvalue computation not possible "
+        "(try firedrake-update --slepc)"
+    )
 __all__ = ["LinearEigenproblem",
            "LinearEigensolver"]
+
+
 class LinearEigenproblem():
     r"""Linear eigenvalue problem A(u; v) = lambda * M (u; v)."""
     def __init__(self, A, M=None, bcs=None):
@@ -21,7 +22,7 @@ def __init__(self, A, M=None, bcs=None):
         :param M: the mass form M(u, v) (optional)
         :param bcs: the boundary conditions (optional)
         """
-        self.A = A  # LHS 
+        self.A = A  # LHS
         args = A.arguments()
         v, u = args[0], args[1]
         if M:
@@ -42,19 +43,15 @@ def dm(self):  # cargo cult
         r"""Return the function space's distributed mesh associated with self.output_space (a cached property)."""
         return self.output_space.dm
 
-class LinearEigensolver(OptionsManager): 
-    r"""Solve a :class:`LinearEigenproblem`."""
-    # DEFAULT_SNES_PARAMETERS = {"snes_type": "newtonls",
-    #                            "snes_linesearch_type": "basic"}
 
-    # # Looser default tolerance for KSP inside SNES.
-    DEFAULT_EPS_PARAMETERS = solving_utils.DEFAULT_KSP_PARAMETERS.copy() 
+class LinearEigensolver(OptionsManager):
+    r"""Solve a :class:`LinearEigenproblem`."""
 
-    DEFAULT_EPS_PARAMETERS = {"eps_gen_non_hermitian": None, 
+    DEFAULT_EPS_PARAMETERS = {"eps_gen_non_hermitian": None,
                               "st_pc_factor_shift_type": "NONZERO",
                               "eps_type": "krylovschur",
                               "eps_largest_imaginary": None,
-                              "eps_tol":1e-10}  # to change
+                              "eps_tol": 1e-10}
 
     def __init__(self, problem, n_evals, *, options_prefix=None, solver_parameters=None):
         r'''
@@ -65,7 +62,7 @@ def __init__(self, problem, n_evals, *, options_prefix=None, solver_parameters=N
         '''
 
         self.es = SLEPc.EPS().create(comm=problem.dm.comm)
-        self._problem = problem 
+        self._problem = problem
         self.n_evals = n_evals
         solver_parameters = flatten_parameters(solver_parameters or {})
         for key in self.DEFAULT_EPS_PARAMETERS:
@@ -73,34 +70,38 @@ def __init__(self, problem, n_evals, *, options_prefix=None, solver_parameters=N
             solver_parameters.setdefault(key, value)
         super().__init__(solver_parameters, options_prefix)
         self.set_from_options(self.es)
-    
+
     def check_es_convergence(self):
         r"""Check the convergence of the eigenvalue problem."""
         r = self.es.getConvergedReason()
         try:
             reason = SLEPc.EPS.ConvergedReasons[r]
         except KeyError:
-            reason = "unknown reason (petsc4py enum incomplete?), try with -eps_converged_reason"
+            reason = ("unknown reason (petsc4py enum incomplete?), "
+                      "try with -eps_converged_reason")
         if r < 0:
-            raise ConvergenceError(r"""Eigenproblem failed to converge after %d iterations.
+            raise ConvergenceError(
+                r"""Eigenproblem failed to converge after %d iterations.
         Reason:
-        %s""" % (self.es.getIterationNumber(), reason))
+        %s""" % (self.es.getIterationNumber(), reason)
+            )
 
     def solve(self):
         r"""Solve the eigenproblem, return the number of converged eigenvalues."""
-        if self._problem.bcs is None: # Neumann BCs 
+        if self._problem.bcs is None:  # Neumann BCs
             self.A_mat = assemble(self._problem.A, bcs=self._problem.bcs).M.handle
-            self.M_mat = assemble(self._problem.M, bcs=self._problem.bcs).M.handle 
+            self.M_mat = assemble(self._problem.M, bcs=self._problem.bcs).M.handle
         else:  # Dirichlet BCs - in this case we are solving Mu = lambda Au, so the eigenvalue needs to be inverted
             self.A_mat = assemble(self._problem.M, bcs=self._problem.bcs, weight=0.).M.handle
             self.M_mat = assemble(self._problem.A, bcs=self._problem.bcs, weight=1.).M.handle
         self.es.setOperators(self.A_mat, self.M_mat)
-        self.es.setDimensions(nev=self.n_evals,ncv=2*self.n_evals)  # SLEPc recommended params
+        # SLEPc recommended params
+        self.es.setDimensions(nev=self.n_evals, ncv=2*self.n_evals)
         with self.inserted_options():
             self.es.solve()
         nconv = self.es.getConverged()
         if nconv == 0:
-            raise ConvergenceError("Did not converge any eigenvalues.") 
+            raise ConvergenceError("Did not converge any eigenvalues.")
         return nconv
 
     def eigenvalue(self, i):
diff --git a/tests/regression/test_eigensolver.py b/tests/regression/test_eigensolver.py
index d3e3527d04..f5f1461059 100644
--- a/tests/regression/test_eigensolver.py
+++ b/tests/regression/test_eigensolver.py
@@ -1,17 +1,15 @@
 # tests/regression/test_helmholtz.py
 
-from os.path import abspath, dirname, join
 import numpy as np
 import pytest
 
 from firedrake import *
 
-cwd = abspath(dirname(__file__))
 
-def evals(n, quadrilateral=False, degree=1, mesh=None):
-    '''We base this test on the 1D Poisson problem with Dirichlet boundary conditions,
-    outlined in part 1 of Daniele Boffi's 'Finite element approximation
-    of eigenvalue problems' Acta Numerica 2010'''
+def evals(n, degree=1, mesh=None):
+    '''We base this test on the 1D Poisson problem with Dirichlet boundary
+    conditions, outlined in part 1 of Daniele Boffi's
+    'Finite element approximation of eigenvalue problems' Acta Numerica 2010'''
     # Create mesh and define function space
     if mesh is None:
         mesh = IntervalMesh(n, 0, pi)
@@ -24,22 +22,25 @@ def evals(n, quadrilateral=False, degree=1, mesh=None):
     # Create eigenproblem with boundary conditions
     bc = DirichletBC(V, 0.0, "on_boundary")
     eigenprob = LinearEigenproblem(a, bcs=bc)
-    
+
     # Create corresponding eigensolver, looking for n eigenvalues
     eigensolver = LinearEigensolver(eigenprob, n)
     ncov = eigensolver.solve()
 
-    # boffi solns 
-    h = pi /n
+    # boffi solns
+    h = pi / n
     true_values = np.zeros(ncov-2)
     estimates = np.zeros(ncov-2)
     for k in range(ncov-2):
         true_val = 6 / h**2
-        true_val *= (1-cos((k+1)*h))/(2+cos((k+1)*h)) # k+1 because we skip the trivial 0 eigenvalue
+        # k+1 because we skip the trivial 0 eigenvalue
+        true_val *= (1-cos((k+1)*h))/(2+cos((k+1)*h))
         true_values[k] = true_val
 
-        estimates[k] = 1/eigensolver.eigenvalue(k).real # takes real part 
-    return sorted(true_values), sorted(estimates) # sort in case order of evals returned differently
+        estimates[k] = 1/eigensolver.eigenvalue(k).real
+    # sort in case order of numerical and analytic values differs.
+    return sorted(true_values), sorted(estimates)
+
 
 @pytest.mark.parametrize(('n', 'quadrilateral', 'degree', 'tolerance'),
                          [(5, False, 1, 1e-13),
@@ -47,5 +48,5 @@ def evals(n, quadrilateral=False, degree=1, mesh=None):
                           (20, False, 1, 1e-13),
                           (30, False, 1, 1e-13)])
 def test_evals_1d(n, quadrilateral, degree, tolerance):
-    true_values, estimates = evals(n, quadrilateral=quadrilateral, degree=degree)
+    true_values, estimates = evals(n, degree=degree)
     assert np.allclose(true_values, estimates, rtol=tolerance)

From 474486d09a37ad4cedbd034d45caf373ea3913b4 Mon Sep 17 00:00:00 2001
From: David Ham <David.Ham@imperial.ac.uk>
Date: Thu, 15 Jun 2023 18:16:35 +0100
Subject: [PATCH 06/22] remove unused parameter

---
 tests/regression/test_eigensolver.py | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/tests/regression/test_eigensolver.py b/tests/regression/test_eigensolver.py
index f5f1461059..e285cf70e6 100644
--- a/tests/regression/test_eigensolver.py
+++ b/tests/regression/test_eigensolver.py
@@ -42,11 +42,11 @@ def evals(n, degree=1, mesh=None):
     return sorted(true_values), sorted(estimates)
 
 
-@pytest.mark.parametrize(('n', 'quadrilateral', 'degree', 'tolerance'),
-                         [(5, False, 1, 1e-13),
-                          (10, False, 1, 1e-13),
-                          (20, False, 1, 1e-13),
-                          (30, False, 1, 1e-13)])
-def test_evals_1d(n, quadrilateral, degree, tolerance):
+@pytest.mark.parametrize(('n', 'degree', 'tolerance'),
+                         [(5, 1, 1e-13),
+                          (10, 1, 1e-13),
+                          (20, 1, 1e-13),
+                          (30, 1, 1e-13)])
+def test_evals_1d(n, degree, tolerance):
     true_values, estimates = evals(n, degree=degree)
     assert np.allclose(true_values, estimates, rtol=tolerance)

From 94f2e3e338059187d9f37d9f9224f1fe149c6357 Mon Sep 17 00:00:00 2001
From: David Ham <David.Ham@imperial.ac.uk>
Date: Fri, 16 Jun 2023 12:54:53 +0100
Subject: [PATCH 07/22] further cleanup and draft new BC approach

---
 firedrake/__init__.py                |   1 +
 firedrake/eigensolver.py             | 104 +++++++++++++++++++--------
 tests/regression/test_eigensolver.py |   3 -
 3 files changed, 74 insertions(+), 34 deletions(-)

diff --git a/firedrake/__init__.py b/firedrake/__init__.py
index 13a1a46f49..369fb85274 100644
--- a/firedrake/__init__.py
+++ b/firedrake/__init__.py
@@ -98,6 +98,7 @@
 from firedrake.ufl_expr import *
 from firedrake.utility_meshes import *
 from firedrake.variational_solver import *
+from firedrake.eigensolver import *
 from firedrake.vector import *
 from firedrake.version import __version__ as ver, __version_info__, check  # noqa: F401
 from firedrake.ensemble import *
diff --git a/firedrake/eigensolver.py b/firedrake/eigensolver.py
index 495c044cad..a18d5f9978 100644
--- a/firedrake/eigensolver.py
+++ b/firedrake/eigensolver.py
@@ -1,3 +1,4 @@
+"""Specify and solve finite element eigenproblems."""
 from firedrake.assemble import assemble
 from firedrake.function import Function
 from firedrake import utils
@@ -6,22 +7,45 @@
 try:
     from slepc4py import SLEPc
 except ImportError:
-    raise ImportError(
-        "Unable to import SLEPc, eigenvalue computation not possible "
-        "(try firedrake-update --slepc)"
-    )
+    SLEPc = None
 __all__ = ["LinearEigenproblem",
            "LinearEigensolver"]
 
 
 class LinearEigenproblem():
-    r"""Linear eigenvalue problem A(u; v) = lambda * M (u; v)."""
-    def __init__(self, A, M=None, bcs=None):
-        r"""
-        :param A: the bilinear form A(u, v)
-        :param M: the mass form M(u, v) (optional)
-        :param bcs: the boundary conditions (optional)
-        """
+    """Linear eigenvalue problem
+
+    The problem has the form::
+
+        A(u, v) = λ * M(u, v)
+
+    Parameters
+    ----------
+    A : ufl.Form
+        the bilinear form A(u, v)
+    M : ufl.Form
+        the mass form M(u, v),  defaults to u * v * dx.
+    bcs : DirichletBC or list of DirichletBC
+        the boundary conditions
+    bc_shift: float
+        the value to shift the boundary condition eigenvalues by.
+
+    Notes
+    -----
+
+    If Dirichlet boundary conditions are supplied then these will result in the
+    eigenproblem having a nullspace spanned by the basis functions with support
+    on the boundary. To facilitate solution, this is shifted by the specified
+    amount. It is the user's responsibility to ensure that the shift is not
+    close to an actual eigenvalue of the system.
+    """
+    def __init__(self, A, M=None, bcs=None, bc_shift=666.0):
+        if not SLEPc:
+            raise ImportError(
+                "Unable to import SLEPc, eigenvalue computation not possible "
+                "(try firedrake-update --slepc)"
+            )
+
         self.A = A  # LHS
         args = A.arguments()
         v, u = args[0], args[1]
@@ -32,20 +56,33 @@ def __init__(self, A, M=None, bcs=None):
             self.M = inner(u, v) * dx
         self.output_space = u.function_space()
         self.bcs = bcs
+        self.bc_shift = bc_shift
 
-    def dirichlet_bcs(self):  # cargo cult
-        r"""Return an iterator over the Dirichlet boundary conditions in self.bcs."""
+    def dirichlet_bcs(self):
+        """Return an iterator over the Dirichlet boundary conditions."""
         for bc in self.bcs:
             yield from bc.dirichlet_bcs()
 
     @utils.cached_property
     def dm(self):  # cargo cult
-        r"""Return the function space's distributed mesh associated with self.output_space (a cached property)."""
+        r"""Return the dm associated with the output space."""
         return self.output_space.dm
 
 
 class LinearEigensolver(OptionsManager):
-    r"""Solve a :class:`LinearEigenproblem`."""
+    r"""Solve a LinearEigenproblem.
+
+    Parameters
+    ----------
+    problem : LinearEigenproblem
+        The eigenproblem to solve.
+    n_evals : int
+        The number of eigenvalues to compute.
+    options_prefix : str
+        The options prefix to use for the eigensolver.
+    solver_parameters : dict
+        PETSc options for the eigenvalue problem.
+    """
 
     DEFAULT_EPS_PARAMETERS = {"eps_gen_non_hermitian": None,
                               "st_pc_factor_shift_type": "NONZERO",
@@ -53,13 +90,8 @@ class LinearEigensolver(OptionsManager):
                               "eps_largest_imaginary": None,
                               "eps_tol": 1e-10}
 
-    def __init__(self, problem, n_evals, *, options_prefix=None, solver_parameters=None):
-        r'''
-        :param problem: :class:`LinearEigenproblem` to solve.
-        :param n_evals: number of eigenvalues to compute.
-        :param options_prefix: options prefix to use for the eigensolver.
-        :param solver_parameters: PETSc options for the eigenvalue problem.
-        '''
+    def __init__(self, problem, n_evals, *, options_prefix=None,
+                 solver_parameters=None):
 
         self.es = SLEPc.EPS().create(comm=problem.dm.comm)
         self._problem = problem
@@ -87,13 +119,17 @@ def check_es_convergence(self):
             )
 
     def solve(self):
-        r"""Solve the eigenproblem, return the number of converged eigenvalues."""
-        if self._problem.bcs is None:  # Neumann BCs
-            self.A_mat = assemble(self._problem.A, bcs=self._problem.bcs).M.handle
-            self.M_mat = assemble(self._problem.M, bcs=self._problem.bcs).M.handle
-        else:  # Dirichlet BCs - in this case we are solving Mu = lambda Au, so the eigenvalue needs to be inverted
-            self.A_mat = assemble(self._problem.M, bcs=self._problem.bcs, weight=0.).M.handle
-            self.M_mat = assemble(self._problem.A, bcs=self._problem.bcs, weight=1.).M.handle
+        r"""Solve the eigenproblem.
+
+        Returns
+        -------
+        int
+            The number of Eigenvalues found.
+        """
+        self.A_mat = assemble(self._problem.A, bcs=self._problem.bcs).M.handle
+        self.M_mat = assemble(self._problem.M, bcs=self._problem.bcs,
+                              weight=self._problem.bc_shift).M.handle
+
         self.es.setOperators(self.A_mat, self.M_mat)
         # SLEPc recommended params
         self.es.setDimensions(nev=self.n_evals, ncv=2*self.n_evals)
@@ -105,11 +141,17 @@ def solve(self):
         return nconv
 
     def eigenvalue(self, i):
-        r"""Return the `i`th eigenvalue of the solved problem."""
+        r"""Return the i-th eigenvalue of the solved problem."""
         return self.es.getEigenvalue(i)
 
     def eigenfunction(self, i):
-        r"""Return the `i`th eigenfunction of the solved problem."""
+        r"""Return the i-th eigenfunction of the solved problem.
+
+        Returns
+        -------
+        (Function, Function)
+            The real and imaginary parts of the eigenfunction.
+        """
         eigenmodes_real = Function(self._problem.output_space)  # fn of V
         eigenmodes_imag = Function(self._problem.output_space)
         with eigenmodes_real.dat.vec_wo as vr:
diff --git a/tests/regression/test_eigensolver.py b/tests/regression/test_eigensolver.py
index e285cf70e6..1585cb7681 100644
--- a/tests/regression/test_eigensolver.py
+++ b/tests/regression/test_eigensolver.py
@@ -1,8 +1,5 @@
-# tests/regression/test_helmholtz.py
-
 import numpy as np
 import pytest
-
 from firedrake import *
 
 

From d82298419665e4c6557ab50a93838ebc9ac13ff0 Mon Sep 17 00:00:00 2001
From: David Ham <David.Ham@imperial.ac.uk>
Date: Fri, 23 Jun 2023 14:01:20 +0100
Subject: [PATCH 08/22] Seemingly working Dirichlet BCs

---
 firedrake/eigensolver.py             |  6 ++++--
 tests/regression/test_eigensolver.py | 11 +++++++----
 2 files changed, 11 insertions(+), 6 deletions(-)

diff --git a/firedrake/eigensolver.py b/firedrake/eigensolver.py
index a18d5f9978..8edcdc00e7 100644
--- a/firedrake/eigensolver.py
+++ b/firedrake/eigensolver.py
@@ -127,8 +127,10 @@ def solve(self):
             The number of Eigenvalues found.
         """
         self.A_mat = assemble(self._problem.A, bcs=self._problem.bcs).M.handle
-        self.M_mat = assemble(self._problem.M, bcs=self._problem.bcs,
-                              weight=self._problem.bc_shift).M.handle
+        self.M_mat = assemble(
+            self._problem.M, bcs=self._problem.bcs,
+            weight=self._problem.bc_shift and 1./self._problem.bc_shift
+        ).M.handle
 
         self.es.setOperators(self.A_mat, self.M_mat)
         # SLEPc recommended params
diff --git a/tests/regression/test_eigensolver.py b/tests/regression/test_eigensolver.py
index 1585cb7681..2ebc8cbea6 100644
--- a/tests/regression/test_eigensolver.py
+++ b/tests/regression/test_eigensolver.py
@@ -26,15 +26,18 @@ def evals(n, degree=1, mesh=None):
 
     # boffi solns
     h = pi / n
-    true_values = np.zeros(ncov-2)
-    estimates = np.zeros(ncov-2)
-    for k in range(ncov-2):
+    true_values = np.zeros(ncov-1)
+    estimates = np.zeros(ncov-1)
+    for k in range(ncov-1):
         true_val = 6 / h**2
         # k+1 because we skip the trivial 0 eigenvalue
         true_val *= (1-cos((k+1)*h))/(2+cos((k+1)*h))
         true_values[k] = true_val
 
-        estimates[k] = 1/eigensolver.eigenvalue(k).real
+        estimates[k] = eigensolver.eigenvalue(k).real
+
+    true_values[-1] = eigenprob.bc_shift
+
     # sort in case order of numerical and analytic values differs.
     return sorted(true_values), sorted(estimates)
 

From e29c8d546f80da1238a2713498e4d145fcbf076c Mon Sep 17 00:00:00 2001
From: David Ham <David.Ham@imperial.ac.uk>
Date: Fri, 23 Jun 2023 14:15:32 +0100
Subject: [PATCH 09/22] tell LaTeX about unicode lambda

---
 docs/source/conf.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/docs/source/conf.py b/docs/source/conf.py
index 47f140d76e..8670269803 100644
--- a/docs/source/conf.py
+++ b/docs/source/conf.py
@@ -262,6 +262,7 @@
 \sphinxDUC{2218}{$\circ$}
 \sphinxDUC{22C5}{$\cdot$}
 \sphinxDUC{25A3}{$\boxdot$}
+\sphinxDUC{03BB}{$\lambda$}
 % Sphinx equivalent of
 % \DeclareUnicodeCharacter{}{}
 

From 59af019b163edc155e82c618c1f93918e69fb228 Mon Sep 17 00:00:00 2001
From: David Ham <David.Ham@imperial.ac.uk>
Date: Fri, 23 Jun 2023 15:56:16 +0100
Subject: [PATCH 10/22] Expose slepc parameters

---
 firedrake/eigensolver.py | 13 ++++++++++---
 1 file changed, 10 insertions(+), 3 deletions(-)

diff --git a/firedrake/eigensolver.py b/firedrake/eigensolver.py
index 8edcdc00e7..523532158d 100644
--- a/firedrake/eigensolver.py
+++ b/firedrake/eigensolver.py
@@ -82,6 +82,12 @@ class LinearEigensolver(OptionsManager):
         The options prefix to use for the eigensolver.
     solver_parameters : dict
         PETSc options for the eigenvalue problem.
+    ncv : int
+        Maximum dimension of the subspace to be used by the solver. See
+        `SLEPc.EPS.setDimensions`.
+    mpd : int
+        Maximum dimension allowed for the projected problem. See
+        `SLEPc.EPS.setDimensions`.
     """
 
     DEFAULT_EPS_PARAMETERS = {"eps_gen_non_hermitian": None,
@@ -91,11 +97,13 @@ class LinearEigensolver(OptionsManager):
                               "eps_tol": 1e-10}
 
     def __init__(self, problem, n_evals, *, options_prefix=None,
-                 solver_parameters=None):
+                 solver_parameters=None, ncv=None, mpd=None):
 
         self.es = SLEPc.EPS().create(comm=problem.dm.comm)
         self._problem = problem
         self.n_evals = n_evals
+        self.ncv = ncv
+        self.mpd = mpd
         solver_parameters = flatten_parameters(solver_parameters or {})
         for key in self.DEFAULT_EPS_PARAMETERS:
             value = self.DEFAULT_EPS_PARAMETERS[key]
@@ -132,9 +140,8 @@ def solve(self):
             weight=self._problem.bc_shift and 1./self._problem.bc_shift
         ).M.handle
 
+        self.es.setDimensions(nev=self.n_evals, ncv=self.ncv, mpd=self.mpd)
         self.es.setOperators(self.A_mat, self.M_mat)
-        # SLEPc recommended params
-        self.es.setDimensions(nev=self.n_evals, ncv=2*self.n_evals)
         with self.inserted_options():
             self.es.solve()
         nconv = self.es.getConverged()

From c7ccd21433705902566d0a7110b225e98563a060 Mon Sep 17 00:00:00 2001
From: David Ham <David.Ham@imperial.ac.uk>
Date: Fri, 23 Jun 2023 15:56:28 +0100
Subject: [PATCH 11/22] port eigensolver demo

---
 .../qgbasinmodes.py.rst                       | 157 +++++++-----------
 1 file changed, 61 insertions(+), 96 deletions(-)

diff --git a/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst b/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst
index 1f4bf8f6c2..0569eebd9d 100644
--- a/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst
+++ b/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst
@@ -6,20 +6,16 @@ Oceanic Basin Modes: Quasi-Geostrophic approach
    This tutorial was contributed by Christine Kaufhold and `Francis
    Poulin <mailto:fpoulin@uwaterloo.ca>`__.
 
-As a continuation of the Quasi-Geostrophic (QG) model described in the
-other tutorial, we will now see how we can use Firedrake to compute
-the spatial structure and frequencies of the freely evolving modes in this system,
-what are referred to as basin modes.
-Oceanic basin modes are low frequency structures that propagate
-zonally in the oceans that alter the dynamics of Western Boundary Currents,
-such as the Gulf Stream. In this particular tutorial we will show how to
-solve the QG eigenvalue problem with no basic state and no dissipative
-forces.
-Unlike the other demo that integrated the equations forward in time, in this
-problem it is necessary to compute the eigenvalues and eigenfunctions
-for a particular differential operator. This requires using
-`PETSc <http://www.mcs.anl.gov/petsc/>`__ matrices
-and eigenvalue solvers in `SLEPc <http://slepc.upv.es>`__.
+As a continuation of the Quasi-Geostrophic (QG) model described in the other
+tutorial, we will now see how we can use Firedrake to compute the spatial
+structure and frequencies of the freely evolving modes in this system, what are
+referred to as basin modes. Oceanic basin modes are low frequency structures
+that propagate zonally in the oceans that alter the dynamics of Western
+Boundary Currents, such as the Gulf Stream. In this particular tutorial we will
+show how to solve the QG eigenvalue problem with no basic state and no
+dissipative forces. Unlike the other demo that integrated the equations forward
+in time, in this problem it is necessary to compute the eigenvalues and
+eigenfunctions for a particular differential operator.
 
 This demo requires SLEPc and slepc4py to be installed.  This is most easily
 achieved by providing the optional `--slepc` flag to either `firedrake-install`
@@ -105,133 +101,102 @@ Firedrake code
 --------------
 
 Using this form, we can now implement this eigenvalue problem in
-Firedrake. We import the Firedrake, PETSc, and SLEPc libraries. ::
-
-   from firedrake import *
-   from firedrake.petsc import PETSc
-   try:
-       from slepc4py import SLEPc
-   except ImportError:
-       import sys
-       warning("Unable to import SLEPc, eigenvalue computation not possible (try firedrake-update --slepc)")
-       sys.exit(0)
+Firedrake. We start by importing Firedrake. ::
 
+  from firedrake import *
 
 We specify the geometry to be a square geometry with :math:`50` cells
 with length :math:`1`. ::
 
-   Lx   = 1.
-   Ly   = 1.
-   n0   = 50
-   mesh = RectangleMesh(n0, n0, Lx, Ly, reorder=None)
+  Lx   = 1.
+  Ly   = 1.
+  n0   = 50
+  mesh = RectangleMesh(n0, n0, Lx, Ly, reorder=None)
 
 Next we define the function spaces within which our solution will
 reside. ::
 
-   Vcg  = FunctionSpace(mesh,'CG',3)
+  Vcg  = FunctionSpace(mesh,'CG',3)
 
 We impose zero Dirichlet boundary conditions, in a strong sense, which
 guarantee that we have no-normal flow at the boundary walls. ::
 
-   bc = DirichletBC(Vcg, 0.0, "on_boundary")
+  bc = DirichletBC(Vcg, 0.0, "on_boundary")
 
 The two non-dimensional parameters are the :math:`\beta` parameter, set
 by the sphericity of the Earth, and the Froude number, the relative
 importance of rotation to stratification. ::
 
-   beta = Constant('1.0')
-   F    = Constant('1.0')
-
-Additionally, we can create some Functions to store the eigenmodes. ::
-
-   eigenmodes_real, eigenmodes_imag = Function(Vcg), Function(Vcg)
+  beta = Constant('1.0')
+  F    = Constant('1.0')
 
 We define the Test Function :math:`\phi` and the Trial Function
 :math:`\psi` in our function space. ::
 
-   phi, psi = TestFunction(Vcg), TrialFunction(Vcg)
+  phi, psi = TestFunction(Vcg), TrialFunction(Vcg)
 
 To build the weak formulation of our equation we need to build two PETSc
 matrices in the form of a generalized eigenvalue problem,
-:math:`A\psi = \lambda M\psi`. We impose the boundary conditions on the
-mass matrix :math:`M`, since that is where we used integration by parts. ::
-
-   a =  beta*phi*psi.dx(0)*dx
-   m = -inner(grad(psi), grad(phi))*dx - F*psi*phi*dx
-   petsc_a = assemble(a).M.handle
-   petsc_m = assemble(m, bcs=bc).M.handle
-
-We can declare how many eigenpairs, eigenfunctions and eigenvalues, we
-want to find ::
+:math:`A\psi = \lambda M\psi` ::
 
-   num_eigenvalues = 1
+  eigenproblem = LinearEigenproblem(
+          A=beta*phi*psi.dx(0)*dx,
+          M=-inner(grad(psi), grad(phi))*dx - F*psi*phi*dx,
+          bcs=bc)
 
-Next we will impose parameters onto our eigenvalue solver. The first is
-specifying that we have an generalized eigenvalue problem that is
-nonhermitian. The second specifies the spectral transform shift factor
-to be non-zero. The third requires we use a Krylov-Schur method,
-which is the default so this is not strictly necessary. Then, we ask for
+Next we program our eigenvalue solver through the PETSc options system. The
+first is specifying that we have an generalized eigenvalue problem that is
+nonhermitian. The second specifies the spectral transform shift factor to be
+non-zero. The third requires we use a Krylov-Schur method. Then, we ask for
 the eigenvalues with the largest imaginary part. Finally, we specify the
-tolerance. ::
+tolerance.
 
-   opts = PETSc.Options()
-   opts.setValue("eps_gen_non_hermitian", None)
-   opts.setValue("st_pc_factor_shift_type", "NONZERO")
-   opts.setValue("eps_type", "krylovschur")
-   opts.setValue("eps_largest_imaginary", None)
-   opts.setValue("eps_tol", 1e-10)
+These options are the defaults used by Firedrake's eigensolver, but they are
+provided here to illustrate the proces of setting parameters::
 
-Finally, we build our eigenvalue solver using SLEPc. We add our PETSc
-matrices into the solver as operators and use setFromOptions() to call
-the PETSc parameters we previously declared. ::
+  opts = {"eps_gen_non_hermitian": None,
+          "st_pc_factor_shift_type": "NONZERO",
+          "eps_type": "krylovschur",
+          "eps_largest_imaginary": None,
+          "eps_tol": 1e-10}
 
-   es = SLEPc.EPS().create(comm=COMM_WORLD)
-   es.setDimensions(num_eigenvalues)
-   es.setOperators(petsc_a, petsc_m)
-   es.setFromOptions()
-   es.solve()
+Finally, we build our eigenvalue solver, specifying in this case that we just
+want to see the first eigenvalue, eigenvector pair::
 
-Additionally we can find the number of converged eigenvalues. ::
+  eigensolver = LinearEigensolver(eigenproblem, n_evals=1,
+                                  solver_parameters=opts)
 
-   nconv = es.getConverged()
+Now solve the system. This returns the number of converged eigenvalues. ::
+
+  nconv = eigensolver.solve()
 
 We now get the real and imaginary parts of the eigenvalue and
 eigenvector for the leading eigenpair (that with the largest in
-magnitude imaginary part).  First we check if we actually managed to
-converge any eigenvalues at all. ::
-
-   if nconv == 0:
-       import sys
-       warning("Did not converge any eigenvalues")
-       sys.exit(0)
-
-If we did, we go ahead and extract them from the SLEPc eigenvalue
-solver::
-
-   vr, vi = petsc_a.getVecs()
+magnitude imaginary part). ::
 
-   lam = es.getEigenpair(0, vr, vi)
+  lam = eigensolver.eigenvalue(0)
 
-and we gather the final eigenfunctions ::
+and we gather the corresponding eigenfunctions ::
 
-   eigenmodes_real.vector()[:], eigenmodes_imag.vector()[:] = vr, vi
+  eigenmode_real, eigenmode_imag = eigensolver.eigenfunction(0)
 
 We can now list and show plots for the eigenvalues and eigenfunctions
 that were found. ::
 
-   print("Leading eigenvalue is:", lam)
+  print("Leading eigenvalue is:", lam)
 
-   try:
-       import matplotlib.pyplot as plt
-       fig, axes = plt.subplots()
-       colors = tripcolor(eigenmodes_real, axes=axes)
-       fig.colorbar(colors)
+  try:
+      import matplotlib.pyplot as plt
+      fig, axes = plt.subplots()
+      colors = tripcolor(eigenmode_real, axes=axes)
+      fig.colorbar(colors)
 
-       fig, axes = plt.subplots()
-       colors = tripcolor(eigenmodes_imag, axes=axes)
-       fig.colorbar(colors)
-   except ImportError:
-       warning("Matplotlib not available, not plotting eigemodes")
+      fig, axes = plt.subplots()
+      colors = tripcolor(eigenmode_imag, axes=axes)
+      fig.colorbar(colors)
+      plt.show()
+  except ImportError:
+      warning("Matplotlib not available, not plotting eigemodes")
 
 Below is a plot of the spatial structure of the real part of one of the eigenmodes computed above.
 

From 29592282443e0f9545778322f216c980341cbd00 Mon Sep 17 00:00:00 2001
From: David Ham <David.Ham@imperial.ac.uk>
Date: Fri, 23 Jun 2023 16:06:15 +0100
Subject: [PATCH 12/22] Add Yian to the contributors

---
 docs/source/team.ini | 1 +
 1 file changed, 1 insertion(+)

diff --git a/docs/source/team.ini b/docs/source/team.ini
index 12388d8be9..f0b79413db 100644
--- a/docs/source/team.ini
+++ b/docs/source/team.ini
@@ -111,4 +111,5 @@ Ben Sepanski: https://bensepanski.github.io
 Jemma Shipton:
 Joseph G. Wallwork: https://www.imperial.ac.uk/people/j.wallwork16
 Florian Wechsung: https://florianwechsung.github.io
+Yian Zeng:
 Fangyi Zhou:

From 65f93a4544b72da3a35a1823eb68bad4716c6706 Mon Sep 17 00:00:00 2001
From: "David A. Ham" <david.ham@imperial.ac.uk>
Date: Mon, 26 Jun 2023 14:33:17 +0100
Subject: [PATCH 13/22] Apply suggestions from code review

Co-authored-by: ksagiyam <46749170+ksagiyam@users.noreply.github.com>
---
 firedrake/eigensolver.py | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/firedrake/eigensolver.py b/firedrake/eigensolver.py
index 523532158d..69b5c90611 100644
--- a/firedrake/eigensolver.py
+++ b/firedrake/eigensolver.py
@@ -22,13 +22,13 @@ class LinearEigenproblem():
     Parameters
     ----------
     A : ufl.Form
-        the bilinear form A(u, v)
+        The bilinear form A(u, v).
     M : ufl.Form
-        the mass form M(u, v),  defaults to u * v * dx.
+        The mass form M(u, v),  defaults to u * v * dx.
     bcs : DirichletBC or list of DirichletBC
-        the boundary conditions
+        The boundary conditions.
     bc_shift: float
-        the value to shift the boundary condition eigenvalues by.
+        The value to shift the boundary condition eigenvalues by.
 
     Notes
     -----
@@ -48,7 +48,7 @@ def __init__(self, A, M=None, bcs=None, bc_shift=666.0):
 
         self.A = A  # LHS
         args = A.arguments()
-        v, u = args[0], args[1]
+        v, u = args
         if M:
             self.M = M
         else:

From c2fb10caa669c874e98ad9a0e9c011b00071d135 Mon Sep 17 00:00:00 2001
From: "David A. Ham" <David.Ham@imperial.ac.uk>
Date: Tue, 27 Jun 2023 17:12:45 +0100
Subject: [PATCH 14/22] WIP eigensolver options

---
 .../qgbasinmodes.py.rst                       | 14 +++-------
 firedrake/eigensolver.py                      | 26 ++++++++++++++++---
 2 files changed, 25 insertions(+), 15 deletions(-)

diff --git a/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst b/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst
index 0569eebd9d..927c7018cb 100644
--- a/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst
+++ b/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst
@@ -146,19 +146,11 @@ matrices in the form of a generalized eigenvalue problem,
 
 Next we program our eigenvalue solver through the PETSc options system. The
 first is specifying that we have an generalized eigenvalue problem that is
-nonhermitian. The second specifies the spectral transform shift factor to be
-non-zero. The third requires we use a Krylov-Schur method. Then, we ask for
-the eigenvalues with the largest imaginary part. Finally, we specify the
-tolerance.
-
-These options are the defaults used by Firedrake's eigensolver, but they are
-provided here to illustrate the proces of setting parameters::
+nonhermitian. Then, we ask for the eigenvalues with the largest imaginary
+part::
 
   opts = {"eps_gen_non_hermitian": None,
-          "st_pc_factor_shift_type": "NONZERO",
-          "eps_type": "krylovschur",
-          "eps_largest_imaginary": None,
-          "eps_tol": 1e-10}
+          "eps_largest_imaginary": None}
 
 Finally, we build our eigenvalue solver, specifying in this case that we just
 want to see the first eigenvalue, eigenvector pair::
diff --git a/firedrake/eigensolver.py b/firedrake/eigensolver.py
index 69b5c90611..7925d6200d 100644
--- a/firedrake/eigensolver.py
+++ b/firedrake/eigensolver.py
@@ -13,7 +13,7 @@
 
 
 class LinearEigenproblem():
-    """Linear eigenvalue problem
+    """Generalised linear eigenvalue problem
 
     The problem has the form::
 
@@ -88,12 +88,30 @@ class LinearEigensolver(OptionsManager):
     mpd : int
         Maximum dimension allowed for the projected problem. See
         `SLEPc.EPS.setDimensions`.
+
+    Notes
+    -----
+
+    Users will typically wish to set solver parameters specifying the symmetry
+    of the eigenproblem and which eigenvalues to search for first.
+
+    The former is set using the options available for `EPSSetProblemType
+    <https://slepc.upv.es/documentation/current/docs/manualpages/EPS/EPSSetProblemType.html>`__. 
+
+    For example if the bilinear form is symmetric (Hermitian in complex mode),
+    one would add this entry to `solver_options`::
+
+        "eps_gen_hermitian": None
+
+    As always when specifying PETSc options, `None` indicates that the option
+    in question is a flag and hence doesn't take an argument.
+
+    The eigenvalues to search for first are specified using the options for
+    `EPSSetWhichEigenPairs`. FIXME
     """
 
-    DEFAULT_EPS_PARAMETERS = {"eps_gen_non_hermitian": None,
-                              "st_pc_factor_shift_type": "NONZERO",
+    DEFAULT_EPS_PARAMETERS = {"st_pc_factor_shift_type": "NONZERO",
                               "eps_type": "krylovschur",
-                              "eps_largest_imaginary": None,
                               "eps_tol": 1e-10}
 
     def __init__(self, problem, n_evals, *, options_prefix=None,

From a9ba32ca1bacf7c6d49af845c69b4ed99a4e3ea8 Mon Sep 17 00:00:00 2001
From: "David A. Ham" <david.ham@imperial.ac.uk>
Date: Tue, 4 Jul 2023 09:21:16 +0100
Subject: [PATCH 15/22] Update tests/regression/test_eigensolver.py

Co-authored-by: ksagiyam <46749170+ksagiyam@users.noreply.github.com>
---
 tests/regression/test_eigensolver.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/tests/regression/test_eigensolver.py b/tests/regression/test_eigensolver.py
index 2ebc8cbea6..ac84d463fb 100644
--- a/tests/regression/test_eigensolver.py
+++ b/tests/regression/test_eigensolver.py
@@ -18,10 +18,10 @@ def evals(n, degree=1, mesh=None):
 
     # Create eigenproblem with boundary conditions
     bc = DirichletBC(V, 0.0, "on_boundary")
-    eigenprob = LinearEigenproblem(a, bcs=bc)
+    eigenprob = LinearEigenproblem(a, bcs=bc, bc_shift=-666.)
 
     # Create corresponding eigensolver, looking for n eigenvalues
-    eigensolver = LinearEigensolver(eigenprob, n)
+    eigensolver = LinearEigensolver(eigenprob, n, solver_parameters={"eps_largest_real": None})
     ncov = eigensolver.solve()
 
     # boffi solns

From 963048b0ad12cd130d787634589d72cd00528f4d Mon Sep 17 00:00:00 2001
From: David Ham <David.Ham@imperial.ac.uk>
Date: Tue, 4 Jul 2023 09:57:34 +0100
Subject: [PATCH 16/22] finish fixing docs.

---
 firedrake/eigensolver.py             | 15 ++++++++++-----
 tests/regression/test_eigensolver.py |  4 +++-
 2 files changed, 13 insertions(+), 6 deletions(-)

diff --git a/firedrake/eigensolver.py b/firedrake/eigensolver.py
index 7925d6200d..ced8a5445b 100644
--- a/firedrake/eigensolver.py
+++ b/firedrake/eigensolver.py
@@ -13,18 +13,18 @@
 
 
 class LinearEigenproblem():
-    """Generalised linear eigenvalue problem
+    """Generalised linear eigenvalue problem.
 
-    The problem has the form::
+    The problem has the form, find `u`, `λ` such that::
 
-        A(u, v) = λ * M(u, v)
+        A(u, v) = λM(u, v)    ∀ v ∈ V
 
     Parameters
     ----------
     A : ufl.Form
         The bilinear form A(u, v).
     M : ufl.Form
-        The mass form M(u, v),  defaults to u * v * dx.
+        The mass form M(u, v),  defaults to inner(u, v) * dx.
     bcs : DirichletBC or list of DirichletBC
         The boundary conditions.
     bc_shift: float
@@ -107,7 +107,12 @@ class LinearEigensolver(OptionsManager):
     in question is a flag and hence doesn't take an argument.
 
     The eigenvalues to search for first are specified using the options for
-    `EPSSetWhichEigenPairs`. FIXME
+    `EPSSetWhichEigenPairs <https://slepc.upv.es/documentation/current/docs/manualpages/EPS/EPSSetWhichEigenpairs.html>`__.
+
+    For example, to look for the eigenvalues with largest real part, one
+    would add this entry to `solver_options`::
+
+        "eps_largest_magnitude": None
     """
 
     DEFAULT_EPS_PARAMETERS = {"st_pc_factor_shift_type": "NONZERO",
diff --git a/tests/regression/test_eigensolver.py b/tests/regression/test_eigensolver.py
index ac84d463fb..04c724e1e4 100644
--- a/tests/regression/test_eigensolver.py
+++ b/tests/regression/test_eigensolver.py
@@ -21,7 +21,9 @@ def evals(n, degree=1, mesh=None):
     eigenprob = LinearEigenproblem(a, bcs=bc, bc_shift=-666.)
 
     # Create corresponding eigensolver, looking for n eigenvalues
-    eigensolver = LinearEigensolver(eigenprob, n, solver_parameters={"eps_largest_real": None})
+    eigensolver = LinearEigensolver(
+        eigenprob, n, solver_parameters={"eps_largest_real": None}
+    )
     ncov = eigensolver.solve()
 
     # boffi solns

From 9ab8a439c33fe5575a30cad68807284bd484a110 Mon Sep 17 00:00:00 2001
From: David Ham <David.Ham@imperial.ac.uk>
Date: Tue, 4 Jul 2023 11:55:49 +0100
Subject: [PATCH 17/22] lint

---
 firedrake/eigensolver.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/firedrake/eigensolver.py b/firedrake/eigensolver.py
index ced8a5445b..636197e944 100644
--- a/firedrake/eigensolver.py
+++ b/firedrake/eigensolver.py
@@ -96,7 +96,7 @@ class LinearEigensolver(OptionsManager):
     of the eigenproblem and which eigenvalues to search for first.
 
     The former is set using the options available for `EPSSetProblemType
-    <https://slepc.upv.es/documentation/current/docs/manualpages/EPS/EPSSetProblemType.html>`__. 
+    <https://slepc.upv.es/documentation/current/docs/manualpages/EPS/EPSSetProblemType.html>`__.
 
     For example if the bilinear form is symmetric (Hermitian in complex mode),
     one would add this entry to `solver_options`::

From 2356bf3a77528332a78e262dfbbeee5e88ed806f Mon Sep 17 00:00:00 2001
From: "David A. Ham" <david.ham@imperial.ac.uk>
Date: Tue, 4 Jul 2023 12:10:34 +0100
Subject: [PATCH 18/22] Update firedrake/eigensolver.py

Co-authored-by: ksagiyam <46749170+ksagiyam@users.noreply.github.com>
---
 firedrake/eigensolver.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/firedrake/eigensolver.py b/firedrake/eigensolver.py
index 636197e944..f6726d9d56 100644
--- a/firedrake/eigensolver.py
+++ b/firedrake/eigensolver.py
@@ -112,7 +112,7 @@ class LinearEigensolver(OptionsManager):
     For example, to look for the eigenvalues with largest real part, one
     would add this entry to `solver_options`::
 
-        "eps_largest_magnitude": None
+        "eps_largest_real": None
     """
 
     DEFAULT_EPS_PARAMETERS = {"st_pc_factor_shift_type": "NONZERO",

From 95e6b53b20456ec0f3312245b603c63b4fb92eb5 Mon Sep 17 00:00:00 2001
From: "David A. Ham" <David.Ham@imperial.ac.uk>
Date: Wed, 5 Jul 2023 15:31:35 +0100
Subject: [PATCH 19/22] change default BC handling

---
 firedrake/eigensolver.py             | 9 ++++++---
 tests/regression/test_eigensolver.py | 2 +-
 2 files changed, 7 insertions(+), 4 deletions(-)

diff --git a/firedrake/eigensolver.py b/firedrake/eigensolver.py
index f6726d9d56..72d4acfe09 100644
--- a/firedrake/eigensolver.py
+++ b/firedrake/eigensolver.py
@@ -39,7 +39,7 @@ class LinearEigenproblem():
     amount. It is the user's responsibility to ensure that the shift is not
     close to an actual eigenvalue of the system.
     """
-    def __init__(self, A, M=None, bcs=None, bc_shift=666.0):
+    def __init__(self, A, M=None, bcs=None, bc_shift=0.0):
         if not SLEPc:
             raise ImportError(
                 "Unable to import SLEPc, eigenvalue computation not possible "
@@ -64,7 +64,7 @@ def dirichlet_bcs(self):
             yield from bc.dirichlet_bcs()
 
     @utils.cached_property
-    def dm(self):  # cargo cult
+    def dm(self):
         r"""Return the dm associated with the output space."""
         return self.output_space.dm
 
@@ -117,7 +117,8 @@ class LinearEigensolver(OptionsManager):
 
     DEFAULT_EPS_PARAMETERS = {"st_pc_factor_shift_type": "NONZERO",
                               "eps_type": "krylovschur",
-                              "eps_tol": 1e-10}
+                              "eps_tol": 1e-10,
+                              "eps_target": 0.0}
 
     def __init__(self, problem, n_evals, *, options_prefix=None,
                  solver_parameters=None, ncv=None, mpd=None):
@@ -131,6 +132,8 @@ def __init__(self, problem, n_evals, *, options_prefix=None,
         for key in self.DEFAULT_EPS_PARAMETERS:
             value = self.DEFAULT_EPS_PARAMETERS[key]
             solver_parameters.setdefault(key, value)
+        if self._problem.bcs:
+            solver_parameters.setdefault("st_type", "sinvert")
         super().__init__(solver_parameters, options_prefix)
         self.set_from_options(self.es)
 
diff --git a/tests/regression/test_eigensolver.py b/tests/regression/test_eigensolver.py
index 04c724e1e4..657e623408 100644
--- a/tests/regression/test_eigensolver.py
+++ b/tests/regression/test_eigensolver.py
@@ -18,7 +18,7 @@ def evals(n, degree=1, mesh=None):
 
     # Create eigenproblem with boundary conditions
     bc = DirichletBC(V, 0.0, "on_boundary")
-    eigenprob = LinearEigenproblem(a, bcs=bc, bc_shift=-666.)
+    eigenprob = LinearEigenproblem(a, bcs=bc, bc_shift=-6666.)
 
     # Create corresponding eigensolver, looking for n eigenvalues
     eigensolver = LinearEigensolver(

From aa488a9e735ddc0401505775c0bf00f8bdbc677f Mon Sep 17 00:00:00 2001
From: "David A. Ham" <David.Ham@imperial.ac.uk>
Date: Wed, 5 Jul 2023 15:47:11 +0100
Subject: [PATCH 20/22] unbreak demo

---
 demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst | 6 ++++--
 firedrake/eigensolver.py                            | 2 +-
 2 files changed, 5 insertions(+), 3 deletions(-)

diff --git a/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst b/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst
index 927c7018cb..8d43d49aaf 100644
--- a/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst
+++ b/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst
@@ -147,10 +147,12 @@ matrices in the form of a generalized eigenvalue problem,
 Next we program our eigenvalue solver through the PETSc options system. The
 first is specifying that we have an generalized eigenvalue problem that is
 nonhermitian. Then, we ask for the eigenvalues with the largest imaginary
-part::
+part. Finally we set the spectral transform to shift with no target::
 
   opts = {"eps_gen_non_hermitian": None,
-          "eps_largest_imaginary": None}
+          "eps_largest_imaginary": None,
+          "st_type": "shift",
+          "eps_target": None}
 
 Finally, we build our eigenvalue solver, specifying in this case that we just
 want to see the first eigenvalue, eigenvector pair::
diff --git a/firedrake/eigensolver.py b/firedrake/eigensolver.py
index 72d4acfe09..ea2d3a5a91 100644
--- a/firedrake/eigensolver.py
+++ b/firedrake/eigensolver.py
@@ -117,7 +117,7 @@ class LinearEigensolver(OptionsManager):
 
     DEFAULT_EPS_PARAMETERS = {"st_pc_factor_shift_type": "NONZERO",
                               "eps_type": "krylovschur",
-                              "eps_tol": 1e-10,
+                              "eps_tol": 1e-10, 
                               "eps_target": 0.0}
 
     def __init__(self, problem, n_evals, *, options_prefix=None,

From 73a7dbcd65c0dd023d856bc7f1c9cc1d20d4c0ad Mon Sep 17 00:00:00 2001
From: "David A. Ham" <David.Ham@imperial.ac.uk>
Date: Wed, 5 Jul 2023 15:58:27 +0100
Subject: [PATCH 21/22] move the factor shift back to the demo

---
 demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst | 3 ++-
 firedrake/eigensolver.py                            | 5 ++---
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst b/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst
index 8d43d49aaf..fb6c299d1e 100644
--- a/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst
+++ b/demos/eigenvalues_QG_basinmodes/qgbasinmodes.py.rst
@@ -152,7 +152,8 @@ part. Finally we set the spectral transform to shift with no target::
   opts = {"eps_gen_non_hermitian": None,
           "eps_largest_imaginary": None,
           "st_type": "shift",
-          "eps_target": None}
+          "eps_target": None,
+          "st_pc_factor_shift_type": "NONZERO"}
 
 Finally, we build our eigenvalue solver, specifying in this case that we just
 want to see the first eigenvalue, eigenvector pair::
diff --git a/firedrake/eigensolver.py b/firedrake/eigensolver.py
index ea2d3a5a91..c883bf96be 100644
--- a/firedrake/eigensolver.py
+++ b/firedrake/eigensolver.py
@@ -115,9 +115,8 @@ class LinearEigensolver(OptionsManager):
         "eps_largest_real": None
     """
 
-    DEFAULT_EPS_PARAMETERS = {"st_pc_factor_shift_type": "NONZERO",
-                              "eps_type": "krylovschur",
-                              "eps_tol": 1e-10, 
+    DEFAULT_EPS_PARAMETERS = {"eps_type": "krylovschur",
+                              "eps_tol": 1e-10,
                               "eps_target": 0.0}
 
     def __init__(self, problem, n_evals, *, options_prefix=None,

From d80d81984f3053b3248d53d5d0f0fa68e233d027 Mon Sep 17 00:00:00 2001
From: "David A. Ham" <David.Ham@imperial.ac.uk>
Date: Mon, 10 Jul 2023 13:55:55 +0100
Subject: [PATCH 22/22] 2d test

---
 tests/regression/test_eigensolver.py | 28 ++++++++++++++++++++++++++++
 1 file changed, 28 insertions(+)

diff --git a/tests/regression/test_eigensolver.py b/tests/regression/test_eigensolver.py
index 657e623408..f58339238f 100644
--- a/tests/regression/test_eigensolver.py
+++ b/tests/regression/test_eigensolver.py
@@ -52,3 +52,31 @@ def evals(n, degree=1, mesh=None):
 def test_evals_1d(n, degree, tolerance):
     true_values, estimates = evals(n, degree=degree)
     assert np.allclose(true_values, estimates, rtol=tolerance)
+
+
+def poisson_eigenvalue_2d(i):
+    mesh = RectangleMesh(10*2**i, 10*2**i, pi, pi)
+    V = FunctionSpace(mesh, "CG", 1)
+    u = TrialFunction(V)
+    v = TestFunction(V)
+
+    bc = DirichletBC(V, 0.0, "on_boundary")
+
+    ep = LinearEigenproblem(inner(grad(u), grad(v)) * dx,
+                            bcs=bc, bc_shift=666.0)
+
+    es = LinearEigensolver(ep, 1, solver_parameters={"eps_gen_hermitian": None,
+                                                     "eps_largest_real": None})
+
+    es.solve()
+    return es.eigenvalue(0)-2.0
+
+
+def test_evals_2d():
+    """2D Eigenvalue convergence test. As with Boffi, we observe that the
+    convergence rate convergest to 2 from above."""
+    errors = np.array([poisson_eigenvalue_2d(i) for i in range(5)])
+
+    convergence = np.log(errors[:-1]/errors[1:])/np.log(2.0)
+
+    assert all(convergence > 2.0)