Mailman 3 November 2003 - SciPy-Dev

band matrices in scipy
by Arnd Baecker Jan. 5, 2004

Jan. 5, 2004

Hi, I am in the progress of wrapping the Lapack dsbev.f routine (for band-matrices). Thanks to the amazing f2py I got pretty far (actually in addition the good docs and Fernandos notes helped a lot!). Once this is finished I would be happy to see this included in scipy ;-). But before this could happen I have a couple of questions (sorry, the text below is long, but I wanted to be explicit ;-): - Is the following way of calling dsbev from python alright, (eg., I omitted ldab on the python side.) Would that parameter be useful to someone so that it is better to expose it (the full details of my few pyf lines are at the end of this mail)? In [7]: MODULENAME.dsbev? Type: fortran String Form: <fortran object at 0x4028ebf0> Namespace: Interactive Docstring: dsbev - Function signature: w,z,info = dsbev(ab,[compute_v,lower]) Required arguments: ab : input rank-2 array('d') with bounds (ldab,*) Optional arguments: compute_v := 1 input int lower := 0 input int Return objects: w : rank-1 array('d') with bounds (n) z : rank-2 array('d') with bounds (ldz,ldz) info : int - How can one supply better doc-strings to this routine? I.e., help(MODULENAME.dsbev) just gives "Help on fortran: <fortran object>" and scipy.info gives the same as ? from Ipython. I would like to add a more detailed description and an example. Is this somehow possible to do from within the .pyf wrapper? (I know that this is a pure f2py question, so maybe I should ask this question as well on the f2py mailing list...) - And one more technical aspect (before I forget about it): Is the wrapping done such that for _large_ matrices no unnecessary copies of arrays are made? - One thing which puzzled me a bit (presumably it is alright) concerns the return object z which contains the eigenvectors if compute_v=1. For compute_v=0 I changed ldz in the wrapper such that it is of length 1 (it is never referenced by dsbev.f). So if someone is just interested in eigenvalues and not eigenvectors a dummy variable for z has to be given. Is this alright (for a more "low" level routine like this)? In principle one could write a wrapper a la scipy.linalg.eig for the band matrices, but for my purposes I don't mind about the (sometimes needed) dummy variable. - Is the following correct? To get this into scipy.linalg.flapack the lines from MODULENAME.pyf have to be added to generic_flapack.pyf and <tchar=s,d> and <type_in_c>* and <type_in> have to be added at the corresponding place. And that's all? - Unittests: I think it would be good to have a few unit-tests for these routines. Presently I just created a sample 10x10 band matrix (see below), converted this to upper band storage and compared the eigenvalues/eigenvectors of scipy.linalg.eig with those of dsbev. In principle one should do this for both upper and lower triangle storage and for single and double. Would something like that be a reasonable unit-test for inclusion in scipy or is there some better approach (better test-matrices etc.)? - Further routines I am thinking of wrapping in this context a) for symmetric band matrices: dsbevd, dsbevx b) for Hermitian band matrices: zhbev, zhbevd, zhbevx and their single precision partners. Any opinions? Many thanks, Arnd ############################# And here the full example plus test: MODULENAME.pyf: # ------ snip here ------------------- ! -*- f90 -*- python module MODULENAME ! in interface ! in :MODULENAME !subroutine dsbev(ab,compute_v,lower) ! in :MODULENAME:dsbev.f subroutine dsbev(ab,compute_v,lower,n,ldab,kd,w,z,ldz,work,info) ! in :MODULENAME:dsbev.f callstatement (*f2py_func)((compute_v?"V":"N"),(lower?"L":"U"),&n,&kd,ab,&ldab,w,z,&ldz,work,&info) callprotoargument char*,char*,int*,int*,double*,int*,double*,double*,int*,double*,int* double precision dimension(ldab,*) :: ab integer optional,intent(in):: compute_v = 1 check(compute_v==1||compute_v==0) compute_v integer optional,intent(in),check(lower==0||lower==1) :: lower = 0 ! FIXME: does one really need the ldab as optional argument ? !integer optional,check(shape(ab,0)==ldab),depend(ab) :: ldab=shape(ab,0) integer intent(hide),depend(ab) :: ldab=shape(ab,0) integer intent(hide),depend(ab) :: n=shape(ab,1) integer intent(hide),depend(ab) :: kd=shape(ab,0)-1 double precision dimension(n),intent(out),depend(n) :: w double precision dimension(ldz,ldz),intent(out),depend(ldz) :: z ! For compute_v=1 z is used and contains the eigenvectors integer intent(hide),depend(n) :: ldz=compute_v?n:1 double precision dimension(ldz,ldz),depend(ldz) :: z double precision dimension(MAX(1,3*n-1)),intent(hide),depend(n) :: work integer intent(out)::info end subroutine dsbev end interface end python module MODULENAME ! This file was auto-generated with f2py (version:2.37.233-1545). ! See http://cens.ioc.ee/projects/f2py2e/ ! and then modified ... # ------ snip here ------------------- Then (dsbev.f from lapack is needed for this) f2py --debug-capi -c MODULENAME.pyf dsbev.f -llapack and a small test is ########### snip here ########### import MODULENAME from Numeric import * from scipy import diag,linalg N=10 fullmat=diag(-1.0*ones(N-1),-1)+diag(-1.0*ones(N-1),1)+diag(1.0*ones(N)) KD=2 # number of diagonals above the diagonal LDAB=KD+1 mat=zeros( (LDAB,N),typecode=Float) for i in arange(LDAB): # Transform to band matrix mat[LDAB-i-1,i:N]=diag(fullmat,i) w,evec,info=MODULENAME.dsbev(mat) w_lin,evec_lin=linalg.eig(fullmat) print sort(w)-sort(w_lin.real) print evec[:,0] print evec_lin[:,0] ########### snip here ###########

3 3

Re: [SciPy-user] MatrixMarket I/O Functions for scipy
by Nils Wagner Nov. 28, 2003

Nov. 28, 2003

------------------- Hi Nils, On Thursday 27 Nov 2003 1:37 pm, Nils Wagner wrote: > Afaik, there exists only a Fortran routine for this purpose. > http://math.nist.gov/~KRemington/mmio/mmiof77.html Since the source code is available, it is (fairly) trivial to wrap the f77 code using f2py and g77. This will allow you to call said functions from within python (+scipy). May I recommend you read the f2py tutorial? It is to be found on <http://www.ifi.uio.no/in228/scripting/doc/python/fc/f2py-tutorial/ tutorial.pdf> Best luck José Hi Jose, A google search on Python MatrixMarket yields http://people.web.psi.ch/geus/pyfemax/pysparse.html It should not be stand-alone solution for myself. I guess, many users will benefit from the feature (I/O of matrices in MatrixMarket or Harwell-Boeing format), if it becomes a part of scipy. Any suggestion ? BTW, the Harwell-Boeing seems to be supported by scipy http://scipy.net/cgi-bin/viewcvsx.cgi/scipy1/sparse/SuperLU/SRC/ zreadhb.c?rev=1.2 How can I invoke this function in scipy ? A small example is appreciated. Nils -- José L Gómez Dans PhD student Tel: +44 114 222 5582 Radar & Communications Group FAX; +44 870 132 2990 Department of Electronic Engineering University of Sheffield UK

1 0

Segfaults with weave
by Travis Oliphant Nov. 27, 2003

Nov. 27, 2003

With Python 2.3.2 on Mandrake 9.2 I ran into trouble with segfaults and weave. I traced the problem to bsddb module in Python 2.3.2. Apparently it segfaults under my configuration when trying to open a file. If you are running Python 2.3.2 could you try import bsddb bsddb.hashopen('anyfile','c') and see if that works for you. I fixed the problem by remmoving _bsddb.so (so that anydbm doesn't use the bsddb module as a db). Now, weave works fine. Could somebody else veryify to see if bsddb is working on Python 2.3.2 for you? -Travis Oliphant

2 1

Machine precision
by Nils Wagner Nov. 27, 2003

Nov. 27, 2003

Dear experts, Is it possible to implement eps as a constant (similar to pi or e) representing the precision of the machine in scipy ? Any suggestion ? Nils

3 2

ATLAS version and failures inscipy.test()
by Nils Wagner Nov. 27, 2003

Nov. 27, 2003

Pearu, I switched back to ATLAS3.4.1 http://www.cs.utk.edu/~rwhaley/ATLAS/errata3.4.1.html#rLLt Now, scipy.test() gives ====================================================================== FAIL: check_normal (scipy.stats.morestats.test_morestats.test_anderson) ---------------------------------------------------------------------- Traceback (most recent call last): File "/usr/local/lib/python2.3/site-packages/scipy/stats/tests/ test_morestats.py", line 49, in check_normal assert(scipy.all(A < crit[-2:])) AssertionError ---------------------------------------------------------------------- Ran 1028 tests in 2.348s FAILED (failures=1) Nils

2 2

Testing SpamAssassin installation
by Joe Cooper Nov. 26, 2003

Nov. 26, 2003

Hi all, I've just installed SpamAssassin, to work alongside SpamBayes on the SciPy mailing lists. Though very few spams are getting through these days, because we don't permit posting by non-subscribers, it is still a burden to sift through the moderated requests to let the 5% legitimate messages through. Anyway, feel free to ignore this message, as it it just a test to be sure we're still delivering. Thanks.

1 0

FAIL:check_simple_complexFAIL:check_random_complex
by Nils Wagner Nov. 26, 2003

Nov. 26, 2003

Pearu, I have used ATLAS3.5.12 as well. However, the failures didn't vanish. So, what is the reason for these strange failures ? gcc ??? Numeric ??? scipy ??? f2py ??? LAPACK ??? ATLAS ??? >>> scipy.__version__ '0.2.1_251.4457' >>> Numeric.__version__ '23.1' nwagner@jane:~> f2py -v 2.37.233-1545 nwagner@jane:~> gcc -v Reading specs from /usr/lib/gcc-lib/i586-suse-linux/3.3.1/specs Konfiguriert mit: ../configure --enable-threads=posix --prefix=/usr --with-local-prefix=/usr/local --infodir=/usr/share/info --mandir=/usr/ share/man --libdir=/usr/lib --enable-languages=c,c++,f77,objc,java,ada --disable-checking --enable-libgcj --with-gxx-include-dir=/usr/include/g ++ --with-slibdir=/lib --with-system-zlib --enable-shared --enable-__cxa_atexit i586-suse-linux Thread model: posix gcc-Version 3.3.1 (SuSE Linux) Is it possible to neglect these failures ? ====================================================================== FAIL: check_random_complex (scipy.linalg.decomp.test_decomp.test_cholesky) ---------------------------------------------------------------------- Traceback (most recent call last): File "/usr/local/lib/python2.3/site-packages/scipy/linalg/tests/ test_decomp.py", line 276, in check_random_complex assert_array_almost_equal(cholesky(a),c) File "/usr/local/lib/python2.3/site-packages/scipy_test/testing.py", line 592, in assert_array_almost_equal assert alltrue(ravel(reduced)),\ AssertionError: Arrays are not almost equal: ====================================================================== FAIL: check_simple_complex (scipy.linalg.decomp.test_decomp.test_cholesky) ---------------------------------------------------------------------- Traceback (most recent call last): File "/usr/local/lib/python2.3/site-packages/scipy/linalg/tests/ test_decomp.py", line 250, in check_simple_complex assert_array_almost_equal(cholesky(a),c) File "/usr/local/lib/python2.3/site-packages/scipy_test/testing.py", line 592, in assert_array_almost_equal assert alltrue(ravel(reduced)),\ AssertionError: Arrays are not almost equal: Nils

2 1

Reason for failures inscipy.test() ?
by Nils Wagner Nov. 25, 2003

Nov. 25, 2003

Dear experts, Can someone reproduce these failures with latest cvs ? ====================================================================== FAIL: check_random_complex (scipy.linalg.decomp.test_decomp.test_cholesky) ---------------------------------------------------------------------- Traceback (most recent call last): File "/usr/local/lib/python2.3/site-packages/scipy/linalg/tests/ test_decomp.py", line 276, in check_random_complex assert_array_almost_equal(cholesky(a),c) File "/usr/local/lib/python2.3/site-packages/scipy_test/testing.py", line 592, in assert_array_almost_equal assert alltrue(ravel(reduced)),\ AssertionError: Arrays are not almost equal: ====================================================================== FAIL: check_simple_complex (scipy.linalg.decomp.test_decomp.test_cholesky) ---------------------------------------------------------------------- Traceback (most recent call last): File "/usr/local/lib/python2.3/site-packages/scipy/linalg/tests/ test_decomp.py", line 250, in check_simple_complex assert_array_almost_equal(cholesky(a),c) File "/usr/local/lib/python2.3/site-packages/scipy_test/testing.py", line 592, in assert_array_almost_equal assert alltrue(ravel(reduced)),\ AssertionError: Arrays are not almost equal: ---------------------------------------------------------------------- Ran 1028 tests in 2.302s FAILED (failures=2) <unittest.TextTestRunner object at 0x41977bcc> >>> I guess these failures are related to ATLAS or gcc, since I was able to run scipy.test() on another system without any failures. They differ in ATLAS and gcc. Any idea or hint ? Nils ATLAS version 3.4.2 built by root on Fre Sep 12 22:56:10 UTC 2003: UNAME : Linux linux 2.4.20-4GB-athlon #1 Mon Mar 17 17:56:47 UTC 2003 i686 unknown unknown GNU/Linux INSTFLG : MMDEF : /home/nwagner/mysoftware/ATLAS/CONFIG/ARCHS/ ATHLONSSE1/gcc/gemm ARCHDEF : /home/nwagner/mysoftware/ATLAS/CONFIG/ARCHS/ ATHLONSSE1/gcc/misc F2CDEFS : -DAdd__ -DStringSunStyle CACHEEDGE: 131072 F77 : /usr/bin/g77, version GNU Fortran (GCC) 3.3 20030226 (prerelease) (SuSE Linux) F77FLAGS : -O CC : /usr/bin/gcc, version gcc (GCC) 3.3 20030226 (prerelease) (SuSE Linux) CC FLAGS : -fomit-frame-pointer -O3 -funroll-all-loops MCC : /usr/bin/gcc, version gcc (GCC) 3.3 20030226 (prerelease) (SuSE Linux) MCCFLAGS : -fomit-frame-pointer -O nwagner@jane:~/cvs/scipy/Lib/linalg>

4 6

Problems with test_fitpack.py fitpack2.py
by Nils Wagner Nov. 25, 2003

Nov. 25, 2003

Dear experts, What is the reason for the following errors ? ===================================================================== ERROR: check_linear_constant (scipy.interpolate.fitpack.test_fitpack.test_LSQBivariateSpline) ---------------------------------------------------------------------- Traceback (most recent call last): File "/usr/local/lib/python2.1/site-packages/scipy/interpolate/tests/test_fitpack.py", line 52, in check_linear_constant lut = LSQBivariateSpline(x,y,z,tx,ty,kx=1,ky=1) File "/usr/local/lib/python2.1/site-packages/scipy/interpolate/fitpack2.py", line 399, in __init__ kx,ky,eps,lwrk2=1) error: (len(y)==m) failed for 3rd argument z ====================================================================== ERROR: check_linear_1d (scipy.interpolate.fitpack.test_fitpack.test_SmoothBivariateSpline) ---------------------------------------------------------------------- Traceback (most recent call last): File "/usr/local/lib/python2.1/site-packages/scipy/interpolate/tests/test_fitpack.py", line 72, in check_linear_1d lut = SmoothBivariateSpline(x,y,z,kx=1,ky=1) File "/usr/local/lib/python2.1/site-packages/scipy/interpolate/fitpack2.py", line 356, in __init__ kx,ky,s=s,eps=eps,lwrk2=1) error: (len(y)==m) failed for 3rd argument z ====================================================================== ERROR: check_linear_constant (scipy.interpolate.fitpack.test_fitpack.test_SmoothBivariateSpline) ---------------------------------------------------------------------- Traceback (most recent call last): File "/usr/local/lib/python2.1/site-packages/scipy/interpolate/tests/test_fitpack.py", line 62, in check_linear_constant lut = SmoothBivariateSpline(x,y,z,kx=1,ky=1) File "/usr/local/lib/python2.1/site-packages/scipy/interpolate/fitpack2.py", line 356, in __init__ kx,ky,s=s,eps=eps,lwrk2=1) error: (len(y)==m) failed for 3rd argument z Nils

4 8

FAILED (failures=2) in scipy.test() with Python2.2
by Nils Wagner Nov. 20, 2003

Nov. 20, 2003

Dear experts, scipy.test() with Python2.2 results in >>> scipy.__version__ '0.2.1_249.4443' >>> Numeric.__version__ '23.1' >>> ====================================================================== FAIL: check_pro_ang1_cv (scipy.special.basic.test_basic.test_cephes) ---------------------------------------------------------------------- Traceback (most recent call last): File "/usr/lib/python2.2/site-packages/scipy/special/tests/test_basic.py", line 384, in check_pro_ang1_cv assert_equal(cephes.pro_ang1_cv(1,1,1,1,0),(1.0,0.0)) File "/usr/lib/python2.2/site-packages/scipy_test/testing.py", line 519, in assert_equal assert desired == actual, msg AssertionError: Items are not equal: DESIRED: (1.0, 0.0) ACTUAL: (0.99999999999999989, 0.0) ====================================================================== FAIL: check_pro_ang1_cv (scipy.special.cephes.test_cephes.test_cephes) ---------------------------------------------------------------------- Traceback (most recent call last): File "/usr/lib/python2.2/site-packages/scipy/special/tests/test_cephes.py", line 355, in check_pro_ang1_cv assert_equal(cephes.pro_ang1_cv(1,1,1,1,0),(1.0,0.0)) File "/usr/lib/python2.2/site-packages/scipy_test/testing.py", line 519, in assert_equal assert desired == actual, msg AssertionError: Items are not equal: DESIRED: (1.0, 0.0) ACTUAL: (0.99999999999999989, 0.0) ---------------------------------------------------------------------- Ran 1195 tests in 3.637s FAILED (failures=2) <unittest.TextTestRunner instance at 0x88fb05c> >>> Nils

1 0