The need for an ARPACK wrapper was raised in #231
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)
(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.
(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.
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.
(eig_vals, eig_vecs) = speigs(A, 10, spectrum='LR')
to get eigenvalues with largest real part.
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.
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.