[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
week.
== 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',
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 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
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 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