[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