The need for an ARPACK wrapper was raised in #231

User Interface

There should probably be an easy and an advanced interface.

Easy Python Interface

The easy interface could look like this:

speigs(A, no_eigs, B=None, spectrum='SM', properties='general', sigma=None)

Simple use

(eig_vals, eig_vecs) = speigs(A, 10)

which will calculate the 10 eigenvalues of A with smallest magnitude and the related eigenvectors. A should be a matrix like object or one that supports matvec(). The python routine should check if A is real or complex, and choose the correct ARPACK routine.

Generalised problem:

(eig_vals, eig_vecs) = speigs(A, 10, B=B)

Same as above, but solves the generalised A*eig_vec = eig_val*B*eig_vec eigenproblem. An inversion of B will be required (I think), so linalg.splu() should be used internally. This interface is perhaps a little messy since the number of eigs is in the "wrong" place. Perhaps it would be better to have a seperate function for the generalised problem.

Spectrum Selection

The user may want a different part of the spectrum. For general matrices, ARPACK gives the choice of (from the ARPACK manual)

          'LM' -> want the NEV eigenvalues of largest magnitude.
          'SM' -> want the NEV eigenvalues of smallest magnitude.
          'LR' -> want the NEV eigenvalues of largest real part.
          'SR' -> want the NEV eigenvalues of smallest real part.
          'LI' -> want the NEV eigenvalues of largest imaginary part.
          'SI' -> want the NEV eigenvalues of smallest imaginary part.

E.g.

(eig_vals, eig_vecs) = speigs(A, 10, spectrum='LR')

to get eigenvalues with largest real part.

Matrix type

Set properties='symmetrical' or somesuch to indicate that the symmetric routines are to be used. This will also prevent imaginary eigenvalues caused by roundoff for a real, symmetric system. We could also try to increase efficiency by using a symmetric solver if matrix inversion is required. The matrix type will however influence the spectrum selection codes used by ARPACK eg. 'LR' (largest real part) is replaced by 'LA' (largest algebraic). Perhaps this is inappropriate for the "simple" wrapper, but symmetric matrices are quite common and the computational savings possible quite significant.

Spectrum Shifting

Setting sigma shifts the spectrum of the orignal system by sigma. This transforms the original eigenvalues to be 1/(original_eig-sigma). IOW, eigenvalues closest to sigma will have the largest magnitude. The spectrum='XX' flag now operates on the shifted spectrum. The ARPACK postprocessing routines will transform the shifted eigenvalues/vectors back to those of the original system.

E.g. spectrum='LM' will get the eigenvalues closest to sigma, spectrum='LR' will get the eigenvalues closest to sigma, but with real part larger than sigma. This mode seems to be most usefull for vector finite element matrices (which is what I'm solving), since you typically have a generalised system with a bunch of uninteresting eigenvalues that are numerically zero (i.e. abs(eig_value) < 10e-14) and then postive real eigenvalues that you actually want.

To solve a shifted-spectrum problem, ARPACK needs us to solve (A-sigma*B)*x = B*v for x, as well as perform B*v seperately. Practically this means we need to solve (A - sigma*B) using linalg.splu(). Here I'm assuming a generalised system. To solve the normal eigenproblem, set B = I (identity matrix).

Advanced Python Interface

Sometimes you may want to use a user specified matrix-vector product or solver. E.g. the elements for the matrix vector product is calculated on the fly to save memory, or iterative solution (potentially using some specialised preconditioner) for shift-invert mode. The advanced interface could also provide more control over the type of ARPACK solver to use. My personal ARPACK python driver function looks like this:

def geneigvals(matvec, sigma_solve, n, sigma, nev, ncv=None):
    """
    Calculate eigenvalues close to sigma for generalised eigen system

    Given a system [A]x = k_i*[M]x where [A] and [M] are matrices and k_i are
    eigenvalues, nev eigenvalues close to sigma are calculated. The user needs
    to provide routines that calculate [M]*x and solve [A]-sigma*[M]*x = b for x.

    Arguments
    =========
    matvec -- Function that provides matrix-vector product, i.e. matvec(x) -> [M]*x
    sigma_solve -- sigma_solve(b) -> x, where [A]-sigma*[M]*x = b
    n -- Matrix dimension of the problem
    sigma -- Eigenvalue spectral shift real value
    nev -- Number of eigenvalues to calculate
    ncv -- Number of Arnoldi basisvectors to use. If None, default to 2*nev 

    Return Values
    =============
    Real array of nev eigenvalues, or complex array if values are complex
    """

IOW, ARPACK only interacts with the matrices through the functions matvec and sigma_solve. This leaves the user free to implement these in any desired way. My interface is missing a spectrum parameter but that should be easy to add.