Skip to content

Commit

Permalink
Reuse the LU transform when solving for left eigenvectors using the s…
Browse files Browse the repository at this point in the history
…parse solver (#277)

* Modify sparse evp solver to reuse the LU transform from the right evp solve.

* Update code to remove blank lines, correct conj(sigma) in comment, and tidy up if condition.

* Change solve_adjoint -> solve_H and evals_adjoint->left_evals for clarity.

* Consolidate solve_H logic

---------

Co-authored-by: Keaton J. Burns <[email protected]>
  • Loading branch information
csskene and kburns authored Jan 8, 2024
1 parent 19a7ada commit 68dcb74
Show file tree
Hide file tree
Showing 3 changed files with 35 additions and 9 deletions.
15 changes: 8 additions & 7 deletions dedalus/core/solvers.py
Original file line number Diff line number Diff line change
Expand Up @@ -250,15 +250,13 @@ def solve_sparse(self, subproblem, N, target, rebuild_matrices=False, left=False
# Solve as sparse general eigenvalue problem
A = sp.L_min
B = - sp.M_min
# Solve for the right eigenvectors
self.eigenvalues, pre_right_evecs = scipy_sparse_eigs(A=A, B=B, N=N, target=target, matsolver=self.matsolver, **kw)
self.right_eigenvectors = self.eigenvectors = sp.pre_right @ pre_right_evecs
# Solve for the right (and optionally left) eigenvectors
eig_output = scipy_sparse_eigs(A=A, B=B, left=left, N=N, target=target, matsolver=self.matsolver, **kw)

if left:
# Solve for the left eigenvectors
# Note: this definition of "left eigenvectors" is consistent with the documentation for scipy.linalg.eig
self.left_eigenvalues, pre_left_evecs = scipy_sparse_eigs(A=A.getH(), B=B.getH(),
N=N, target=np.conjugate(target),
matsolver=self.matsolver, **kw)
self.eigenvalues, pre_right_evecs, self.left_eigenvalues, pre_left_evecs = eig_output
self.right_eigenvectors = self.eigenvectors = sp.pre_right @ pre_right_evecs
self.left_eigenvectors = sp.pre_left.H @ pre_left_evecs
self.modified_left_eigenvectors = (sp.M_min @ sp.pre_right_pinv).H @ pre_left_evecs
# Check that eigenvalues match
Expand All @@ -278,6 +276,9 @@ def solve_sparse(self, subproblem, N, target, rebuild_matrices=False, left=False
norms = np.diag(pre_left_evecs.T.conj() @ sp.M_min @ pre_right_evecs)
self.left_eigenvectors /= np.conj(norms)
self.modified_left_eigenvectors /= np.conj(norms)
else:
self.eigenvalues, pre_right_evecs = eig_output
self.right_eigenvectors = self.eigenvectors = sp.pre_right @ pre_right_evecs

def set_state(self, index, subsystem):
"""
Expand Down
11 changes: 11 additions & 0 deletions dedalus/libraries/matsolvers.py
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,9 @@ def __init__(self, matrix, solver=None):
def solve(self, vector):
pass

def solve_H(self,vector):
raise NotImplementedError("%s has not implemented 'solve_H' method" % type(self))


@add_solver
class DummySolver(SolverBase):
Expand Down Expand Up @@ -145,6 +148,14 @@ def __init__(self, matrix, solver=None):
def solve(self, vector):
return self.LU.solve(vector, trans=self.trans)

def solve_H(self,vector):
if self.trans == "N":
return self.LU.solve(vector, trans="H")
elif self.trans == "H":
return self.LU.solve(vector)
elif self.trans == "T":
return np.conj(self.LU.solve(np.conj(vector)))


@add_solver
class SuperluNaturalFactorized(_SuperluFactorizedBase):
Expand Down
18 changes: 16 additions & 2 deletions dedalus/tools/array.py
Original file line number Diff line number Diff line change
Expand Up @@ -395,7 +395,7 @@ def drop_empty_rows(mat):
return mat[nonempty_rows]


def scipy_sparse_eigs(A, B, N, target, matsolver, **kw):
def scipy_sparse_eigs(A, B, left, N, target, matsolver, **kw):
"""
Perform targeted eigenmode search using the scipy/ARPACK sparse solver
for the reformulated generalized eigenvalue problem
Expand All @@ -410,6 +410,8 @@ def scipy_sparse_eigs(A, B, N, target, matsolver, **kw):
Sparse matrices for generalized eigenvalue problem
N : int
Number of eigenmodes to return
left: boolean
Whether to solve for the left eigenvectors or not
target : complex
Target σ for eigenvalue search
matsolver : matrix solver class
Expand All @@ -427,7 +429,19 @@ def matvec(x):
evals, evecs = spla.eigs(D, k=N, which='LM', sigma=None, **kw)
# Rectify eigenvalues
evals = 1 / evals + target
return evals, evecs
if left:
# Build sparse linear operator representing (A^H - conj(σ)B^H)^I B^H = C^-H B^H = left_D
# Note: left_D is not equal to D^H
def matvec_left(x):
return solver.solve_H(B.H.dot(x))
left_D = spla.LinearOperator(dtype=A.dtype, shape=A.shape, matvec=matvec_left)
# Solve using scipy sparse algorithm
left_evals, left_evecs = spla.eigs(left_D, k=N, which='LM', sigma=None, **kw)
# Rectify left eigenvalues
left_evals = 1 / left_evals + np.conj(target)
return evals, evecs, left_evals, left_evecs
else:
return evals, evecs


def interleave_matrices(matrices):
Expand Down

0 comments on commit 68dcb74

Please sign in to comment.