[SciPy-user] RREF - upper triangular

Nils Wagner nwagner at mecha.uni-stuttgart.de
Thu Feb 23 02:17:57 CST 2006


arnd.baecker at web.de wrote:
>On Wed, 22 Feb 2006, Alan G Isaac wrote:
>
>  
>>On Wed, 22 Feb 2006, Nils Wagner apparently wrote:
>>    
>>>here is a code. However it doesn't make use of LAPACK. So
>>>it might be to slow.
>>>      
>>This exchange suggests to me that Robert thinks that
>>wrapping the LAPACK code is approximately trivial but that
>>Nils does not see it that way.  Apologies for any
>>misunderstanding.
>>    
>
>But in this particular case previous exposure to wrapping lapack
>functions is granted:
>http://aspn.activestate.com/ASPN/Mail/Message/scipy-user/2016895
>(for dgehrd)
>
>The number of arguments of DGEQP3 is the same as for dgehrd,
>so on this level it is not more difficult.
>
>  
>>Anyway, coming from GAUSS as a *user*, I would not know
>>where to begin if I needed to wrap a new LAPACK function.
>>    
>
>Thanks to f2py, wrapping Fortran code is (with a bit of effort)
>trivial in many cases. For complicated functions
>requiring many arguments the wrapper can become longish.
>Fortunately, many things can be learnt from
>looking at ``scipy/Lib/linalg/generic_flapack.pyf``
>In particular, the documentation at
>http://cens.ioc.ee/projects/f2py2e/
>is excellent.
>I also found  the f2py notes by Fernando Perez very helpful,
>http://cens.ioc.ee/pipermail/f2py-users/2003-April/000472.html
>
>Let me try to give some general remarks on how to start
>(the real authority on all this is of course Pearu, so
>please correct me if I got things wrong here):
>- first find a routine which will do the job you want:
>  - If the lapack documentation is installed properly
>    on Linux you could do
>      apropos  keyword
>  - www.netlib.org provides a nice decision tree
>- make sure that that it does not exist in scipy:
>    from scipy.lib import lapack
>    lapack.clapack.<TAB>         (assuming Ipython)
>    lapack.clapack.<routine_name>
>
>  Remark: routines starting with c/z are for double/single complex
>  and routines for d/s for double/single real numbers.
>  The calling sequence for c/z and d/s are (I think always) the same and
>  sometimes they are also the same for the real and complex case.
>- Then one has to download the fortran file for the lapack routine
>  of interest.
>- Generate wrapper by calling
>    f2py -m wrap_lap -h wrap_lap.pyf lapack_routine.f
>
>- Generate library
>    f2py  -c wrap_lap.pyf  lapack_routine.f -latlas -llapack -lblas
>
>- You can use this by
>
>    import wrap_lap
>
>  Note, that this is not yet polished (this is the part on
>  which has to spent some effort ;-), i.e. one
>  has to tell which variables are input, which are output
>  and which are optional. In addition temporary
>  storage has to be provided with the right dimensions
>  as described in the documentation part of the lapack routine.
>
>Concrete (and very simple) example (non-lapack):
>
>Wrapping Hermite polynomials
>----------------------------
>
>Download code (found after hours of googling ;-)::
>
>  http://cdm.unimo.it/home/matematica/funaro.daniele/splib.txt
>
>and extract ``hermite.f``
>
>Generate wrapper framework::
>
>  # only run the following line _once_
>  # (and never again, otherwise the hand-modified hermite.pyf
>  #  goes down the drains)
>  f2py -m hermite -h hermite.pyf hermite.f
>
>Then modify ``hermite.pyf``
>
>Create the module::
>
>  f2py2.3  -c hermite.pyf  hermite.f
>
>  # add this if you want:
>  -DF2PY_REPORT_ON_ARRAY_COPY=1 -DF2PY_REPORT_ATEXIT
>
>Simple test::
>
>  import hermite
>  hermite.vahepo(2,2.0)
>  import scipy
>  scipy.special.hermite(2)(2.0)
>
>A more complicated example about
>how to wrap routines for band matrices can be found at
>http://www.physik.tu-dresden.de/~baecker/comp_talks.html
>under "Python and Co - some recent developments".
>
>It would be nice if someone could add notes
>on all this to the scipy wiki.
>
>Best, Arnd
>
>_______________________________________________
>SciPy-user mailing list
>SciPy-user at scipy.net
>http://www.scipy.net/mailman/listinfo/scipy-user
>  
 Hi Arnd,

If it is so trivial to add the rank revealing QR decomposition to scipy,
I would be glad if the several steps on building a wrapper for DGEQP3
could be
documented on the Wiki as an exercise.
The best way to learn something is to have a well-elaborated example. :-)

Nils



More information about the SciPy-user mailing list