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

[lorenzo bolla]

hi all,
I've attached a patch for the problem here:
http://scipy.org/scipy/scipy/ticket/231
maybe, some of the developer will look into it and commit it to cvs, if
correct.
hth,
lorenzo.


On 8/23/07, lorenzo bolla <lbolla at gmail.com> wrote:
>
> 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 at 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://mail.scipy.org/pipermail/scipy-user/attachments/20070830/99067919/attachment.html>