# [Scipy-svn] r3902 - in trunk/scipy/splinalg/eigen/arpack: . tests

scipy-svn@scip... scipy-svn@scip...
Wed Feb 6 17:53:43 CST 2008

```Author: hagberg
Date: 2008-02-06 17:53:38 -0600 (Wed, 06 Feb 2008)
New Revision: 3902

Modified:
trunk/scipy/splinalg/eigen/arpack/arpack.py
trunk/scipy/splinalg/eigen/arpack/tests/test_arpack.py
Log:
Allow starting vector for arpack eigensolver.  Clean up type handling.

Modified: trunk/scipy/splinalg/eigen/arpack/arpack.py
===================================================================
--- trunk/scipy/splinalg/eigen/arpack/arpack.py	2008-02-06 21:11:37 UTC (rev 3901)
+++ trunk/scipy/splinalg/eigen/arpack/arpack.py	2008-02-06 23:53:38 UTC (rev 3902)
@@ -51,7 +51,7 @@
_ndigits = {'f':5, 'd':12, 'F':5, 'D':12}

-def eigen(A, k=6, M=None, sigma=None, which='LM',
+def eigen(A, k=6, M=None, sigma=None, which='LM', v0=None,
ncv=None, maxiter=None, tol=0,
return_eigenvectors=True):
"""Find k eigenvalues and eigenvectors of the square matrix A.
@@ -91,6 +91,8 @@
(Not implemented)
Find eigenvalues near sigma.  Shift spectrum by sigma.

+    v0 : array
+        Starting vector for iteration.

ncv : integer
The number of Lanczos vectors generated
@@ -126,18 +128,17 @@

"""
A = aslinearoperator(A)
-
-    if M is not None:
-        raise NotImplementedError("generalized eigenproblem not supported yet")
-
if A.shape[0] != A.shape[1]:
raise ValueError('expected square matrix (shape=%s)' % shape)
-
n = A.shape[0]

+    # guess type
+    typ = A.dtype.char
+    if typ not in 'fdFD':
+        raise ValueError("matrix type must be 'f', 'd', 'F', or 'D'")
+
if M is not None:
raise NotImplementedError("generalized eigenproblem not supported yet")
-
if sigma is not None:
raise NotImplementedError("shifted eigenproblem not supported yet")

@@ -148,15 +149,14 @@
ncv=min(ncv,n)
if maxiter==None:
maxiter=n*10
+    # assign starting vector
+    if v0 is not None:
+        resid=v0
+        info=1
+    else:
+        resid = np.zeros(n,typ)
+        info=0

-    # guess type
-    resid = np.zeros(n,'f')
-    try:
-        typ = A.dtype.char
-    except AttributeError:
-        typ = A.matvec(resid).dtype.char
-    if typ not in 'fdFD':
-        raise ValueError("matrix type must be 'f', 'd', 'F', or 'D'")

# some sanity checks
if k <= 0:
@@ -175,20 +175,18 @@
ltr = _type_conv[typ]
eigsolver = _arpack.__dict__[ltr+'naupd']
eigextract = _arpack.__dict__[ltr+'neupd']
-    matvec = A.matvec

v = np.zeros((n,ncv),typ) # holds Ritz vectors
-    resid = np.zeros(n,typ) # residual
workd = np.zeros(3*n,typ) # workspace
workl = np.zeros(3*ncv*ncv+6*ncv,typ) # workspace
iparam = np.zeros(11,'int') # problem parameters
ipntr = np.zeros(14,'int') # pointers into workspaces
-    info = 0
ido = 0

if typ in 'FD':
rwork = np.zeros(ncv,typ.lower())

+    # set solver mode and parameters
# only supported mode is 1: Ax=lx
ishfts = 1
mode1 = 1
@@ -207,12 +205,15 @@
eigsolver(ido,bmat,which,k,tol,resid,v,iparam,ipntr,
workd,workl,rwork,info)

-        if (ido == -1 or ido == 1):
-            # compute y = A * x
-            xslice = slice(ipntr[0]-1, ipntr[0]-1+n)
-            yslice = slice(ipntr[1]-1, ipntr[1]-1+n)
-            workd[yslice]=matvec(workd[xslice])
-        else: # done
+        xslice = slice(ipntr[0]-1, ipntr[0]-1+n)
+        yslice = slice(ipntr[1]-1, ipntr[1]-1+n)
+        if ido == -1:
+            # initialization
+            workd[yslice]=A.matvec(workd[xslice])
+        elif ido == 1:
+            # compute y=Ax
+            workd[yslice]=A.matvec(workd[xslice])
+        else:
break

if  info < -1 :
@@ -300,8 +301,8 @@
return d

-def eigen_symmetric(A, k=6, M=None, sigma=None, which='LM',
-                    ncv=None, maxiter=None, tol=0, v0=None,
+def eigen_symmetric(A, k=6, M=None, sigma=None, which='LM', v0=None,
+                    ncv=None, maxiter=None, tol=0,
return_eigenvectors=True):
"""Find k eigenvalues and eigenvectors of the real symmetric
square matrix A.
@@ -341,6 +342,8 @@
(Not implemented)
Find eigenvalues near sigma.  Shift spectrum by sigma.

+    v0 : array
+        Starting vector for iteration.

ncv : integer
The number of Lanczos vectors generated
@@ -375,32 +378,33 @@
--------
"""
A = aslinearoperator(A)
-
if A.shape[0] != A.shape[1]:
raise ValueError('expected square matrix (shape=%s)' % shape)
-
n = A.shape[0]

+    # guess type
+    typ = A.dtype.char
+    if typ not in 'fd':
+        raise ValueError("matrix must be real valued (type must be 'f' or 'd')")
+
if M is not None:
raise NotImplementedError("generalized eigenproblem not supported yet")
if sigma is not None:
raise NotImplementedError("shifted eigenproblem not supported yet")
+
if ncv is None:
ncv=2*k+1
ncv=min(ncv,n)
if maxiter==None:
maxiter=n*10
+    # assign starting vector
+    if v0 is not None:
+        resid=v0
+        info=1
+    else:
+        resid = np.zeros(n,typ)
+        info=0

-
-    # guess type
-    resid = np.zeros(n,'f')
-    try:
-        typ = A.dtype.char
-    except AttributeError:
-        typ = A.matvec(resid).dtype.char
-    if typ not in 'fd':
-        raise ValueError("matrix must be real valued (type must be 'f' or 'd')")
-
# some sanity checks
if k <= 0:
raise ValueError("k must be positive, k=%d"%k)
@@ -418,17 +422,16 @@
ltr = _type_conv[typ]
eigsolver = _arpack.__dict__[ltr+'saupd']
eigextract = _arpack.__dict__[ltr+'seupd']
-    matvec = A.matvec

+    # set output arrays, parameters, and workspace
v = np.zeros((n,ncv),typ)
-    resid = np.zeros(n,typ)
workd = np.zeros(3*n,typ)
workl = np.zeros(ncv*(ncv+8),typ)
iparam = np.zeros(11,'int')
ipntr = np.zeros(11,'int')
-    info = 0
ido = 0

+    # set solver mode and parameters
# only supported mode is 1: Ax=lx
ishfts = 1
mode1 = 1
@@ -439,12 +442,17 @@

while True:
ido,resid,v,iparam,ipntr,info =\
-            eigsolver(ido,bmat,which,k,tol,resid,v,iparam,ipntr,
-               workd,workl,info)
-        if (ido == -1 or ido == 1):
-            xslice = slice(ipntr[0]-1, ipntr[0]-1+n)
-            yslice = slice(ipntr[1]-1, ipntr[1]-1+n)
-            workd[yslice]=matvec(workd[xslice])
+            eigsolver(ido,bmat,which,k,tol,resid,v,
+                      iparam,ipntr,workd,workl,info)
+
+        xslice = slice(ipntr[0]-1, ipntr[0]-1+n)
+        yslice = slice(ipntr[1]-1, ipntr[1]-1+n)
+        if ido == -1:
+            # initialization
+            workd[yslice]=A.matvec(workd[xslice])
+        elif ido == 1:
+            # compute y=Ax
+            workd[yslice]=A.matvec(workd[xslice])
else:
break

Modified: trunk/scipy/splinalg/eigen/arpack/tests/test_arpack.py
===================================================================
--- trunk/scipy/splinalg/eigen/arpack/tests/test_arpack.py	2008-02-06 21:11:37 UTC (rev 3901)
+++ trunk/scipy/splinalg/eigen/arpack/tests/test_arpack.py	2008-02-06 23:53:38 UTC (rev 3902)
@@ -8,7 +8,7 @@
from scipy.testing import *

from numpy import array,real,imag,finfo,concatenate,\
-    column_stack,argsort,dot,round,conj,sort
+    column_stack,argsort,dot,round,conj,sort,random
from scipy.splinalg.eigen.arpack import eigen_symmetric,eigen

@@ -89,10 +89,10 @@
high=range(len(eval))[-h:]
return eval[low+high]

-    def eval_evec(self,d,typ,k,which):
+    def eval_evec(self,d,typ,k,which,**kwds):
a=d['mat'].astype(typ)
exact_eval=self.get_exact_eval(d,typ,k,which)
-        eval,evec=eigen_symmetric(a,k,which=which)
+        eval,evec=eigen_symmetric(a,k,which=which,**kwds)
# check eigenvalues
assert_array_almost_equal(eval,exact_eval,decimal=_ndigits[typ])
# check eigenvectors A*evec=eval*evec
@@ -101,12 +101,20 @@
eval[i]*evec[:,i],
decimal=_ndigits[typ])

-    def test_symmetric(self):
+    def test_symmetric_modes(self):
k=2
for typ in 'fd':
for which in ['LM','SM','BE']:
self.eval_evec(self.symmetric[0],typ,k,which)

+    def test_starting_vector(self):
+        k=2
+        for typ in 'fd':
+            A=self.symmetric[0]['mat']
+            n=A.shape[0]
+            v0 = random.rand(n).astype(typ)
+            self.eval_evec(self.symmetric[0],typ,k,which='LM',v0=v0)
+

class TestEigenComplexSymmetric(TestArpack):

@@ -141,7 +149,7 @@
decimal=_ndigits[typ])

-#     def test_complex_symmetric(self):
+#     def test_complex_symmetric_modes(self):
#         k=2
#         for typ in 'FD':
#             for which in ['LM','SM','LR','SR']:
@@ -168,14 +176,14 @@
return ind[:k]

-    def eval_evec(self,d,typ,k,which):
+    def eval_evec(self,d,typ,k,which,**kwds):
a=d['mat'].astype(typ)
# get exact eigenvalues
exact_eval=d['eval'].astype(typ.upper())
ind=self.sort_choose(exact_eval,typ,k,which)
exact_eval=exact_eval[ind]
# compute eigenvalues
-        eval,evec=eigen(a,k,which=which)
+        eval,evec=eigen(a,k,which=which,**kwds)
ind=self.sort_choose(eval,typ,k,which)
eval=eval[ind]
evec=evec[:,ind]
@@ -188,7 +196,7 @@
decimal=_ndigits[typ])

-    def test_nonsymmetric(self):
+    def test_nonsymmetric_modes(self):
k=2
for typ in 'fd':
for which in ['LI','LR','LM','SM','SR','SI']:
@@ -197,8 +205,17 @@

+    def test_starting_vector(self):
+        k=2
+        for typ in 'fd':
+            A=self.symmetric[0]['mat']
+            n=A.shape[0]
+            v0 = random.rand(n).astype(typ)
+            self.eval_evec(self.symmetric[0],typ,k,which='LM',v0=v0)

+
+
class TestEigenComplexNonSymmetric(TestArpack):

def sort_choose(self,eval,typ,k,which):
@@ -241,7 +258,7 @@
decimal=_ndigits[typ])

-#     def test_complex_nonsymmetric(self):
+#     def test_complex_nonsymmetric_modes(self):
#         k=2
#         for typ in 'FD':
#             for which in ['LI','LR','LM','SI','SR','SM']:

```