Function scipy_sparse_eigs

dedalus/tools/array.py:398–444

Perform targeted eigenmode search using the scipy/ARPACK sparse solver for the reformulated generalized eigenvalue problem A.x = λ B.x ==> (A - σB)^I B.x = (1/(λ-σ)) x for eigenvalues λ near the target σ. Parameters ---------- A, B : scipy sparse matrices

(A, B, left, 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

        A.x = λ B.x  ==>  (A - σB)^I B.x = (1/(λ-σ)) x

    for eigenvalues λ near the target σ.

    Parameters
    ----------
    A, B : scipy sparse matrices
        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
        Class implementing solve method for solving sparse systems.

    Other keyword options passed to scipy.sparse.linalg.eigs.
    """
    # Build sparse linear operator representing (A - σB)^I B = C^I B = D
    C = A - target * B
    solver = matsolver(C)
    def matvec(x):
        return solver.solve(B.dot(x))
    D = spla.LinearOperator(dtype=A.dtype, shape=A.shape, matvec=matvec)
    # Solve using scipy sparse algorithm
    evals, evecs = spla.eigs(D, k=N, which='LM', sigma=None, **kw)
    # Rectify eigenvalues
    evals = 1 / evals + target
    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.conj().T.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):