[Scipy-svn] r3902 - in trunk/scipy/splinalg/eigen/arpack: . tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Wed Feb 6 18:53:43 EST 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