[SciPy-dev] help: wrapping generalized symmetric evp functions

[Robert Cimrman]

Pearu Peterson wrote:
> 
> Robert Cimrman wrote:
>> I am now trying to add the following LAPACK functions into scipy:
>> ssygv, dsygv, chegv, zhegv, to be able to solve 'Ax = lambda Bx' with 
>> symmetric or Hermitian matrices efficiently (cf. 
>> http://www.netlib.org/lapack/lug/node34.html).
> 
> New wrappers to lapack functions should be defined in scipy.lib.lapack.
> Eventually we should get rid of lapack wrappers in scipy.linalg.
> Grep for flapack_*.pyf.src files to see in which file the signature
> for a given lapack function should be inserted.
> 
> Btw, scipy.lib.lapack already has wrappers for sygv and hegv functions.

*** to Pearu:
Great, thanks!

Does it mean that scipy.linalg should be updated to use the wrappers in 
scipy.lib.lapack?

*** to all:
IMHO scipy.linalg.decomp deserves a major refitting in the spirit of 
what is now going on with scipy.sparse - the arguments of the functions 
with same functionality have different names (both comparing to sparse 
equivalents and within decomp.py), there are too many eig* functions, etc.

There could be even proper solver classes (LinearSolver, 
EigenvalueSolver, etc.) with a common interface in each solver family. 
Below, as an example, is the code we use in sfepy. As you can see, the 
solver classes have only the __init__ method (can be used e.g. to 
pre-factorize a matrix in the case of a linear solver) and the __call__ 
method (application of the solver). Would it be interesting to have 
something in that spirit in scipy?

Then, in the meantime, in lobpcg, I will use directly something like:
import scipy.lib.lapack as ll
ll.get_lapack_funcs( ['hegv'] )

best regards,
r.

----
class Solver( Struct ):
     def __init__( self, conf, **kwargs ):
         Struct.__init__( self, conf = conf, **kwargs )
     def __call__( self, **kwargs ):
         print 'called an abstract Solver instance!'
         raise ValueError

class LinearSolver( Solver ):
     def __init__( self, conf, mtx = None, status = None, **kwargs ):
         Solver.__init__( self, conf = conf, mtx = mtx, status = status,
                          **kwargs )
     def __call__( self, rhs, conf = None, mtx = None, status = None ):
         print 'called an abstract LinearSolver instance!'
         raise ValueError

class NonlinearSolver( Solver ):
     def __init__( self, conf, evaluator = None, linSolver = None,
                   status = None, **kwargs ):
         Solver.__init__( self, conf = conf, evaluator = evaluator,
                          linSolver = linSolver, status = status,
                          **kwargs )
     def __call__( self, state0, conf = None, evaluator = None,
                   linSolver = None, status = None ):
         print 'called an abstract NonlinearSolver instance!'
         raise ValueError

class EigenvalueSolver( Solver ):
     def __init__( self, conf, mtxA = None, mtxB = None, nEigs = None,
                   eigenvectors = None, status = None ):
         Solver.__init__( self, conf = conf, mtxA = mtxA, mtxB = mtxB,
                          nEigs = nEigs, eigenvectors = eigenvectors,
                          status = status )
     def __call__( self, mtxA, mtxB = None, nEigs = None,
                   eigenvectors = None, status = None, conf = None ):
         print 'called an abstract EigenvalueSolver instance!'
         raise ValueError