[SciPy-dev] newscipy, band matrix routines
Arnd Baecker
arnd.baecker at web.de
Wed Nov 9 03:56:03 CST 2005
Hi,
a while ago I created f2py wrappers for
- dsbev
- dsbevd
- dgbtrf
- dgbtrs
- dgbrfs
which deal with eigenvalue/vector computation,
LU decomposition, ... for band matrices.
The present version of this, including examples/tests and
a little bit of documentation is at
http://www.physik.tu-dresden.de/~baecker/python/Band2.zip
I would like to get this into newscipy but we need some assistance:
- should the wrappers be added to generic_flapack.pyf ?
(+adding the new functions to setup_linalg.py so that the
corresponding single precision wrappers can be skipped)
- Are the only scipy specific tricks the use of
<tchar=s,d,c,z> and <type_in>
to support the different types?
- For testing/development purpuses, is there
a way to put this into a separate sandbox
or is it just easier to use all the tools which
are provided by Lib/linalg/?
- it would be great if an expert (Pearu - are you listening? ;-)
has a look at the wrapper (either already now, or
when the above steps are done) - see also the questions below.
Any further guidance/comments/suggestions are welcome!
Best, Arnd
Some more technical questions
Concerning dsbevx:
a) abstol parameter:
from dsbev.f:
"Eigenvalues will be computed most accurately when ABSTOL is
set to twice the underflow threshold 2*DLAMCH('S'), not zero."
Maybe the easiest (best?) is to wrap DLAMCH as well (easy)
and let the user provide the value, if necessary?
((Even nicer would be to compute this value internally -
is this possible with f2py ? If not, I am not sure if adding this
is worth the effort...))
b) Presently I suppress the array q
Q (output) DOUBLE PRECISION array, dimension (LDQ, N)
If JOBZ = 'V', the N-by-N orthogonal matrix used in the
reduction to tridiagonal form.
If JOBZ = 'N', the array Q is not referenced.
Should one better return this as well and therefore replace
double precision dimension(ldq,ldq),intent(hide),depend(ldq) :: q
by
double precision dimension(ldq,ldq),intent(out),depend(ldq) :: q
?
c) Should one make vl,vu,il,iu, optional as well ?
Background: dsbevx allows to compute
range: 0 = 'A': all eigenvalues will be found;
1 = 'V': all eigenvalues in the half-open interval (vl,vu]
will be found;
2 = 'I': the il-th through iu-th eigenvalues will be found.
so we may have that either vl,vu or il,iu or neither of
the two pairs are used.
If one makes all vl,vu,il,iu, optional, one could set
set the defaults to
il=0 iu=10
vl=0.0 vu=20.0 # this of course absolutely
# dependent on the problem,
# i.e in general nonsense...
More information about the Scipy-dev
mailing list