[Scipy-svn] r6989 - in trunk/scipy/sparse/linalg/eigen/arpack: . tests

scipy-svn@scip... scipy-svn@scip...
Sat Dec 4 14:25:25 CST 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