[SciPy-dev] band matrices in scipy
Arnd Baecker
arnd.baecker at web.de
Mon Jan 5 03:51:43 CST 2004
Hi,
I continued working on the wrapper for dsbev, dsbevd and dsbevx.
The present status can be obtained from
http://www.physik.tu-dresden.de/~baecker/python/Band.zip
Before I continue with the last changes so that it could
be included in scipy it would be nice if someone (Pearu ?;-)
with more knowledge of f2py has a quick look at it.
Concerning dsbevx I have the following questions:
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...
And one more remark/question concerning the
copying of arrays when using f2py
- sorry if I re-state things already said before, but
I just want to make sure that I understand things correctly here.
First the option -DF2PY_REPORT_ON_ARRAY_COPY=1 is very helpful
(thanks Pearu for pointing this out).
With this I get the following picture (please correct me if I am wrong!):
If one allocates a matrix in python with
mat=zeros( (N,M))
this is not Fortran contiguous, so a copy is done.
However, declaring
mat=transpose(zeros(M,N))
one obtains a Fortran contiguous array and no copy is done.
If this is correct, would you say that
the second variant is a good way or is there a better one?
Many thanks,
Arnd
More information about the Scipy-dev
mailing list