[SciPy-user] bug in ARPACK from scipy.sandbox?

lorenzo bolla lbolla@gmail....
Thu Aug 23 03:42:19 CDT 2007


sorry for the noise, but I think I've found the bug...
this is what I changed in arpack.py to get the correct results (see the test
files attached).

should we commit the change to the CVS?

--------------------------------------------------------------

$> diff -c arpack_original.py
/usr/local/lib/python2.5/site-packages/scipy/sandbox/arpack/arpack.py

*** arpack_original.py  Thu Aug 23 10:30:44 2007
---
/usr/local/lib/python2.5/site-packages/scipy/sandbox/arpack/arpack.py
Thu Aug 23 10:39:59 2007
***************
*** 201,216 ****
          dr=sb.zeros(k+1,typ)
          di=sb.zeros(k+1,typ)
          zr=sb.zeros((n,k+1),typ)
!         dr,di,z,info=\
              eigextract(rvec,howmny,sselect,sigmar,sigmai,workev,
                     bmat,which,k,tol,resid,v,iparam,ipntr,
                     workd,workl,info)
!
          # make eigenvalues complex
!         d=dr+1.0j*di
          # futz with the eigenvectors:
          # complex are stored as real,imaginary in consecutive columns
!         z=zr.astype(typ.upper())
          for i in range(k): # fix c.c. pairs
              if di[i] > 0 :
                  z[:,i]=zr[:,i]+1.0j*zr[:,i+1]
--- 201,216 ----
          dr=sb.zeros(k+1,typ)
          di=sb.zeros(k+1,typ)
          zr=sb.zeros((n,k+1),typ)
!         dr,di,zr,info=\
              eigextract(rvec,howmny,sselect,sigmar,sigmai,workev,
                     bmat,which,k,tol,resid,v,iparam,ipntr,
                     workd,workl,info)
!
          # make eigenvalues complex
!         d=(dr+1.0j*di)[:k]
          # futz with the eigenvectors:
          # complex are stored as real,imaginary in consecutive columns
!         z=zr.astype(typ.upper())[:,:k]
          for i in range(k): # fix c.c. pairs
              if di[i] > 0 :
                  z[:,i]=zr[:,i]+1.0j*zr[:,i+1]

--------------------------------------------------------------


lorenzo.


On 8/23/07, lorenzo bolla <lbolla@gmail.com> wrote:
>
> I have problems with the ARPACK wrappers in scipy.sandbox.
> take a look at this snippet of code.
>
>
> ---------------------------------------------------------------------------------------------------------
>
> In [295]: A = numpy.array([[1,2,3,4],[0,2,3,4],[0,0,3,4],[0,0,0,4]],
> dtype=float)
>
> In [296]: A
> Out[296]:
> array([[ 1.,  2.,  3.,  4.],
>        [ 0.,  2.,  3.,  4.],
>        [ 0.,  0.,  3.,  4.],
>        [ 0.,  0.,  0.,  4.]])
>
> In [297]: [w,v] = arpack.eigen(A,2)
> In [298]: w
> Out[298]: array([ 4.+0.j,  3.+0.j,  0.+0.j])
>
> --> WRONG: I get 3 eigenvalues instead of two!
>
> In [299]: v
> Out[299]:
> array([[ 0.+0.j,  0.+0.j,  0.+0.j],
>        [ 0.+0.j,  0.+0.j,  0.+0.j],
>        [ 0.+0.j,  0.+0.j,  0.+0.j],
>        [ 0.+0.j,  0.+0.j,  0.+0.j]])
> --> WRONG: all the eigenvectors are null!
>
> In [300]: [w,v] = arpack.eigen(A.astype(numpy.complex),2)
>
> In [301]: w
> Out[301]: array([ 4. -2.41126563e-15j,  3. +1.34425147e-15j])
> --> CORRECT: casting the matrix to complex type gives the correct result
> and the correct numbers of eigenvalues
>
> In [302]: v
> Out[302]:
> array([[ -1.37221970e-01-0.75187207j,  -7.50019180e-01-0.32694452j],
>        [ -1.02916477e-01-0.56390405j,  -5.00012787e-01-0.21796301j],
>        [ -5.14582387e-02-0.28195202j,  -1.66670929e-01-0.07265434j ],
>        [ -1.28645597e-02-0.07048801j,   2.49800181e-16+0.j        ]])
> --> MAYBE: and the eigenvectors are not null, but...
>
> In [303]: [w,v] = arpack.eigen(A.astype(numpy.complex128),2)
>
> In [304]: w
> Out[304]: array([ 4. +7.28583860e-16j,  3. +2.23881966e-16j])
>
> In [305]: v
> Out[305]:
> array([[ -6.65958925e-01 -3.75020242e-01j,
>           8.08076904e-01 +1.28192062e-01j],
>        [ -4.99469194e-01 -2.81265182e-01j,
>           5.38717936e-01 +8.54613743e-02j],
>        [ - 2.49734597e-01 -1.40632591e-01j,
>           1.79572645e-01 +2.84871248e-02j],
>        [ -6.24336492e-02 -3.51581477e-02j,
>          -5.20417043e-16 -2.39391840e-16j]])
>  --> WRONG: casting to a complex128 changes the values of the
> eigenvectors!!!
>
>
> ---------------------------------------------------------------------------------------------------------
>
> in any case, the result for the eigenvectors are different than Matlab
> (while the eigenvalues are ok):
>
> v =
>
>    -8.181818181818171e-001    7.642914835078907e-001
>    -5.454545454545460e-001    5.732186126309180e-001
>    -1.818181818181832e-001    2.866093063154587e-001
>    -6.938893903907228e-018    7.165232657886386e-002
>
>
> w =
>
>     3.000000000000010e+000                         0
>                          0    3.999999999999995e+000
>
> Can someone explain me what's wrong?
> Thanks in advance,
> lorenzo.
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://projects.scipy.org/pipermail/scipy-user/attachments/20070823/ee2fd83e/attachment.html 
-------------- next part --------------
import numpy
from scipy import sparse
from scipy.sandbox import arpack

A = numpy.array([[1.,2,3,4],[0,2,3,4],[0,0,3,4],[0,0,0,4]])
[w,v] = arpack.eigen(A.astype('f'),2)
print w
print v

[w,v] = arpack.eigen(A.astype('d'),2)
print w
print v

[w,v] = arpack.eigen(A.astype('F'),2)
print w
print v
print numpy.abs(v)
print numpy.angle(v)

[w,v] = arpack.eigen(A.astype('D'),2)
print w
print v
print numpy.abs(v)
print numpy.angle(v)
-------------- next part --------------
[ 3.99998975+0.j  3.00000763+0.j]
[[  7.64292002e-01+0.j   8.18181396e-01+0.j]
 [  5.73218405e-01+0.j   5.45454800e-01+0.j]
 [  2.86608338e-01+0.j   1.81819022e-01+0.j]
 [  7.16515183e-02+0.j   4.24683094e-07+0.j]]
[ 4.+0.j  3.+0.j]
[[  7.64291484e-01+0.j   8.18181818e-01+0.j]
 [  5.73218613e-01+0.j   5.45454545e-01+0.j]
 [  2.86609306e-01+0.j   1.81818182e-01+0.j]
 [  7.16523266e-02+0.j  -1.11022302e-15+0.j]]
[ 4.00000143 -2.02525530e-06j  3.00000048 +3.33409253e-06j]
[[ -6.78678036e-01 +3.51479203e-01j   6.90856159e-01 -4.38337117e-01j]
 [ -5.09008467e-01 +2.63609469e-01j   4.60570931e-01 -2.92224437e-01j]
 [ -2.54504293e-01 +1.31804958e-01j   1.53523907e-01 -9.74078998e-02j]
 [ -6.36259913e-02 +3.29513997e-02j   3.53902578e-08 +1.07567757e-07j]]
[[  7.64291525e-01   8.18181932e-01]
 [  5.73218584e-01   5.45454562e-01]
 [  2.86609471e-01   1.81818292e-01]
 [  7.16523677e-02   1.13239977e-07]]
[[ 2.6637373  -0.56539017]
 [ 2.66373706 -0.56538951]
 [ 2.66373658 -0.56538761]
 [ 2.66373396  1.25294697]]
[ 4. +6.69343053e-15j  3. -9.84463195e-15j]
[[ -6.78176623e-01 +3.52445654e-01j   6.90719756e-01 -4.38551828e-01j]
 [ -5.08632467e-01 +2.64334241e-01j   4.60479838e-01 -2.92367885e-01j]
 [ -2.54316234e-01 +1.32167120e-01j   1.53493279e-01 -9.74559617e-02j]
 [ -6.35790584e-02 +3.30417801e-02j   9.28944421e-16 -1.13017234e-15j]]
[[  7.64291484e-01   8.18181818e-01]
 [  5.73218613e-01   5.45454545e-01]
 [  2.86609306e-01   1.81818182e-01]
 [  7.16523266e-02   1.46295156e-15]]
[[ 2.66231271 -0.56570106]
 [ 2.66231271 -0.56570106]
 [ 2.66231271 -0.56570106]
 [ 2.66231271 -0.88281419]]


More information about the SciPy-user mailing list