[Scipy-svn] r6989 - in trunk/scipy/sparse/linalg/eigen/arpack: . tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Sat Dec 4 15:25:25 EST 2010
Author: ptvirtan
Date: 2010-12-04 14:25:25 -0600 (Sat, 04 Dec 2010)
New Revision: 6989
Removed:
trunk/scipy/sparse/linalg/eigen/arpack/speigs.py
trunk/scipy/sparse/linalg/eigen/arpack/tests/test_speigs.py
Modified:
trunk/scipy/sparse/linalg/eigen/arpack/__init__.py
Log:
DEP: sparse/arpack: remove a duplicate ARPACK interface
Modified: trunk/scipy/sparse/linalg/eigen/arpack/__init__.py
===================================================================
--- trunk/scipy/sparse/linalg/eigen/arpack/__init__.py 2010-12-04 20:25:14 UTC (rev 6988)
+++ trunk/scipy/sparse/linalg/eigen/arpack/__init__.py 2010-12-04 20:25:25 UTC (rev 6989)
@@ -1,3 +1,2 @@
from info import __doc__
from arpack import *
-import speigs
Deleted: trunk/scipy/sparse/linalg/eigen/arpack/speigs.py
===================================================================
--- trunk/scipy/sparse/linalg/eigen/arpack/speigs.py 2010-12-04 20:25:14 UTC (rev 6988)
+++ trunk/scipy/sparse/linalg/eigen/arpack/speigs.py 2010-12-04 20:25:25 UTC (rev 6989)
@@ -1,224 +0,0 @@
-import numpy as np
-import _arpack
-
-__all___=['ArpackException','ARPACK_eigs', 'ARPACK_gen_eigs']
-
-class ArpackException(RuntimeError):
- ARPACKErrors = { 0: """Normal exit.""",
- 3: """No shifts could be applied during a cycle of the
- Implicitly restarted Arnoldi iteration. One possibility
- is to increase the size of NCV relative to NEV.""",
- -1: """N must be positive.""",
- -2: """NEV must be positive.""",
- -3: """NCV-NEV >= 2 and less than or equal to N.""",
- -4: """The maximum number of Arnoldi update iteration
- must be greater than zero.""",
- -5: """WHICH must be one of 'LM', 'SM', 'LR', 'SR', 'LI', 'SI'""",
- -6: """BMAT must be one of 'I' or 'G'.""",
- -7: """Length of private work array is not sufficient.""",
- -8: """Error return from LAPACK eigenvalue calculation;""",
- -9: """Starting vector is zero.""",
- -10: """IPARAM(7) must be 1,2,3,4.""",
- -11: """IPARAM(7) = 1 and BMAT = 'G' are incompatable.""",
- -12: """IPARAM(1) must be equal to 0 or 1.""",
- -9999: """Could not build an Arnoldi factorization.
- IPARAM(5) returns the size of the current Arnoldi
- factorization.""",
- }
- def __init__(self, info):
- self.info = info
- def __str__(self):
- try: return self.ARPACKErrors[self.info]
- except KeyError: return "Unknown ARPACK error"
-
-def check_init(n, nev, ncv):
- assert(nev <= n-4) # ARPACK seems to cause a segfault otherwise
- if ncv is None:
- ncv = min(2*nev+1, n-1)
- maxitr = max(n, 1000) # Maximum number of iterations
- return ncv, maxitr
-
-def init_workspaces(n,nev,ncv):
- ipntr = np.zeros(14, np.int32) # Pointers into memory structure used by F77 calls
- d = np.zeros((ncv, 3), np.float64, order='FORTRAN') # Temp workspace
- # Temp workspace/error residuals upon iteration completion
- resid = np.zeros(n, np.float64)
- workd = np.zeros(3*n, np.float64) # workspace
- workl = np.zeros(3*ncv*ncv+6*ncv, np.float64) # more workspace
- # Storage for the Arnoldi basis vectors
- v = np.zeros((n, ncv), dtype=np.float64, order='FORTRAN')
- return (ipntr, d, resid, workd, workl, v)
-
-def init_debug():
- # Causes various debug info to be printed by ARPACK
- _arpack.debug.ndigit = -3
- _arpack.debug.logfil = 6
- _arpack.debug.mnaitr = 0
- _arpack.debug.mnapps = 0
- _arpack.debug.mnaupd = 1
- _arpack.debug.mnaup2 = 0
- _arpack.debug.mneigh = 0
- _arpack.debug.mneupd = 1
-
-def init_postproc_workspace(n, nev, ncv):
- # Used as workspace and to return eigenvectors if requested. Not touched if
- # eigenvectors are not requested
- workev = np.zeros(3*ncv, np.float64) # More workspace
- select = np.zeros(ncv, np.int32) # Used as internal workspace since dneupd
- # parameter HOWMNY == 'A'
- return (workev, select)
-
-def postproc(n, nev, ncv, sigmar, sigmai, bmat, which,
- tol, resid, v, iparam, ipntr, workd, workl, info):
- workev, select = init_postproc_workspace(n, nev, ncv)
- ierr = 0
- # Postprocess the Arnouldi vectors to extract eigenvalues/vectors
- # If dneupd's first paramter is 'True' the eigenvectors are also calculated,
- # 'False' only the eigenvalues
- dr,di,z,info = _arpack.dneupd(
- True, 'A', select, sigmar, sigmai, workev, bmat, which, nev, tol, resid, v,
- iparam, ipntr, workd, workl, info)
-
- if np.abs(di[:-1]).max() == 0: dr = dr[:-1]
- else: dr = dr[:-1] + 1j*di[:-1]
- return (dr, z[:,:-1])
-
-
-def ARPACK_eigs(matvec, n, nev, which='SM', ncv=None, tol=1e-14):
- """
- Calculate eigenvalues for system with matrix-vector product matvec, dimension n
-
- Arguments
- =========
- matvec -- Function that provides matrix-vector product, i.e. matvec(x) -> A*x
- n -- Matrix dimension of the problem
- nev -- Number of eigenvalues to calculate
- which -- Spectrum selection. See details below. Defaults to 'SM'
- ncv -- Number of Arnoldi basisvectors to use. If None, default to 2*nev+1
- tol -- Numerical tollerance for Arnouldi iteration convergence. Defaults to 1e-14
-
- Spectrum Selection
- ==================
- which can take one of several values:
-
- 'LM' -> Request eigenvalues with largest magnitude.
- 'SM' -> Request eigenvalues with smallest magnitude.
- 'LR' -> Request eigenvalues with largest real part.
- 'SR' -> Request eigenvalues with smallest real part.
- 'LI' -> Request eigenvalues with largest imaginary part.
- 'SI' -> Request eigenvalues with smallest imaginary part.
-
- Return Values
- =============
- (eig_vals, eig_vecs) where eig_vals are the requested eigenvalues and
- eig_vecs the corresponding eigenvectors. If all the eigenvalues are real,
- eig_vals is a real array but if some eigenvalues are complex it is a
- complex array.
-
- """
- bmat = 'I' # Standard eigenproblem
- ncv, resid, iparam, ipntr, v, workd, workl, info = ARPACK_iteration(
- matvec, lambda x: x, n, bmat, which, nev, tol, ncv, mode=1)
- return postproc(n, nev, ncv, 0., 0., bmat, which, tol,
- resid, v, iparam, ipntr, workd, workl, info)
-
-def ARPACK_gen_eigs(matvec, sigma_solve, n, sigma, nev, which='LR', ncv=None, tol=1e-14):
- """
- Calculate eigenvalues close to sigma for generalised eigen system
-
- Given a system [A]x = k_i*[M]x where [A] and [M] are matrices and k_i are
- eigenvalues, nev eigenvalues close to sigma are calculated. The user needs
- to provide routines that calculate [M]*x and solve [A]-sigma*[M]*x = b for x.
-
- Arguments
- =========
- matvec -- Function that provides matrix-vector product, i.e. matvec(x) -> [M]*x
- sigma_solve -- sigma_solve(b) -> x, where [A]-sigma*[M]*x = b
- n -- Matrix dimension of the problem
- sigma -- Eigenvalue spectral shift real value
- nev -- Number of eigenvalues to calculate
- which -- Spectrum selection. See details below. Defaults to 'LR'
- ncv -- Number of Arnoldi basisvectors to use. If None, default to 2*nev+1
- tol -- Numerical tollerance for Arnouldi iteration convergence. Defaults to 1e-14
-
- Spectrum Shift
- ==============
-
- The spectrum of the orignal system is shifted by sigma. This transforms the
- original eigenvalues to be 1/(original_eig-sigma) in the shifted
- system. ARPACK then operates on the shifted system, transforming it back to
- the original system in a postprocessing step.
-
- The spectrum shift causes eigenvalues close to sigma to become very large
- in the transformed system. This allows quick convergence for these
- eigenvalues. This is particularly useful if a system has a number of
- trivial zero-eigenvalues that are to be ignored.
-
- Spectrum Selection
- ==================
- which can take one of several values:
-
- 'LM' -> Request spectrum shifted eigenvalues with largest magnitude.
- 'SM' -> Request spectrum shifted eigenvalues with smallest magnitude.
- 'LR' -> Request spectrum shifted eigenvalues with largest real part.
- 'SR' -> Request spectrum shifted eigenvalues with smallest real part.
- 'LI' -> Request spectrum shifted eigenvalues with largest imaginary part.
- 'SI' -> Request spectrum shifted eigenvalues with smallest imaginary part.
-
- The effect on the actual system is:
- 'LM' -> Eigenvalues closest to sigma on the complex plane
- 'LR' -> Eigenvalues with real part > sigma, provided they exist
-
-
- Return Values
- =============
- (eig_vals, eig_vecs) where eig_vals are the requested eigenvalues and
- eig_vecs the corresponding eigenvectors. If all the eigenvalues are real,
- eig_vals is a real array but if some eigenvalues are complex it is a
- complex array. The eigenvalues and vectors correspond to the original
- system, not the shifted system. The shifted system is only used interally.
-
- """
- bmat = 'G' # Generalised eigenproblem
- ncv, resid, iparam, ipntr, v, workd, workl, info = ARPACK_iteration(
- matvec, sigma_solve, n, bmat, which, nev, tol, ncv, mode=3)
- sigmar = sigma
- sigmai = 0.
- return postproc(n, nev, ncv, sigmar, sigmai, bmat, which, tol,
- resid, v, iparam, ipntr, workd, workl, info)
-
-def ARPACK_iteration(matvec, sigma_solve, n, bmat, which, nev, tol, ncv, mode):
- ncv, maxitr = check_init(n, nev, ncv)
- ipntr, d, resid, workd, workl, v = init_workspaces(n,nev,ncv)
- #init_debug()
- ishfts = 1 # Some random arpack parameter
- # Some random arpack parameter (I think it tells ARPACK to solve the
- # general eigenproblem using shift-invert
- iparam = np.zeros(11, np.int32) # Array with assorted extra paramters for F77 call
- iparam[[0,2,6]] = ishfts, maxitr, mode
- ido = 0 # Communication variable used by ARPACK to tell the user what to do
- info = 0 # Used for error reporting
- # Arnouldi iteration.
- while True:
- ido,resid,v,iparam,ipntr,info = _arpack.dnaupd(
- ido, bmat, which, nev, tol, resid, v, iparam, ipntr, workd, workl, info)
- if ido == -1 or ido == 1 and mode not in (3,4):
- # Perform y = inv[A - sigma*M]*M*x
- x = workd[ipntr[0]-1:ipntr[0]+n-1]
- Mx = matvec(x) # Mx = [M]*x
- workd[ipntr[1]-1:ipntr[1]+n-1] = sigma_solve(Mx)
- elif ido == 1: # Perform y = inv[A - sigma*M]*M*x using saved M*x
- # Mx = [M]*x where it was saved by ARPACK
- Mx = workd[ipntr[2]-1:ipntr[2]+n-1]
- workd[ipntr[1]-1:ipntr[1]+n-1] = sigma_solve(Mx)
- elif ido == 2: # Perform y = M*x
- x = workd[ipntr[0]-1:ipntr[0]+n-1]
- workd[ipntr[1]-1:ipntr[1]+n-1] = matvec(x)
- else: # Finished, or error
- break
- if info == 1:
- warn.warn("Maximum number of iterations taken: %s"%iparam[2])
- elif info != 0:
- raise ArpackException(info)
-
- return (ncv, resid, iparam, ipntr, v, workd, workl, info)
Deleted: trunk/scipy/sparse/linalg/eigen/arpack/tests/test_speigs.py
===================================================================
--- trunk/scipy/sparse/linalg/eigen/arpack/tests/test_speigs.py 2010-12-04 20:25:14 UTC (rev 6988)
+++ trunk/scipy/sparse/linalg/eigen/arpack/tests/test_speigs.py 2010-12-04 20:25:25 UTC (rev 6989)
@@ -1,52 +0,0 @@
-#!/usr/bin/env python
-
-from numpy.testing import run_module_suite, TestCase, assert_array_almost_equal
-
-from scipy.sparse.linalg.interface import aslinearoperator
-from scipy.sparse.linalg.eigen.arpack import speigs
-
-import numpy as np
-
-class TestEigs(TestCase):
- def test(self):
- maxn=15 # Dimension of square matrix to be solved
- # Use a PDP^-1 factorisation to construct matrix with known
- # eigenvalues/vectors. Used random eigenvectors initially.
- P = np.mat(np.random.random((maxn,)*2))
- P /= map(np.linalg.norm, P.T) # Normalise the eigenvectors
- D = np.mat(np.zeros((maxn,)*2))
- D[range(maxn), range(maxn)] = (np.arange(maxn, dtype=float)+1)/np.sqrt(maxn)
- A = P*D*np.linalg.inv(P)
- vals = np.array(D.diagonal())[0]
- vecs = P
- uv_sortind = vals.argsort()
- vals = vals[uv_sortind]
- vecs = vecs[:,uv_sortind]
-
- A=aslinearoperator(A)
- matvec = A.matvec
- #= lambda x: np.asarray(A*x)[0]
- nev=4
- eigvs = speigs.ARPACK_eigs(matvec, A.shape[0], nev=nev)
- calc_vals = eigvs[0]
- # Ensure the calculated eigenvectors have the same sign as the reference values
- calc_vecs = eigvs[1] / [np.sign(x[0]) for x in eigvs[1].T]
- assert_array_almost_equal(calc_vals, vals[0:nev], decimal=7)
- assert_array_almost_equal(calc_vecs, np.array(vecs)[:,0:nev], decimal=7)
-
-
-# class TestGeneigs(TestCase):
-# def test(self):
-# import pickle
-# from scipy.sparse.linalg import dsolve
-# A,B = pickle.load(file('mats.pickle'))
-# sigma = 27.
-# sigma_solve = dsolve.splu(A - sigma*B).solve
-# w = ARPACK_gen_eigs(B.matvec, sigma_solve, B.shape[0], sigma, 10)[0]
-# assert_array_almost_equal(w,
-# [27.346442255386375, 49.100299170945405, 56.508474856551544, 56.835800191692492,
-# 65.944215785041365, 66.194792400328367, 78.003788872725238, 79.550811647295944,
-# 94.646308846854879, 95.30841709116271], decimal=11)
-
-if __name__ == "__main__":
- run_module_suite()
More information about the Scipy-svn
mailing list