[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']:



More information about the Scipy-svn mailing list