[SciPy-user] eigenvalues of sparse matrix

[Travis Oliphant]

Pearu Peterson wrote:

>
>
> On Thu, 28 Oct 2004, David Grant wrote:
>
>> Travis Oliphant wrote:
>>
>>>
>>>>
>>>> Does anyone know how much work it would be to get ARPACK support in 
>>>> scipy? Even just minimal support.  This is totally new ground for 
>>>> me, so I'd like to know if I'm biting off more than I can chew.
>>>
>>>
>>>
>>>
>>> Not a lot for simple support if you use f2py.
>>>
>>> My suggestion would be to write an f2py interface for the desired 
>>> routines and we can build up from there.
>>>
>>>
>> I assume your suggestion does not involve using this at all: 
>> http://jrfonseca.dyndns.org/work/phd/ ?
>> I believe Josa Fonseca was going to use f2py, but he didn't.
>>
>> Here is from his web page:
>>
>> "These bindings weren't generated with any automatic binding 
>> generation tool. Even though I initially tried both PyFortran and 
>> f2py, both showed to be inappropriate to handle the specificity of 
>> the ARPACK API. ARPACK uses a /reverse communication interface/ where 
>> basically the API sucessively returns to caller which must take some 
>> update steps, and re-call the API with most arguments untouched. The 
>> intelligent (and silent) argument conversions made by the above tools 
>> made very difficult to implement and debug the most simple example. 
>> Also, for large-scale problems we wouldn't want any kind of array 
>> conversion/transposing happening behind the scenes as that would 
>> completely kill performance."
>
>
> I cannot agree with these conclusions, at least when using f2py. The 
> array conversion can be easily avoided behind the scenes (because the 
> main loop is in Fortran) and transposing is not even an issue as 
> user-specific parts of ARPACK API contains only 1-dimensional arrays.
>
I am in agreement with Pearu.  Look at the iterative package in 
scipy.linalg.  The underlying Fortran code uses reverse communication 
and this was managed with f2py-generated wrappers.

>> "Therefore the bindings are faithfull to the original FORTRAN API and 
>> are not very friendly as is. Nevertheless a Python wrapper to these 
>> calls can easily be made, where all kind of type conversions and 
>> dummy-safe actions can be made."
>
>
>> So who is right?
>
>
> I looked into ARPACK API and I think f2py could be effectively used to 
> wrap
> ARPACK. For example, when considering SIMPLE/dssimp.f case then
> I would write a Fortran subroutine that contains dsaupd main-loop and 
> call to dseupd (the subroutine av would be provided as a callback 
> argument) and then wrap this subroutine with f2py.

My suggestion is to do something very similar as is done in 
scipy/linalg/iterative.py   wherein the object passed to the 
Python-level code simply requires a matvec method (like all sparse 
arrays have), or is a normal 2-d array.  Then, call python-wrapped 
dsaupd in a while loop which does the required thing after returning 
from every dsaupd call.  

This is the most flexible interface and deserves a place in scipy.

Pearu's suggestion is also useful and may be a bit faster (though I'm 
not sure it would be noticeably faster). 

-Travis