[Scipy-tickets] [SciPy] #231: A wrapper for ARPACK

SciPy scipy-tickets at scipy.net
Thu Oct 12 07:03:43 CDT 2006

#231: A wrapper for ARPACK
 Reporter:  nils          |        Owner:  somebody
     Type:  enhancement   |       Status:  new     
 Priority:  normal        |    Milestone:          
Component:  scipy.sparse  |      Version:          
 Severity:  normal        |   Resolution:          
 Keywords:                |  
Comment (by NeilenMarais):

 aric: Great news about tracking down the cause of the wonky results!

 I'm busy with RealWork (TM) ATM but will take a look again early next

 == Intent(inout) vs intent(inplace) ==

 If the numpy array passed to the fortran function is of exactly the right
 type (i.e. same dtype and fortran ordering for multi-dimensional arrays),
 both intent(inout) and intent(inplace) just passes a pointer to the Fortan
 routine. The Fortran routine then operates directly on the numpy array,
 and no data is copied. This is obviously the fastest/most memory efficient
 way to go.

 If the numpy array is of the wrong type, intent(inout) will raise an error
 (not sure if it catches all possible cases), whereas intent(inplace) will
 copy/convert the input array and again convert and copy the fortran output
 over the original input array.

 Which one to choose? Unless performance is __really__ critical and you
 want to be 100% sure that no copying is happening behind your back,
 intent(inplace) ought to prevent problems.

 == Easy Python Interface ==

 There should probably be an easy and an advanced interface. The easy
 interface could look like this:

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

 === 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.


 (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 Inteface ==

 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
     eigenvalues, nev eigenvalues close to sigma are calculated. The user
     to provide routines that calculate [M]*x and solve [A]-sigma*[M]*x = b
 for x.

     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

     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 inteface is missing a spectrum parameter but that should
 be easy to add.

 == The end ==

 Phew, quite a marathon piece. Hope it can help someone :) Anyway, if we
 can get the f2py wrappers to perform reliably, I'd be more than willing to
 spend some time on the  higher level python interfaces.

Ticket URL: <http://projects.scipy.org/scipy/scipy/ticket/231#comment:7>
SciPy <http://www.scipy.org/>
SciPy is open-source software for mathematics, science, and engineering.

More information about the Scipy-tickets mailing list