[Scipy-svn] r6272 - trunk/scipy/sparse/linalg/eigen/arpack
scipy-svn at scipy.org
scipy-svn at scipy.org
Fri Mar 26 01:35:26 EDT 2010
Author: cdavid
Date: 2010-03-26 00:35:25 -0500 (Fri, 26 Mar 2010)
New Revision: 6272
Modified:
trunk/scipy/sparse/linalg/eigen/arpack/arpack.py
Log:
REF: abstract eigen value/vector extraction into *Params classes.
Modified: trunk/scipy/sparse/linalg/eigen/arpack/arpack.py
===================================================================
--- trunk/scipy/sparse/linalg/eigen/arpack/arpack.py 2010-03-26 05:35:16 UTC (rev 6271)
+++ trunk/scipy/sparse/linalg/eigen/arpack/arpack.py 2010-03-26 05:35:25 UTC (rev 6272)
@@ -123,7 +123,7 @@
ltr = _type_conv[self.tp]
self._arpack_solver = _arpack.__dict__[ltr + 'saupd']
- self.extract = _arpack.__dict__[ltr + 'seupd']
+ self._arpack_extract = _arpack.__dict__[ltr + 'seupd']
self.ipntr = np.zeros(11, "int")
@@ -152,6 +152,26 @@
if self.iparam[4] < self.k:
warnings.warn("Only %d/%d eigenvectors converged" % (self.iparam[4], self.k))
+ def extract(self, return_eigenvectors):
+ rvec = return_eigenvectors
+ ierr = 0
+ howmny = 'A' # return all eigenvectors
+ sselect = np.zeros(self.ncv, 'int') # unused
+ sigma = 0.0 # no shifts, not implemented
+
+ d, z, info = self._arpack_extract(rvec, howmny, sselect, sigma, self.bmat,
+ self.which, self.k, self.tol, self.resid, self.v,
+ self.iparam[0:7], self.ipntr, self.workd[0:2*self.n],
+ self.workl,ierr)
+
+ if ierr != 0:
+ raise RuntimeError("Error info=%d in arpack" % params.info)
+
+ if return_eigenvectors:
+ return d, z
+ else:
+ return d
+
class _UnsymmetricArpackParams(_ArpackParams):
def __init__(self, n, k, tp, matvec, sigma=None,
ncv=None, v0=None, maxiter=None, which="LM", tol=0):
@@ -166,7 +186,7 @@
ltr = _type_conv[self.tp]
self._arpack_solver = _arpack.__dict__[ltr + 'naupd']
- self.extract = _arpack.__dict__[ltr + 'neupd']
+ self._arpack_extract = _arpack.__dict__[ltr + 'neupd']
self.ipntr = np.zeros(14, "int")
@@ -203,6 +223,88 @@
elif self.info == -1:
warnings.warn("Maximum number of iterations taken: %s" % self.iparam[2])
+ def extract(self, return_eigenvectors):
+ k, n = self.k, self.n
+
+ ierr = 0
+ howmny = 'A' # return all eigenvectors
+ sselect = np.zeros(self.ncv, 'int') # unused
+ sigmai = 0.0 # no shifts, not implemented
+ sigmar = 0.0 # no shifts, not implemented
+ workev = np.zeros(3 * self.ncv, self.tp)
+
+ if self.tp in 'fd':
+ dr = np.zeros(k+1, self.tp)
+ di = np.zeros(k+1, self.tp)
+ zr = np.zeros((n, k+1), self.tp)
+ dr, di, zr, self.info=\
+ self._arpack_extract(return_eigenvectors,
+ howmny, sselect, sigmar, sigmai, workev,
+ self.bmat, self.which, k, self.tol, self.resid,
+ self.v, self.iparam, self.ipntr,
+ self.workd, self.workl, self.info)
+
+ # The ARPACK nonsymmetric real and double interface (s,d)naupd return
+ # eigenvalues and eigenvectors in real (float,double) arrays.
+
+ # Build complex eigenvalues from real and imaginary parts
+ d = dr + 1.0j * di
+
+ # Arrange the eigenvectors: complex eigenvectors are stored as
+ # real,imaginary in consecutive columns
+ z = zr.astype(self.tp.upper())
+ eps = np.finfo(self.tp).eps
+ i = 0
+ while i<=k:
+ # check if complex
+ if abs(d[i].imag) > eps:
+ # assume this is a complex conjugate pair with eigenvalues
+ # in consecutive columns
+ z[:,i] = zr[:,i] + 1.0j * zr[:,i+1]
+ z[:,i+1] = z[:,i].conjugate()
+ i +=1
+ i += 1
+
+ # Now we have k+1 possible eigenvalues and eigenvectors
+ # Return the ones specified by the keyword "which"
+ nreturned = self.iparam[4] # number of good eigenvalues returned
+ if nreturned == k: # we got exactly how many eigenvalues we wanted
+ d = d[:k]
+ z = z[:,:k]
+ else: # we got one extra eigenvalue (likely a cc pair, but which?)
+ # cut at approx precision for sorting
+ rd = np.round(d, decimals = _ndigits[self.tp])
+ if self.which in ['LR','SR']:
+ ind = np.argsort(rd.real)
+ elif self.which in ['LI','SI']:
+ # for LI,SI ARPACK returns largest,smallest abs(imaginary) why?
+ ind = np.argsort(abs(rd.imag))
+ else:
+ ind = np.argsort(abs(rd))
+ if self.which in ['LR','LM','LI']:
+ d = d[ind[-k:]]
+ z = z[:,ind[-k:]]
+ if self.which in ['SR','SM','SI']:
+ d = d[ind[:k]]
+ z = z[:,ind[:k]]
+
+ else:
+ # complex is so much simpler...
+ d, z, self.info =\
+ self._arpack_extract(return_eigenvectors,
+ howmny, sselect, sigmar, workev,
+ self.bmat, self.which, k, self.tol, self.resid,
+ self.v, self.iparam, self.ipntr,
+ self.workd, self.workl, self.rwork, ierr)
+
+ if ierr != 0:
+ raise RuntimeError("Error info=%d in arpack" % info)
+
+ if return_eigenvectors:
+ return d, z
+ else:
+ return d
+
def eigen(A, k=6, M=None, sigma=None, which='LM', v0=None,
ncv=None, maxiter=None, tol=0,
return_eigenvectors=True):
@@ -294,88 +396,8 @@
while not params.converged:
params.iterate()
- # now extract eigenvalues and (optionally) eigenvectors
- rvec = return_eigenvectors
- ierr = 0
- howmny = 'A' # return all eigenvectors
- sselect = np.zeros(params.ncv, 'int') # unused
- sigmai = 0.0 # no shifts, not implemented
- sigmar = 0.0 # no shifts, not implemented
- workev = np.zeros(3 * params.ncv, params.tp)
+ return params.extract(return_eigenvectors)
- if params.tp in 'fd':
- dr = np.zeros(k+1, params.tp)
- di = np.zeros(k+1, params.tp)
- zr = np.zeros((n, k+1), params.tp)
- dr, di, zr, params.info=\
- params.extract(rvec, howmny, sselect, sigmar, sigmai, workev,
- params.bmat, params.which, k, params.tol, params.resid,
- params.v, params.iparam, params.ipntr,
- params.workd, params.workl, params.info)
-
- # The ARPACK nonsymmetric real and double interface (s,d)naupd return
- # eigenvalues and eigenvectors in real (float,double) arrays.
-
- # Build complex eigenvalues from real and imaginary parts
- d = dr + 1.0j * di
-
- # Arrange the eigenvectors: complex eigenvectors are stored as
- # real,imaginary in consecutive columns
- z = zr.astype(params.tp.upper())
- eps = np.finfo(params.tp).eps
- i = 0
- while i<=k:
- # check if complex
- if abs(d[i].imag) > eps:
- # assume this is a complex conjugate pair with eigenvalues
- # in consecutive columns
- z[:,i] = zr[:,i] + 1.0j * zr[:,i+1]
- z[:,i+1] = z[:,i].conjugate()
- i +=1
- i += 1
-
- # Now we have k+1 possible eigenvalues and eigenvectors
- # Return the ones specified by the keyword "which"
- nreturned = params.iparam[4] # number of good eigenvalues returned
- if nreturned == k: # we got exactly how many eigenvalues we wanted
- d = d[:k]
- z = z[:,:k]
- else: # we got one extra eigenvalue (likely a cc pair, but which?)
- # cut at approx precision for sorting
- rd = np.round(d, decimals = _ndigits[params.tp])
- if params.which in ['LR','SR']:
- ind = np.argsort(rd.real)
- elif which in ['LI','SI']:
- # for LI,SI ARPACK returns largest,smallest abs(imaginary) why?
- ind = np.argsort(abs(rd.imag))
- else:
- ind = np.argsort(abs(rd))
- if params.which in ['LR','LM','LI']:
- d = d[ind[-k:]]
- z = z[:,ind[-k:]]
- if params.which in ['SR','SM','SI']:
- d = d[ind[:k]]
- z = z[:,ind[:k]]
-
-
- else:
- # complex is so much simpler...
- d, z, params.info =\
- params.extract(rvec, howmny, sselect, sigmar, workev,
- params.bmat, params.which, k, params.tol, params.resid,
- params.v, params.iparam, params.ipntr,
- params.workd, params.workl, params.rwork, ierr)
-
-
-
- if ierr != 0:
- raise RuntimeError("Error info=%d in arpack"%info)
- return None
- if return_eigenvectors:
- return d,z
- return d
-
-
def eigen_symmetric(A, k=6, M=None, sigma=None, which='LM', v0=None,
ncv=None, maxiter=None, tol=0,
return_eigenvectors=True):
@@ -467,24 +489,8 @@
while not params.converged:
params.iterate()
- # now extract eigenvalues and (optionally) eigenvectors
- rvec = return_eigenvectors
- ierr = 0
- howmny = 'A' # return all eigenvectors
- sselect = np.zeros(params.ncv, 'int') # unused
- sigma = 0.0 # no shifts, not implemented
+ return params.extract(return_eigenvectors)
- d, z, info = params.extract(rvec, howmny, sselect, sigma, params.bmat,
- params.which, k, params.tol, params.resid, params.v,
- params.iparam[0:7], params.ipntr, params.workd[0:2*n],
- params.workl,ierr)
-
- if ierr != 0:
- raise RuntimeError("Error info=%d in arpack" % params.info)
- if return_eigenvectors:
- return d, z
- return d
-
def svd(A, k=6):
"""Compute a few singular values/vectors for a sparse matrix using ARPACK.
More information about the Scipy-svn
mailing list