[SciPy-user] bug in ARPACK from scipy.sandbox?
lorenzo bolla
lbolla at gmail.com
Thu Aug 30 03:53:41 EDT 2007
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>
More information about the SciPy-User
mailing list