[Scipy-svn] r3897 - in trunk/scipy: sandbox splinalg splinalg/eigen splinalg/eigen/arpack splinalg/eigen/arpack/tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Wed Feb 6 12:57:19 EST 2008
Author: hagberg
Date: 2008-02-06 11:56:43 -0600 (Wed, 06 Feb 2008)
New Revision: 3897
Added:
trunk/scipy/splinalg/eigen/arpack/
trunk/scipy/splinalg/eigen/arpack/ARPACK/
trunk/scipy/splinalg/eigen/arpack/README
trunk/scipy/splinalg/eigen/arpack/__init__.py
trunk/scipy/splinalg/eigen/arpack/arpack.py
trunk/scipy/splinalg/eigen/arpack/arpack.pyf.src
trunk/scipy/splinalg/eigen/arpack/info.py
trunk/scipy/splinalg/eigen/arpack/setup.py
trunk/scipy/splinalg/eigen/arpack/speigs.py
trunk/scipy/splinalg/eigen/arpack/tests/
Removed:
trunk/scipy/sandbox/arpack/
trunk/scipy/splinalg/eigen/arpack/ARPACK/
trunk/scipy/splinalg/eigen/arpack/README
trunk/scipy/splinalg/eigen/arpack/__init__.py
trunk/scipy/splinalg/eigen/arpack/arpack.py
trunk/scipy/splinalg/eigen/arpack/arpack.pyf.src
trunk/scipy/splinalg/eigen/arpack/info.py
trunk/scipy/splinalg/eigen/arpack/setup.py
trunk/scipy/splinalg/eigen/arpack/speigs.py
trunk/scipy/splinalg/eigen/arpack/tests/
Modified:
trunk/scipy/splinalg/__init__.py
trunk/scipy/splinalg/eigen/arpack/tests/test_arpack.py
trunk/scipy/splinalg/setup.py
Log:
Move arpack sparse eigenvalue solver from sandbox to splinalg.eigen
Also:
Fixes #554 and addresses #231 (the part related to the patch for returning
wrong number of eigenvalues).
Improves docstrings, code comments, implementation notes
Uses new LinearOperator
Modified: trunk/scipy/splinalg/__init__.py
===================================================================
--- trunk/scipy/splinalg/__init__.py 2008-02-06 01:18:26 UTC (rev 3896)
+++ trunk/scipy/splinalg/__init__.py 2008-02-06 17:56:43 UTC (rev 3897)
@@ -5,7 +5,9 @@
from isolve import *
from dsolve import *
from interface import *
+from eigen import *
+
__all__ = filter(lambda s:not s.startswith('_'),dir())
from scipy.testing.pkgtester import Tester
test = Tester().test
Copied: trunk/scipy/splinalg/eigen/arpack (from rev 3876, trunk/scipy/sandbox/arpack)
Copied: trunk/scipy/splinalg/eigen/arpack/ARPACK (from rev 3896, trunk/scipy/sandbox/arpack/ARPACK)
Deleted: trunk/scipy/splinalg/eigen/arpack/README
===================================================================
--- trunk/scipy/sandbox/arpack/README 2008-01-30 17:05:38 UTC (rev 3876)
+++ trunk/scipy/splinalg/eigen/arpack/README 2008-02-06 17:56:43 UTC (rev 3897)
@@ -1,98 +0,0 @@
-This is the ARPACK package from
-http://www.caam.rice.edu/software/ARPACK/
-
-Specifically the files are from
-http://www.caam.rice.edu/software/ARPACK/SRC/arpack96.tar.gz
-with the patch
-http://www.caam.rice.edu/software/ARPACK/SRC/patch.tar.gz
-
-The ARPACK README is at
-http://www.caam.rice.edu/software/ARPACK/SRC/readme.arpack
-
----
-
-ARPACK is a collection of Fortran77 subroutines designed to solve large
-scale eigenvalue problems.
-
-The package is designed to compute a few eigenvalues and corresponding
-eigenvectors of a general n by n matrix A. It is most appropriate for large
-sparse or structured matrices A where structured means that a matrix-vector
-product w <- Av requires order n rather than the usual order n**2 floating
-point operations. This software is based upon an algorithmic variant of the
-Arnoldi process called the Implicitly Restarted Arnoldi Method (IRAM). When
-the matrix A is symmetric it reduces to a variant of the Lanczos process
-called the Implicitly Restarted Lanczos Method (IRLM). These variants may be
-viewed as a synthesis of the Arnoldi/Lanczos process with the Implicitly
-Shifted QR technique that is suitable for large scale problems. For many
-standard problems, a matrix factorization is not required. Only the action
-of the matrix on a vector is needed. ARPACK software is capable of solving
-large scale symmetric, nonsymmetric, and generalized eigenproblems from
-significant application areas. The software is designed to compute a few (k)
-eigenvalues with user specified features such as those of largest real part
-or largest magnitude. Storage requirements are on the order of n*k locations.
-No auxiliary storage is required. A set of Schur basis vectors for the desired
-k-dimensional eigen-space is computed which is numerically orthogonal to working
-precision. Numerically accurate eigenvectors are available on request.
-
-Important Features:
-
- o Reverse Communication Interface.
- o Single and Double Precision Real Arithmetic Versions for Symmetric,
- Non-symmetric, Standard or Generalized Problems.
- o Single and Double Precision Complex Arithmetic Versions for Standard
- or Generalized Problems.
- o Routines for Banded Matrices - Standard or Generalized Problems.
- o Routines for The Singular Value Decomposition.
- o Example driver routines that may be used as templates to implement
- numerous Shift-Invert strategies for all problem types, data types
- and precision.
-
----
-
-The ARPACK license is BSD-like.
-http://www.caam.rice.edu/software/ARPACK/RiceBSD.doc
-
----
-
-Rice BSD Software License
-Permits source and binary redistribution of the software ARPACK and
-P_ARPACK for both non-commercial and commercial use.
-
- Copyright (©) 2001, Rice University
- Developed by D.C. Sorensen, R.B. Lehoucq, C. Yang, and K. Maschhoff.
- All rights reserved.
-
-
- Redistribution and use in source and binary forms, with or without
- modification, are permitted provided that the following conditions are
- met:
- . Redistributions of source code must retain the above copyright notice,
- this list of conditions and the following disclaimer.
- . Redistributions in binary form must reproduce the above copyright
- notice, this list of conditions and the following disclaimer in the
- documentation and/or other materials provided with the distribution.
- . If you modify the source for these routines we ask that you change the
- name of the routine and comment the changes made to the original.
- . Written notification is provided to the developers of intent to use
- this software. Also, we ask that use of ARPACK is properly cited in
- any resulting publications or software documentation.
- . Neither the name of Rice University (RICE) nor the names of its
- contributors may be used to endorse or promote products derived from
- this software without specific prior written permission.
-
-
-THIS SOFTWARE IS PROVIDED BY RICE AND CONTRIBUTORS "AS IS" AND ANY EXPRESS
-OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
-WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
-DISCLAIMED. IN NO EVENT SHALL RICE OR CONTRIBUTORS BE LIABLE FOR ANY
-DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
-(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
-SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
-CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
-LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY
-OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
-DAMAGE.
-
-
-
-
Copied: trunk/scipy/splinalg/eigen/arpack/README (from rev 3896, trunk/scipy/sandbox/arpack/README)
Deleted: trunk/scipy/splinalg/eigen/arpack/__init__.py
===================================================================
--- trunk/scipy/sandbox/arpack/__init__.py 2008-01-30 17:05:38 UTC (rev 3876)
+++ trunk/scipy/splinalg/eigen/arpack/__init__.py 2008-02-06 17:56:43 UTC (rev 3897)
@@ -1,3 +0,0 @@
-from info import __doc__
-from arpack import *
-import speigs
Copied: trunk/scipy/splinalg/eigen/arpack/__init__.py (from rev 3896, trunk/scipy/sandbox/arpack/__init__.py)
Deleted: trunk/scipy/splinalg/eigen/arpack/arpack.py
===================================================================
--- trunk/scipy/sandbox/arpack/arpack.py 2008-01-30 17:05:38 UTC (rev 3876)
+++ trunk/scipy/splinalg/eigen/arpack/arpack.py 2008-02-06 17:56:43 UTC (rev 3897)
@@ -1,379 +0,0 @@
-"""
-arpack - Scipy module to find a few eigenvectors and eigenvalues of a matrix
-
-Uses ARPACK: http://www.caam.rice.edu/software/ARPACK/
-
-"""
-__all___=['eigen','eigen_symmetric']
-
-import _arpack
-import numpy as sb
-import warnings
-
-# inspired by iterative.py
-# so inspired, in fact, that some of it was copied directly
-try:
- False, True
-except NameError:
- False, True = 0, 1
-
-_type_conv = {'f':'s', 'd':'d', 'F':'c', 'D':'z'}
-
-class get_matvec:
- methname = 'matvec'
- def __init__(self, obj, *args):
- self.obj = obj
- self.args = args
- if isinstance(obj, sb.matrix):
- self.callfunc = self.type1m
- return
- if isinstance(obj, sb.ndarray):
- self.callfunc = self.type1
- return
- meth = getattr(obj,self.methname,None)
- if not callable(meth):
- raise ValueError, "Object must be an array "\
- "or have a callable %s attribute." % (self.methname,)
-
- self.obj = meth
- self.callfunc = self.type2
-
- def __call__(self, x):
- return self.callfunc(x)
-
- def type1(self, x):
- return sb.dot(self.obj, x)
-
- def type1m(self, x):
- return sb.dot(self.obj.A, x)
-
- def type2(self, x):
- return self.obj(x,*self.args)
-
-
-def eigen(A,k=6,M=None,ncv=None,which='LM',
- maxiter=None,tol=0, return_eigenvectors=True):
- """ Return k eigenvalues and eigenvectors of the matrix A.
-
- Solves A * x[i] = w[i] * x[i], the standard eigenvalue problem for
- w[i] eigenvalues with corresponding eigenvectors x[i].
-
- Inputs:
-
- A -- A matrix, array or an object with matvec(x) method to perform
- the matrix vector product A * x. The sparse matrix formats
- in scipy.sparse are appropriate for A.
-
- k -- The number of eigenvalue/eigenvectors desired
-
- M -- (Not implemented)
- A symmetric positive-definite matrix for the generalized
- eigenvalue problem A * x = w * M * x
-
- Outputs:
-
- w -- An array of k eigenvalues
-
- v -- An array of k eigenvectors, k[i] is the eigenvector corresponding
- to the eigenvector w[i]
-
- Optional Inputs:
-
- ncv -- Number of Lanczos vectors generated, ncv must be greater than k
- and is recommended to be ncv > 2*k
-
- which -- String specifying which eigenvectors to compute.
- Compute the k eigenvalues of:
- 'LM' - largest magnitude.
- 'SM' - smallest magnitude.
- 'LR' - largest real part.
- 'SR' - smallest real part.
- 'LI' - largest imaginary part.
- 'SI' - smallest imaginary part.
-
- maxiter -- Maximum number of Arnoldi update iterations allowed
-
- tol -- Relative accuracy for eigenvalues (stopping criterion)
-
- return_eigenvectors -- True|False, return eigenvectors
-
- """
- try:
- n,ny=A.shape
- n==ny
- except:
- raise AttributeError("matrix is not square")
- if M is not None:
- raise NotImplementedError("generalized eigenproblem not supported yet")
-
- # some defaults
- if ncv is None:
- ncv=2*k+1
- ncv=min(ncv,n)
- if maxiter==None:
- maxiter=n*10
-
- # guess type
- resid = sb.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:
- raise ValueError("k must be positive, k=%d"%k)
- if k == n:
- raise ValueError("k must be less than rank(A), k=%d"%k)
- if maxiter <= 0:
- raise ValueError("maxiter must be positive, maxiter=%d"%maxiter)
- whiches=['LM','SM','LR','SR','LI','SI']
- if which not in whiches:
- raise ValueError("which must be one of %s"%' '.join(whiches))
- if ncv > n or ncv < k:
- raise ValueError("ncv must be k<=ncv<=n, ncv=%s"%ncv)
-
- # assign solver and postprocessor
- ltr = _type_conv[typ]
- eigsolver = _arpack.__dict__[ltr+'naupd']
- eigextract = _arpack.__dict__[ltr+'neupd']
- matvec = get_matvec(A)
-
- v = sb.zeros((n,ncv),typ) # holds Ritz vectors
- resid = sb.zeros(n,typ) # residual
- workd = sb.zeros(3*n,typ) # workspace
- workl = sb.zeros(3*ncv*ncv+6*ncv,typ) # workspace
- iparam = sb.zeros(11,'int') # problem parameters
- ipntr = sb.zeros(14,'int') # pointers into workspaces
- info = 0
- ido = 0
-
- if typ in 'FD':
- rwork = sb.zeros(ncv,typ.lower())
-
- # only supported mode is 1: Ax=lx
- ishfts = 1
- mode1 = 1
- bmat = 'I'
- iparam[0] = ishfts
- iparam[2] = maxiter
- iparam[6] = mode1
-
- while True:
- if typ in 'fd':
- ido,resid,v,iparam,ipntr,info =\
- eigsolver(ido,bmat,which,k,tol,resid,v,iparam,ipntr,
- workd,workl,info)
- else:
- ido,resid,v,iparam,ipntr,info =\
- 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
- break
-
- if info < -1 :
- raise RuntimeError("Error info=%d in arpack"%info)
- return None
- if info == -1:
- warnings.warn("Maximum number of iterations taken: %s"%iparam[2])
-# if iparam[3] != k:
-# warnings.warn("Only %s eigenvalues converged"%iparam[3])
-
-
- # now extract eigenvalues and (optionally) eigenvectors
- rvec = return_eigenvectors
- ierr = 0
- howmny = 'A' # return all eigenvectors
- sselect = sb.zeros(ncv,'int') # unused
- sigmai = 0.0 # no shifts, not implemented
- sigmar = 0.0 # no shifts, not implemented
- workev = sb.zeros(3*ncv,typ)
-
- if typ in 'fd':
- dr=sb.zeros(k+1,typ)
- di=sb.zeros(k+1,typ)
- zr=sb.zeros((n,k+1),typ)
- dr,di,z,info=\
- eigextract(rvec,howmny,sselect,sigmar,sigmai,workev,
- bmat,which,k,tol,resid,v,iparam,ipntr,
- workd,workl,info)
-
- # make eigenvalues complex
- d=dr+1.0j*di
- # futz with the eigenvectors:
- # complex are stored as real,imaginary in consecutive columns
- z=zr.astype(typ.upper())
- for i in range(k): # fix c.c. pairs
- if di[i] > 0 :
- z[:,i]=zr[:,i]+1.0j*zr[:,i+1]
- z[:,i+1]=z[:,i].conjugate()
-
- else:
- d,z,info =\
- eigextract(rvec,howmny,sselect,sigmar,workev,
- bmat,which,k,tol,resid,v,iparam,ipntr,
- workd,workl,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,ncv=None,which='LM',
- maxiter=None,tol=0, return_eigenvectors=True):
- """ Return k eigenvalues and eigenvectors of the real symmetric matrix A.
-
- Solves A * x[i] = w[i] * x[i], the standard eigenvalue problem for
- w[i] eigenvalues with corresponding eigenvectors x[i].
- A must be real and symmetric.
- See eigen() for nonsymmetric or complex symmetric (Hermetian) matrices.
-
- Inputs:
-
- A -- A symmetric matrix, array or an object with matvec(x) method
- to perform the matrix vector product A * x.
- The sparse matrix formats in scipy.sparse are appropriate for A.
-
- k -- The number of eigenvalue/eigenvectors desired
-
- M -- (Not implemented)
- A symmetric positive-definite matrix for the generalized
- eigenvalue problem A * x = w * M * x
-
- Outputs:
-
- w -- An real array of k eigenvalues
-
- v -- An array of k real eigenvectors, k[i] is the eigenvector corresponding
- to the eigenvector w[i]
-
- Optional Inputs:
-
- ncv -- Number of Lanczos vectors generated, ncv must be greater than k
- and is recommended to be ncv > 2*k
-
- which -- String specifying which eigenvectors to compute.
- Compute the k
- 'LA' - largest (algebraic) eigenvalues.
- 'SA' - smallest (algebraic) eigenvalues.
- 'LM' - largest (in magnitude) eigenvalues.
- 'SM' - smallest (in magnitude) eigenvalues.
- 'BE' - eigenvalues, half from each end of the
- spectrum. When NEV is odd, compute one more from the
- high end than from the low end.
-
- maxiter -- Maximum number of Arnoldi update iterations allowed
-
- tol -- Relative accuracy for eigenvalues (stopping criterion)
-
- return_eigenvectors -- True|False, return eigenvectors
-
- """
- try:
- n,ny=A.shape
- n==ny
- except:
- raise AttributeError("matrix is not square")
- if M is not None:
- raise NotImplementedError("generalized eigenproblem not supported yet")
- if ncv is None:
- ncv=2*k+1
- ncv=min(ncv,n)
- if maxiter==None:
- maxiter=n*10
-
-
- # guess type
- resid = sb.zeros(n,'f')
- try:
- typ = A.dtype.char
- except AttributeError:
- typ = A.matvec(resid).dtype.char
- if typ not in 'fd':
- raise ValueError("matrix type must be 'f' or 'd'")
-
- # some sanity checks
- if k <= 0:
- raise ValueError("k must be positive, k=%d"%k)
- if k == n:
- raise ValueError("k must be less than rank(A), k=%d"%k)
- if maxiter <= 0:
- raise ValueError("maxiter must be positive, maxiter=%d"%maxiter)
- whiches=['LM','SM','LA','SA','BE']
- if which not in whiches:
- raise ValueError("which must be one of %s"%' '.join(whiches))
- if ncv > n or ncv < k:
- raise ValueError("ncv must be k<=ncv<=n, ncv=%s"%ncv)
-
- # assign solver and postprocessor
- ltr = _type_conv[typ]
- eigsolver = _arpack.__dict__[ltr+'saupd']
- eigextract = _arpack.__dict__[ltr+'seupd']
- matvec = get_matvec(A)
-
- v = sb.zeros((n,ncv),typ)
- resid = sb.zeros(n,typ)
- workd = sb.zeros(3*n,typ)
- workl = sb.zeros(ncv*(ncv+8),typ)
- iparam = sb.zeros(11,'int')
- ipntr = sb.zeros(11,'int')
- info = 0
- ido = 0
-
- # only supported mode is 1: Ax=lx
- ishfts = 1
- mode1 = 1
- bmat='I'
- iparam[0] = ishfts
- iparam[2] = maxiter
- iparam[6] = mode1
-
-
- 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])
- else:
- break
-
- if info < -1 :
- raise RuntimeError("Error info=%d in arpack"%info)
- return None
- if info == -1:
- warnings.warn("Maximum number of iterations taken: %s"%iparam[2])
-
- # now extract eigenvalues and (optionally) eigenvectors
- rvec = return_eigenvectors
- ierr = 0
- howmny = 'A' # return all eigenvectors
- sselect = sb.zeros(ncv,'int') # unused
- sigma = 0.0 # no shifts, not implemented
-
- d,z,info =\
- eigextract(rvec,howmny,sselect,sigma,
- bmat,which, k,tol,resid,v,iparam[0:7],ipntr,
- workd[0:2*n],workl,ierr)
-
- if ierr != 0:
- raise RuntimeError("Error info=%d in arpack"%info)
- return None
- if return_eigenvectors:
- return d,z
- return d
Copied: trunk/scipy/splinalg/eigen/arpack/arpack.py (from rev 3896, trunk/scipy/sandbox/arpack/arpack.py)
===================================================================
--- trunk/scipy/sandbox/arpack/arpack.py 2008-02-06 01:18:26 UTC (rev 3896)
+++ trunk/scipy/splinalg/eigen/arpack/arpack.py 2008-02-06 17:56:43 UTC (rev 3897)
@@ -0,0 +1,474 @@
+"""
+Find a few eigenvectors and eigenvalues of a matrix.
+
+
+Uses ARPACK: http://www.caam.rice.edu/software/ARPACK/
+
+"""
+# Wrapper implementation notes
+#
+# ARPACK Entry Points
+# -------------------
+# The entry points to ARPACK are
+# - (s,d)seupd : single and double precision symmetric matrix
+# - (s,d,c,z)neupd: single,double,complex,double complex general matrix
+# This wrapper puts the *neupd (general matrix) interfaces in eigen()
+# and the *seupd (symmetric matrix) in eigen_symmetric().
+# There is no Hermetian complex/double complex interface.
+# To find eigenvalues of a Hermetian matrix you
+# must use eigen() and not eigen_symmetric()
+# It might be desirable to handle the Hermetian case differently
+# and, for example, return real eigenvalues.
+
+# Number of eigenvalues returned and complex eigenvalues
+# ------------------------------------------------------
+# The ARPACK nonsymmetric real and double interface (s,d)naupd return
+# eigenvalues and eigenvectors in real (float,double) arrays.
+# Since the eigenvalues and eigenvectors are, in general, complex
+# ARPACK puts the real and imaginary parts in consecutive entries
+# in real-valued arrays. This wrapper puts the real entries
+# into complex data types and attempts to return the requested eigenvalues
+# and eigenvectors.
+
+
+# Solver modes
+# ------------
+# ARPACK and handle shifted and shift-inverse computations
+# for eigenvalues by providing a shift (sigma) and a solver.
+# This is currently not implemented
+
+__docformat__ = "restructuredtext en"
+
+__all___=['eigen','eigen_symmetric']
+
+import warnings
+
+import _arpack
+import numpy as np
+from scipy.splinalg.interface import aslinearoperator
+
+_type_conv = {'f':'s', 'd':'d', 'F':'c', 'D':'z'}
+_ndigits = {'f':5, 'd':12, 'F':5, 'D':12}
+
+
+def eigen(A, k=6, M=None, sigma=None, which='LM',
+ ncv=None, maxiter=None, tol=0,
+ return_eigenvectors=True):
+ """Find k eigenvalues and eigenvectors of the square matrix A.
+
+ Solves A * x[i] = w[i] * x[i], the standard eigenvalue problem for
+ w[i] eigenvalues with corresponding eigenvectors x[i].
+
+
+ Parameters
+ ----------
+ A : A : matrix, array, or object with matvec(x) method
+ An N x N matrix, array, or an object with matvec(x) method to perform
+ the matrix vector product A * x. The sparse matrix formats
+ in scipy.sparse are appropriate for A.
+
+ k : integer
+ The number of eigenvalues and eigenvectors desired
+
+ Returns
+ -------
+ w : array
+ Array of k eigenvalues
+
+ v : array
+ An array of k eigenvectors
+ The v[i] is the eigenvector corresponding to the eigenvector w[i]
+
+ Other Parameters
+ ----------------
+
+ M : matrix or array
+ (Not implemented)
+ A symmetric positive-definite matrix for the generalized
+ eigenvalue problem A * x = w * M * x
+
+ sigma : real or complex
+ (Not implemented)
+ Find eigenvalues near sigma. Shift spectrum by sigma.
+
+
+ ncv : integer
+ The number of Lanczos vectors generated
+ ncv must be greater than k; it is recommended that ncv > 2*k
+
+ which : string
+ Which k eigenvectors and eigenvalues to find:
+ - 'LM' : largest magnitude
+ - 'SM' : smallest magnitude
+ - 'LR' : largest real part
+ - 'SR' : smallest real part
+ - 'LI' : largest imaginary part
+ - 'SI' : smallest imaginary part
+
+ maxiter : integer
+ Maximum number of Arnoldi update iterations allowed
+
+ tol : float
+ Relative accuracy for eigenvalues (stopping criterion)
+
+ return_eigenvectors : boolean
+ Return eigenvectors (True) in addition to eigenvalues
+
+ See Also
+ --------
+ eigen_symmetric : eigenvalues and eigenvectors for symmetric matrix A
+
+ Notes
+ -----
+
+ Examples
+ --------
+
+ """
+ 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]
+
+ if M is not None:
+ raise NotImplementedError("generalized eigenproblem not supported yet")
+
+ if sigma is not None:
+ raise NotImplementedError("shifted eigenproblem not supported yet")
+
+
+ # some defaults
+ if ncv is None:
+ ncv=2*k+1
+ ncv=min(ncv,n)
+ if maxiter==None:
+ maxiter=n*10
+
+ # 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:
+ raise ValueError("k must be positive, k=%d"%k)
+ if k == n:
+ raise ValueError("k must be less than rank(A), k=%d"%k)
+ if maxiter <= 0:
+ raise ValueError("maxiter must be positive, maxiter=%d"%maxiter)
+ whiches=['LM','SM','LR','SR','LI','SI']
+ if which not in whiches:
+ raise ValueError("which must be one of %s"%' '.join(whiches))
+ if ncv > n or ncv < k:
+ raise ValueError("ncv must be k<=ncv<=n, ncv=%s"%ncv)
+
+ # assign solver and postprocessor
+ 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())
+
+ # only supported mode is 1: Ax=lx
+ ishfts = 1
+ mode1 = 1
+ bmat = 'I'
+ iparam[0] = ishfts
+ iparam[2] = maxiter
+ iparam[6] = mode1
+
+ while True:
+ if typ in 'fd':
+ ido,resid,v,iparam,ipntr,info =\
+ eigsolver(ido,bmat,which,k,tol,resid,v,iparam,ipntr,
+ workd,workl,info)
+ else:
+ ido,resid,v,iparam,ipntr,info =\
+ 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
+ break
+
+ if info < -1 :
+ raise RuntimeError("Error info=%d in arpack"%info)
+ return None
+ if info == -1:
+ warnings.warn("Maximum number of iterations taken: %s"%iparam[2])
+# if iparam[3] != k:
+# warnings.warn("Only %s eigenvalues converged"%iparam[3])
+
+
+ # now extract eigenvalues and (optionally) eigenvectors
+ rvec = return_eigenvectors
+ ierr = 0
+ howmny = 'A' # return all eigenvectors
+ sselect = np.zeros(ncv,'int') # unused
+ sigmai = 0.0 # no shifts, not implemented
+ sigmar = 0.0 # no shifts, not implemented
+ workev = np.zeros(3*ncv,typ)
+
+ if typ in 'fd':
+ dr=np.zeros(k+1,typ)
+ di=np.zeros(k+1,typ)
+ zr=np.zeros((n,k+1),typ)
+ dr,di,zr,info=\
+ eigextract(rvec,howmny,sselect,sigmar,sigmai,workev,
+ bmat,which,k,tol,resid,v,iparam,ipntr,
+ workd,workl,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(typ.upper())
+ eps=np.finfo(typ).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 eigenvalues and eigenvectors
+ # Return the ones specified by the keyword "which"
+
+ # cut at approx precision for sorting
+ rd=np.round(d,decimals=_ndigits[typ])
+ if 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 which in ['LR','LM','LI']:
+ d=d[ind[-k:]]
+ z=z[:,ind[-k:]]
+ if which in ['SR','SM','SI']:
+ d=d[ind[:k]]
+ z=z[:,ind[:k]]
+
+
+ else:
+ # complex is so much simpler...
+ d,z,info =\
+ eigextract(rvec,howmny,sselect,sigmar,workev,
+ bmat,which,k,tol,resid,v,iparam,ipntr,
+ workd,workl,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',
+ ncv=None, maxiter=None, tol=0, v0=None,
+ return_eigenvectors=True):
+ """Find k eigenvalues and eigenvectors of the real symmetric
+ square matrix A.
+
+ Solves A * x[i] = w[i] * x[i], the standard eigenvalue problem for
+ w[i] eigenvalues with corresponding eigenvectors x[i].
+
+
+ Parameters
+ ----------
+ A : matrix or array with real entries or object with matvec(x) method
+ An N x N real symmetric matrix or array or an object with matvec(x)
+ method to perform the matrix vector product A * x. The sparse
+ matrix formats in scipy.sparse are appropriate for A.
+
+ k : integer
+ The number of eigenvalues and eigenvectors desired
+
+ Returns
+ -------
+ w : array
+ Array of k eigenvalues
+
+ v : array
+ An array of k eigenvectors
+ The v[i] is the eigenvector corresponding to the eigenvector w[i]
+
+ Other Parameters
+ ----------------
+ M : matrix or array
+ (Not implemented)
+ A symmetric positive-definite matrix for the generalized
+ eigenvalue problem A * x = w * M * x
+
+
+ sigma : real
+ (Not implemented)
+ Find eigenvalues near sigma. Shift spectrum by sigma.
+
+
+ ncv : integer
+ The number of Lanczos vectors generated
+ ncv must be greater than k; it is recommended that ncv > 2*k
+
+ which : string
+ Which k eigenvectors and eigenvalues to find:
+ - 'LA' : Largest (algebraic) eigenvalues
+ - 'SA' : Smallest (algebraic) eigenvalues
+ - 'LM' : Largest (in magnitude) eigenvalues
+ - 'SM' : Smallest (in magnitude) eigenvalues
+ - 'BE' : Half (k/2) from each end of the spectrum
+ When k is odd, return one more (k/2+1) from the high end
+
+ maxiter : integer
+ Maximum number of Arnoldi update iterations allowed
+
+ tol : float
+ Relative accuracy for eigenvalues (stopping criterion)
+
+ return_eigenvectors : boolean
+ Return eigenvectors (True) in addition to eigenvalues
+
+ See Also
+ --------
+ eigen : eigenvalues and eigenvectors for a general (nonsymmetric) matrix A
+
+ Notes
+ -----
+
+ Examples
+ --------
+ """
+ A = aslinearoperator(A)
+
+ if A.shape[0] != A.shape[1]:
+ raise ValueError('expected square matrix (shape=%s)' % shape)
+
+ n = A.shape[0]
+
+ 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
+
+
+ # 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)
+ if k == n:
+ raise ValueError("k must be less than rank(A), k=%d"%k)
+ if maxiter <= 0:
+ raise ValueError("maxiter must be positive, maxiter=%d"%maxiter)
+ whiches=['LM','SM','LA','SA','BE']
+ if which not in whiches:
+ raise ValueError("which must be one of %s"%' '.join(whiches))
+ if ncv > n or ncv < k:
+ raise ValueError("ncv must be k<=ncv<=n, ncv=%s"%ncv)
+
+ # assign solver and postprocessor
+ ltr = _type_conv[typ]
+ eigsolver = _arpack.__dict__[ltr+'saupd']
+ eigextract = _arpack.__dict__[ltr+'seupd']
+ matvec = A.matvec
+
+ 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
+
+ # only supported mode is 1: Ax=lx
+ ishfts = 1
+ mode1 = 1
+ bmat='I'
+ iparam[0] = ishfts
+ iparam[2] = maxiter
+ iparam[6] = mode1
+
+ 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])
+ else:
+ break
+
+ if info < -1 :
+ raise RuntimeError("Error info=%d in arpack"%info)
+ return None
+ if info == -1:
+ warnings.warn("Maximum number of iterations taken: %s"%iparam[2])
+
+ # now extract eigenvalues and (optionally) eigenvectors
+ rvec = return_eigenvectors
+ ierr = 0
+ howmny = 'A' # return all eigenvectors
+ sselect = np.zeros(ncv,'int') # unused
+ sigma = 0.0 # no shifts, not implemented
+
+ d,z,info =\
+ eigextract(rvec,howmny,sselect,sigma,
+ bmat,which, k,tol,resid,v,iparam[0:7],ipntr,
+ workd[0:2*n],workl,ierr)
+
+ if ierr != 0:
+ raise RuntimeError("Error info=%d in arpack"%info)
+ return None
+ if return_eigenvectors:
+ return d,z
+ return d
Deleted: trunk/scipy/splinalg/eigen/arpack/arpack.pyf.src
===================================================================
--- trunk/scipy/sandbox/arpack/arpack.pyf.src 2008-01-30 17:05:38 UTC (rev 3876)
+++ trunk/scipy/splinalg/eigen/arpack/arpack.pyf.src 2008-02-06 17:56:43 UTC (rev 3897)
@@ -1,207 +0,0 @@
-! -*- f90 -*-
-! Note: the context of this file is case sensitive.
-
-python module _arpack ! in
- <_rd=real,double precision>
- <_cd=complex,double complex>
- interface ! in :_arpack
- subroutine <s,d>saupd(ido,bmat,n,which,nev,tol,resid,ncv,v,ldv,iparam,ipntr,workd,workl,lworkl,info) ! in :_arpack:src/ssaupd.f
- integer intent(in,out):: ido
- character*1 :: bmat
- integer optional,check(len(resid)>=n),depend(resid) :: n=len(resid)
- character*2 :: which
- integer :: nev
- <_rd> :: tol
- <_rd> dimension(n),intent(in,out) :: resid
- integer optional,check(shape(v,1)==ncv),depend(v) :: ncv=shape(v,1)
- <_rd> dimension(ldv,ncv),intent(in,out) :: v
- integer optional,check(shape(v,0)==ldv),depend(v) :: ldv=shape(v,0)
- integer dimension(11),intent(in,out) :: iparam
- integer dimension(11),intent(in,out) :: ipntr
- <_rd> dimension(3 * n),depend(n),intent(inout) :: workd
- <_rd> dimension(lworkl),intent(inout) :: workl
- integer optional,check(len(workl)>=lworkl),depend(workl) :: lworkl=len(workl)
- integer intent(in,out):: info
- end subroutine <s,d>saupd
-
- subroutine <s,d>seupd(rvec,howmny,select,d,z,ldz,sigma,bmat,n,which,nev,tol,resid,ncv,v,ldv,iparam,ipntr,workd,workl,lworkl,info) ! in :_arpack:src/sseupd.f
- logical :: rvec
- character :: howmny
- logical dimension(ncv) :: select
- <_rd> dimension(nev),intent(out),depend(nev) :: d
- <_rd> dimension(n,nev),intent(out),depend(nev) :: z
- integer optional,check(shape(z,0)==ldz),depend(z) :: ldz=shape(z,0)
- <_rd> :: sigma
- character :: bmat
- integer optional,check(len(resid)>=n),depend(resid) :: n=len(resid)
- character*2 :: which
- integer :: nev
- <_rd> :: tol
- <_rd> dimension(n) :: resid
- integer optional,check(len(select)>=ncv),depend(select) :: ncv=len(select)
- <_rd> dimension(ldv,ncv),depend(ncv) :: v
- integer optional,check(shape(v,0)==ldv),depend(v) :: ldv=shape(v,0)
- integer dimension(7) :: iparam
- integer dimension(11) :: ipntr
- <_rd> dimension(2 * n),depend(n) :: workd
- <_rd> dimension(lworkl) :: workl
- integer optional,check(len(workl)>=lworkl),depend(workl) :: lworkl=len(workl)
- integer intent(in,out):: info
- end subroutine <s,d>seupd
-
- subroutine <s,d>naupd(ido,bmat,n,which,nev,tol,resid,ncv,v,ldv,iparam,ipntr,workd,workl,lworkl,info) ! in :_arpack:src/snaupd.f
- integer intent(in,out):: ido
- character*1 :: bmat
- integer optional,check(len(resid)>=n),depend(resid) :: n=len(resid)
- character*2 :: which
- integer :: nev
- <_rd> :: tol
- <_rd> dimension(n),intent(in,out) :: resid
- integer optional,check(shape(v,1)==ncv),depend(v) :: ncv=shape(v,1)
- <_rd> dimension(ldv,ncv),intent(in,out) :: v
- integer optional,check(shape(v,0)==ldv),depend(v) :: ldv=shape(v,0)
- integer dimension(11),intent(in,out) :: iparam
- integer dimension(14),intent(in,out) :: ipntr
- <_rd> dimension(3 * n),depend(n),intent(inout) :: workd
- <_rd> dimension(lworkl),intent(inout) :: workl
- integer optional,check(len(workl)>=lworkl),depend(workl) :: lworkl=len(workl)
- integer intent(in,out):: info
- end subroutine <s,d>naupd
-
- subroutine <s,d>neupd(rvec,howmny,select,dr,di,z,ldz,sigmar,sigmai,workev,bmat,n,which,nev,tol,resid,ncv,v,ldv,iparam,ipntr,workd,workl,lworkl,info) ! in ARPACK/SRC/sneupd.f
- logical :: rvec
- character :: howmny
- logical dimension(ncv) :: select
- <_rd> dimension(nev + 1),depend(nev),intent(out) :: dr
- <_rd> dimension(nev + 1),depend(nev),intent(out) :: di
- <_rd> dimension(n,nev+1),depend(n,nev),intent(out) :: z
- integer optional,check(shape(z,0)==ldz),depend(z) :: ldz=shape(z,0)
- <_rd> :: sigmar
- <_rd> :: sigmai
- <_rd> dimension(3 * ncv),depend(ncv) :: workev
- character :: bmat
- integer optional,check(len(resid)>=n),depend(resid) :: n=len(resid)
- character*2 :: which
- integer :: nev
- <_rd> :: tol
- <_rd> dimension(n) :: resid
- integer optional,check(len(select)>=ncv),depend(select) :: ncv=len(select)
- <_rd> dimension(n,ncv),depend(n,ncv) :: v
- integer optional,check(shape(v,0)==ldv),depend(v) :: ldv=shape(v,0)
- integer dimension(11) :: iparam
- integer dimension(14) :: ipntr
- <_rd> dimension(3 * n),depend(n):: workd
- <_rd> dimension(lworkl) :: workl
- integer optional,check(len(workl)>=lworkl),depend(workl) :: lworkl=len(workl)
- integer intent(in,out):: info
- end subroutine <s,d>neupd
-
- subroutine <c,z>naupd(ido,bmat,n,which,nev,tol,resid,ncv,v,ldv,iparam,ipntr,workd,workl,lworkl,rwork,info) ! in :_arpack:src/snaupd.f
- integer intent(in,out):: ido
- character*1 :: bmat
- integer optional,check(len(resid)>=n),depend(resid) :: n=len(resid)
- character*2 :: which
- integer :: nev
- <_rd> :: tol
- <_cd> dimension(n),intent(in,out) :: resid
- integer optional,check(shape(v,1)==ncv),depend(v) :: ncv=shape(v,1)
- <_cd> dimension(ldv,ncv),intent(in,out) :: v
- integer optional,check(shape(v,0)==ldv),depend(v) :: ldv=shape(v,0)
- integer dimension(11),intent(in,out) :: iparam
- integer dimension(14),intent(in,out) :: ipntr
- <_cd> dimension(3 * n),depend(n),intent(inout) :: workd
- <_cd> dimension(lworkl),intent(inout) :: workl
- integer optional,check(len(workl)>=lworkl),depend(workl) :: lworkl=len(workl)
- <_rd> dimension(ncv),depend(ncv),intent(inout) :: rwork
- integer intent(in,out):: info
- end subroutine <c,z>naupd
-
- subroutine <c,z>neupd(rvec,howmny,select,d,z,ldz,sigma,workev,bmat,n,which,nev,tol,resid,ncv,v,ldv,iparam,ipntr,workd,workl,lworkl,rwork,info) ! in :_arpack:src/sneupd.f
- logical :: rvec
- character :: howmny
- logical dimension(ncv) :: select
- <_cd> dimension(nev),depend(nev),intent(out) :: d
- <_cd> dimension(n,nev), depend(nev),intent(out) :: z
- integer optional,check(shape(z,0)==ldz),depend(z) :: ldz=shape(z,0)
- <_cd> :: sigma
- <_cd> dimension(3 * ncv),depend(ncv) :: workev
- character :: bmat
- integer optional,check(len(resid)>=n),depend(resid) :: n=len(resid)
- character*2 :: which
- integer :: nev
- <_rd> :: tol
- <_cd> dimension(n) :: resid
- integer optional,check(len(select)>=ncv),depend(select) :: ncv=len(select)
- <_cd> dimension(ldv,ncv),depend(ncv) :: v
- integer optional,check(shape(v,0)==ldv),depend(v) :: ldv=shape(v,0)
- integer dimension(11) :: iparam
- integer dimension(14) :: ipntr
- <_cd> dimension(3 * n),depend(n) :: workd
- <_cd> dimension(lworkl) :: workl
- integer optional,check(len(workl)>=lworkl),depend(workl) :: lworkl=len(workl)
- <_rd> dimension(ncv),depend(ncv) :: rwork
- integer intent(in,out):: info
- end subroutine <c,z>neupd
- integer :: logfil
- integer :: ndigit
- integer :: mgetv0
- integer :: msaupd
- integer :: msaup2
- integer :: msaitr
- integer :: mseigt
- integer :: msapps
- integer :: msgets
- integer :: mseupd
- integer :: mnaupd
- integer :: mnaup2
- integer :: mnaitr
- integer :: mneigh
- integer :: mnapps
- integer :: mngets
- integer :: mneupd
- integer :: mcaupd
- integer :: mcaup2
- integer :: mcaitr
- integer :: mceigh
- integer :: mcapps
- integer :: mcgets
- integer :: mceupd
- integer :: nopx
- integer :: nbx
- integer :: nrorth
- integer :: nitref
- integer :: nrstrt
- real :: tsaupd
- real :: tsaup2
- real :: tsaitr
- real :: tseigt
- real :: tsgets
- real :: tsapps
- real :: tsconv
- real :: tnaupd
- real :: tnaup2
- real :: tnaitr
- real :: tneigh
- real :: tngets
- real :: tnapps
- real :: tnconv
- real :: tcaupd
- real :: tcaup2
- real :: tcaitr
- real :: tceigh
- real :: tcgets
- real :: tcapps
- real :: tcconv
- real :: tmvopx
- real :: tmvbx
- real :: tgetv0
- real :: titref
- real :: trvec
- common /debug/ logfil,ndigit,mgetv0,msaupd,msaup2,msaitr,mseigt,msapps,msgets,mseupd,mnaupd,mnaup2,mnaitr,mneigh,mnapps,mngets,mneupd,mcaupd,mcaup2,mcaitr,mceigh,mcapps,mcgets,mceupd
- common /timing/ nopx,nbx,nrorth,nitref,nrstrt,tsaupd,tsaup2,tsaitr,tseigt,tsgets,tsapps,tsconv,tnaupd,tnaup2,tnaitr,tneigh,tngets,tnapps,tnconv,tcaupd,tcaup2,tcaitr,tceigh,tcgets,tcapps,tcconv,tmvopx,tmvbx,tgetv0,titref,trvec
-
- end interface
-end python module _arpack
-
-! This file was auto-generated with f2py (version:2_3198).
-! See http://cens.ioc.ee/projects/f2py2e/
Copied: trunk/scipy/splinalg/eigen/arpack/arpack.pyf.src (from rev 3896, trunk/scipy/sandbox/arpack/arpack.pyf.src)
Deleted: trunk/scipy/splinalg/eigen/arpack/info.py
===================================================================
--- trunk/scipy/sandbox/arpack/info.py 2008-01-30 17:05:38 UTC (rev 3876)
+++ trunk/scipy/splinalg/eigen/arpack/info.py 2008-02-06 17:56:43 UTC (rev 3897)
@@ -1,14 +0,0 @@
-"""
-Eigenvalue solver using iterative methods.
-
-Find k eigenvectors and eigenvalues of a matrix A using the
-Arnoldi/Lanczos iterative methods from ARPACK.
-
-These methods are most useful for large sparse matrices.
-
- - eigen(A,k)
- - eigen_symmetric(A,k)
-
-"""
-global_symbols = []
-postpone_import = 1
Copied: trunk/scipy/splinalg/eigen/arpack/info.py (from rev 3896, trunk/scipy/sandbox/arpack/info.py)
===================================================================
--- trunk/scipy/sandbox/arpack/info.py 2008-02-06 01:18:26 UTC (rev 3896)
+++ trunk/scipy/splinalg/eigen/arpack/info.py 2008-02-06 17:56:43 UTC (rev 3897)
@@ -0,0 +1,21 @@
+"""
+Eigenvalue solver using iterative methods.
+
+Find k eigenvectors and eigenvalues of a matrix A using the
+Arnoldi/Lanczos iterative methods from ARPACK.
+
+These methods are most useful for large sparse matrices.
+
+ - eigen(A,k)
+ - eigen_symmetric(A,k)
+
+Reference
+---------
+ - http://www.caam.rice.edu/
+software/ARPACK/
+ - http://www.caam.rice.edu/software/ARPACK/UG/ug.html
+ - http://books.google.com/books?hl=en&id=4E9PY7NT8a0C&dq=arpack+users+guide
+
+"""
+global_symbols = []
+postpone_import = 1
Deleted: trunk/scipy/splinalg/eigen/arpack/setup.py
===================================================================
--- trunk/scipy/sandbox/arpack/setup.py 2008-01-30 17:05:38 UTC (rev 3876)
+++ trunk/scipy/splinalg/eigen/arpack/setup.py 2008-02-06 17:56:43 UTC (rev 3897)
@@ -1,36 +0,0 @@
-#!/usr/bin/env python
-
-from os.path import join
-
-def configuration(parent_package='',top_path=None):
- from numpy.distutils.system_info import get_info, NotFoundError
- from numpy.distutils.misc_util import Configuration
-
- lapack_opt = get_info('lapack_opt')
-
- if not lapack_opt:
- raise NotFoundError,'no lapack/blas resources found'
-
- config = Configuration('arpack', parent_package, top_path)
-
- arpack_sources=[join('ARPACK','SRC', '*.f')]
- arpack_sources.extend([join('ARPACK','UTIL', '*.f')])
-# arpack_sources.extend([join('ARPACK','BLAS', '*.f')])
- arpack_sources.extend([join('ARPACK','LAPACK', '*.f')])
-
- config.add_library('arpack', sources=arpack_sources,
- include_dirs=[join('ARPACK', 'SRC')])
-
-
- config.add_extension('_arpack',
- sources='arpack.pyf.src',
- libraries=['arpack'],
- extra_info = lapack_opt
- )
-
- config.add_data_dir('tests')
- return config
-
-if __name__ == '__main__':
- from numpy.distutils.core import setup
- setup(**configuration(top_path='').todict())
Copied: trunk/scipy/splinalg/eigen/arpack/setup.py (from rev 3896, trunk/scipy/sandbox/arpack/setup.py)
===================================================================
--- trunk/scipy/sandbox/arpack/setup.py 2008-02-06 01:18:26 UTC (rev 3896)
+++ trunk/scipy/splinalg/eigen/arpack/setup.py 2008-02-06 17:56:43 UTC (rev 3897)
@@ -0,0 +1,37 @@
+#!/usr/bin/env python
+
+from os.path import join
+
+def configuration(parent_package='',top_path=None):
+ from numpy.distutils.system_info import get_info, NotFoundError
+ from numpy.distutils.misc_util import Configuration
+
+ config = Configuration('arpack',parent_package,top_path)
+
+ lapack_opt = get_info('lapack_opt')
+
+ if not lapack_opt:
+ raise NotFoundError,'no lapack/blas resources found'
+
+ config = Configuration('arpack', parent_package, top_path)
+
+ arpack_sources=[join('ARPACK','SRC', '*.f')]
+ arpack_sources.extend([join('ARPACK','UTIL', '*.f')])
+ arpack_sources.extend([join('ARPACK','LAPACK', '*.f')])
+
+ config.add_library('arpack', sources=arpack_sources,
+ include_dirs=[join('ARPACK', 'SRC')])
+
+
+ config.add_extension('_arpack',
+ sources='arpack.pyf.src',
+ libraries=['arpack'],
+ extra_info = lapack_opt
+ )
+
+ config.add_data_dir('tests')
+ return config
+
+if __name__ == '__main__':
+ from numpy.distutils.core import setup
+ setup(**configuration(top_path='').todict())
Deleted: trunk/scipy/splinalg/eigen/arpack/speigs.py
===================================================================
--- trunk/scipy/sandbox/arpack/speigs.py 2008-01-30 17:05:38 UTC (rev 3876)
+++ trunk/scipy/splinalg/eigen/arpack/speigs.py 2008-02-06 17:56:43 UTC (rev 3897)
@@ -1,225 +0,0 @@
-import numpy as N
-import _arpack
-import warnings
-
-__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 = N.zeros(14, N.int32) # Pointers into memory structure used by F77 calls
- d = N.zeros((ncv, 3), N.float64, order='FORTRAN') # Temp workspace
- # Temp workspace/error residuals upon iteration completion
- resid = N.zeros(n, N.float64)
- workd = N.zeros(3*n, N.float64) # workspace
- workl = N.zeros(3*ncv*ncv+6*ncv, N.float64) # more workspace
- # Storage for the Arnoldi basis vectors
- v = N.zeros((n, ncv), dtype=N.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 = N.zeros(3*ncv, N.float64) # More workspace
- select = N.zeros(ncv, N.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 N.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 = N.zeros(11, N.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)
Copied: trunk/scipy/splinalg/eigen/arpack/speigs.py (from rev 3896, trunk/scipy/sandbox/arpack/speigs.py)
===================================================================
--- trunk/scipy/sandbox/arpack/speigs.py 2008-02-06 01:18:26 UTC (rev 3896)
+++ trunk/scipy/splinalg/eigen/arpack/speigs.py 2008-02-06 17:56:43 UTC (rev 3897)
@@ -0,0 +1,225 @@
+import numpy as np
+import _arpack
+import warnings
+
+__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)
Copied: trunk/scipy/splinalg/eigen/arpack/tests (from rev 3896, trunk/scipy/sandbox/arpack/tests)
Modified: trunk/scipy/splinalg/eigen/arpack/tests/test_arpack.py
===================================================================
--- trunk/scipy/sandbox/arpack/tests/test_arpack.py 2008-02-06 01:18:26 UTC (rev 3896)
+++ trunk/scipy/splinalg/eigen/arpack/tests/test_arpack.py 2008-02-06 17:56:43 UTC (rev 3897)
@@ -1,33 +1,26 @@
#!/usr/bin/env python
-__usage__ = """
-First ensure that scipy core modules are installed.
-Build interface to arpack
- python setup.py build
-Run tests locally:
- python tests/test_arpack.py [-l<int>] [-v<int>]
-
-"""
-
from scipy.testing import *
-from scipy.sandbox.arpack import *
-import numpy
+from scipy.splinalg.eigen.arpack import eigen_symmetric,eigen
from scipy.linalg import eig,eigh,norm
+from numpy import array,abs,take,concatenate,dot
+
+
class TestEigenNonsymmetric(TestCase):
def get_a1(self,typ):
- mat=numpy.array([[-2., -8., 1., 2., -5.],
- [ 6., 6., 0., 2., 1.],
- [ 0., 4., -2., 11., 0.],
- [ 1., 6., 1., 0., -4.],
- [ 2., -6., 4., 9., -3]],typ)
+ mat=array([[-2., -8., 1., 2., -5.],
+ [ 6., 6., 0., 2., 1.],
+ [ 0., 4., -2., 11., 0.],
+ [ 1., 6., 1., 0., -4.],
+ [ 2., -6., 4., 9., -3]],typ)
- w=numpy.array([-2.21691+8.59661*1j,-2.21691-8.59661*1j,\
- 4.45961+3.80078*1j, 4.45961-3.80078*1j,\
- -5.48541+0j],typ.upper())
+ w=array([-2.21691+8.59661*1j,-2.21691-8.59661*1j,\
+ 4.45961+3.80078*1j, 4.45961-3.80078*1j,\
+ -5.48541+0j],typ.upper())
return mat,w
@@ -35,9 +28,9 @@
a,aw = self.get_a1(typ)
w,v = eigen(a,k,which='LM')
for i in range(k):
- assert_array_almost_equal(sb.dot(a,v[:,i]),w[i]*v[:,i],decimal=5)
- exact=numpy.abs(aw)
- num=numpy.abs(w)
+ assert_array_almost_equal(dot(a,v[:,i]),w[i]*v[:,i],decimal=5)
+ exact=abs(aw)
+ num=abs(w)
exact.sort()
num.sort()
assert_array_almost_equal(num[-k:],exact[-k:],decimal=5)
@@ -46,9 +39,9 @@
a,aw = self.get_a1(typ)
w,v = eigen(a,k,which='SM')
for i in range(k):
- assert_array_almost_equal(sb.dot(a,v[:,i]),w[i]*v[:,i],decimal=5)
- exact=numpy.abs(aw)
- num=numpy.abs(w)
+ assert_array_almost_equal(dot(a,v[:,i]),w[i]*v[:,i],decimal=5)
+ exact=abs(aw)
+ num=abs(w)
exact.sort()
num.sort()
assert_array_almost_equal(num[:k],exact[:k],decimal=5)
@@ -58,9 +51,9 @@
a,aw = self.get_a1(typ)
w,v = eigen(a,k,which='LR')
for i in range(k):
- assert_array_almost_equal(sb.dot(a,v[:,i]),w[i]*v[:,i],decimal=5)
- exact=numpy.real(aw)
- num=numpy.real(w)
+ assert_array_almost_equal(dot(a,v[:,i]),w[i]*v[:,i],decimal=5)
+ exact=aw.real
+ num=w.real
exact.sort()
num.sort()
assert_array_almost_equal(num[-k:],exact[-k:],decimal=5)
@@ -69,9 +62,9 @@
a,aw = self.get_a1(typ)
w,v = eigen(a,k,which='SR')
for i in range(k):
- assert_array_almost_equal(sb.dot(a,v[:,i]),w[i]*v[:,i],decimal=5)
- exact=numpy.real(aw)
- num=numpy.real(w)
+ assert_array_almost_equal(dot(a,v[:,i]),w[i]*v[:,i],decimal=5)
+ exact=aw.real
+ num=w.real
exact.sort()
num.sort()
assert_array_almost_equal(num[:k],exact[:k],decimal=5)
@@ -81,11 +74,11 @@
a,aw = self.get_a1(typ)
w,v = eigen(a,k,which='LI')
for i in range(k):
- assert_array_almost_equal(sb.dot(a,v[:,i]),w[i]*v[:,i],decimal=5)
+ assert_array_almost_equal(dot(a,v[:,i]),w[i]*v[:,i],decimal=5)
print w
print aw
- exact=numpy.imag(aw)
- num=numpy.imag(w)
+ exact=aw.imag
+ num=w.imag
exact.sort()
num.sort()
assert_array_almost_equal(num[-k:],exact[-k:],decimal=5)
@@ -94,9 +87,9 @@
a,aw = self.get_a1(typ)
w,v = eigen(a,k,which='SI')
for i in range(k):
- assert_array_almost_equal(sb.dot(a,v[:,i]),w[i]*v[:,i],decimal=5)
- exact=numpy.imag(aw)
- num=numpy.imag(w)
+ assert_array_almost_equal(dot(a,v[:,i]),w[i]*v[:,i],decimal=5)
+ exact=aw.imag
+ num=w.imag
exact.sort()
num.sort()
print num
@@ -121,15 +114,15 @@
class TestEigenComplexNonsymmetric(TestCase):
def get_a1(self,typ):
- mat=numpy.array([[-2., -8., 1., 2., -5.],
- [ 6., 6., 0., 2., 1.],
- [ 0., 4., -2., 11., 0.],
- [ 1., 6., 1., 0., -4.],
- [ 2., -6., 4., 9., -3]],typ)
+ mat=array([[-2., -8., 1., 2., -5.],
+ [ 6., 6., 0., 2., 1.],
+ [ 0., 4., -2., 11., 0.],
+ [ 1., 6., 1., 0., -4.],
+ [ 2., -6., 4., 9., -3]],typ)
- w=numpy.array([-2.21691+8.59661*1j,-2.21691-8.59661*1j,\
- 4.45961+3.80078*1j, 4.45961-3.80078*1j,\
- -5.48541+0j],typ.upper())
+ w=array([-2.21691+8.59661*1j,-2.21691-8.59661*1j,\
+ 4.45961+3.80078*1j, 4.45961-3.80078*1j,\
+ -5.48541+0j],typ.upper())
return mat,w
@@ -137,9 +130,9 @@
a,aw = self.get_a1(typ)
w,v = eigen(a,k,which='LM')
for i in range(k):
- assert_array_almost_equal(sb.dot(a,v[:,i]),w[i]*v[:,i],decimal=5)
- exact=numpy.abs(aw)
- num=numpy.abs(w)
+ assert_array_almost_equal(dot(a,v[:,i]),w[i]*v[:,i],decimal=5)
+ exact=abs(aw)
+ num=abs(w)
exact.sort()
num.sort()
assert_array_almost_equal(num,exact[-k:],decimal=5)
@@ -148,9 +141,9 @@
a,aw = self.get_a1(typ)
w,v = eigen(a,k,which='SM')
for i in range(k):
- assert_array_almost_equal(sb.dot(a,v[:,i]),w[i]*v[:,i],decimal=5)
- exact=numpy.abs(aw)
- num=numpy.abs(w)
+ assert_array_almost_equal(dot(a,v[:,i]),w[i]*v[:,i],decimal=5)
+ exact=abs(aw)
+ num=abs(w)
exact.sort()
num.sort()
assert_array_almost_equal(num,exact[:k],decimal=5)
@@ -160,9 +153,9 @@
a,aw = self.get_a1(typ)
w,v = eigen(a,k,which='LR')
for i in range(k):
- assert_array_almost_equal(sb.dot(a,v[:,i]),w[i]*v[:,i],decimal=5)
- exact=numpy.real(aw)
- num=numpy.real(w)
+ assert_array_almost_equal(dot(a,v[:,i]),w[i]*v[:,i],decimal=5)
+ exact=aw.real
+ num=w.real
exact.sort()
num.sort()
assert_array_almost_equal(num,exact[-k:],decimal=5)
@@ -171,9 +164,9 @@
a,aw = self.get_a1(typ)
w,v = eigen(a,k,which='SR')
for i in range(k):
- assert_array_almost_equal(sb.dot(a,v[:,i]),w[i]*v[:,i],decimal=5)
- exact=numpy.real(aw)
- num=numpy.real(w)
+ assert_array_almost_equal(dot(a,v[:,i]),w[i]*v[:,i],decimal=5)
+ exact=aw.real
+ num=w.real
exact.sort()
num.sort()
assert_array_almost_equal(num,exact[:k],decimal=5)
@@ -183,9 +176,9 @@
a,aw = self.get_a1(typ)
w,v = eigen(a,k,which='LI')
for i in range(k):
- assert_array_almost_equal(sb.dot(a,v[:,i]),w[i]*v[:,i],decimal=5)
- exact=numpy.imag(aw)
- num=numpy.imag(w)
+ assert_array_almost_equal(dot(a,v[:,i]),w[i]*v[:,i],decimal=5)
+ exact=aw.imag
+ num=w.imag
exact.sort()
num.sort()
assert_array_almost_equal(num,exact[-k:],decimal=5)
@@ -194,23 +187,23 @@
a,aw = self.get_a1(typ)
w,v = eigen(a,k,which='SI')
for i in range(k):
- assert_array_almost_equal(sb.dot(a,v[:,i]),w[i]*v[:,i],decimal=5)
- exact=numpy.imag(aw)
- num=numpy.imag(w)
+ assert_array_almost_equal(dot(a,v[:,i]),w[i]*v[:,i],decimal=5)
+ exact=aw.imag
+ num=w.imag
exact.sort()
num.sort()
assert_array_almost_equal(num,exact[:k],decimal=5)
- def test_type(self):
- k=2
- for typ in 'FD':
- self.large_magnitude(typ,k)
- self.small_magnitude(typ,k)
- self.large_real(typ,k)
- self.small_real(typ,k)
- self.large_imag(typ,k)
- self.small_imag(typ,k)
+# def test_type(self):
+# k=2
+# for typ in 'FD':
+# self.large_magnitude(typ,k)
+# self.small_magnitude(typ,k)
+# self.large_real(typ,k)
+# self.small_real(typ,k)
+# self.large_imag(typ,k)
+# self.small_imag(typ,k)
@@ -218,13 +211,13 @@
class TestEigenSymmetric(TestCase):
def get_a1(self,typ):
- mat_a1=numpy.array([[ 2., 0., 0., -1., 0., -1.],
- [ 0., 2., 0., -1., 0., -1.],
- [ 0., 0., 2., -1., 0., -1.],
- [-1., -1., -1., 4., 0., -1.],
- [ 0., 0., 0., 0., 1., -1.],
- [-1., -1., -1., -1., -1., 5.]],
- typ)
+ mat_a1=array([[ 2., 0., 0., -1., 0., -1.],
+ [ 0., 2., 0., -1., 0., -1.],
+ [ 0., 0., 2., -1., 0., -1.],
+ [-1., -1., -1., 4., 0., -1.],
+ [ 0., 0., 0., 0., 1., -1.],
+ [-1., -1., -1., -1., -1., 5.]],
+ typ)
w = [0,1,2,2,5,6] # eigenvalues of a1
return mat_a1,w
@@ -249,28 +242,28 @@
w,v = eigen_symmetric(a,k,which='LM')
ew,ev = eigh(a)
ind=ew.argsort()
- assert_array_almost_equal(w,numpy.take(ew,ind[-k:]))
+ assert_array_almost_equal(w,take(ew,ind[-k:]))
for i in range(k):
- assert_array_almost_equal(sb.dot(a,v[:,i]),w[i]*v[:,i])
+ assert_array_almost_equal(dot(a,v[:,i]),w[i]*v[:,i])
def small_eigenvectors(self,typ,k):
a,aw = self.get_a1(typ)
w,v = eigen_symmetric(a,k,which='SM',tol=1e-7)
ew,ev = eigh(a)
ind=ew.argsort()
- assert_array_almost_equal(w,numpy.take(ew,ind[:k]))
+ assert_array_almost_equal(w,take(ew,ind[:k]))
for i in range(k):
- assert_array_almost_equal(sb.dot(a,v[:,i]),w[i]*v[:,i])
+ assert_array_almost_equal(dot(a,v[:,i]),w[i]*v[:,i])
def end_eigenvectors(self,typ,k):
a,aw = self.get_a1(typ)
w,v = eigen_symmetric(a,k,which='BE')
ew,ev = eigh(a)
ind=ew.argsort()
- exact=numpy.concatenate(([ind[:k/2],ind[-k/2:]]))
- assert_array_almost_equal(w,numpy.take(ew,exact))
+ exact=concatenate(([ind[:k/2],ind[-k/2:]]))
+ assert_array_almost_equal(w,take(ew,exact))
for i in range(k):
- assert_array_almost_equal(sb.dot(a,v[:,i]),w[i]*v[:,i])
+ assert_array_almost_equal(dot(a,v[:,i]),w[i]*v[:,i])
def test_eigenvectors(self):
k=2
@@ -290,21 +283,21 @@
class TestEigenComplexSymmetric(TestCase):
def get_a1(self,typ):
- mat_a1=numpy.array([[ 2., 0., 0., -1., 0., -1.],
- [ 0., 2., 0., -1., 0., -1.],
- [ 0., 0., 2., -1., 0., -1.],
- [-1., -1., -1., 4., 0., -1.],
- [ 0., 0., 0., 0., 1., -1.],
- [-1., -1., -1., -1., -1., 5.]],
- typ)
- w = numpy.array([0+0j,1+0j,2+0j,2+0j,5+0j,6+0j]) # eigenvalues of a1
+ mat_a1=array([[ 2., 0., 0., -1., 0., -1.],
+ [ 0., 2., 0., -1., 0., -1.],
+ [ 0., 0., 2., -1., 0., -1.],
+ [-1., -1., -1., 4., 0., -1.],
+ [ 0., 0., 0., 0., 1., -1.],
+ [-1., -1., -1., -1., -1., 5.]],
+ typ)
+ w = array([0+0j,1+0j,2+0j,2+0j,5+0j,6+0j]) # eigenvalues of a1
return mat_a1,w
def large_magnitude(self,typ,k):
a,aw = self.get_a1(typ)
w,v = eigen(a,k,which='LM')
for i in range(k):
- assert_array_almost_equal(sb.dot(a,v[:,i]),w[i]*v[:,i],decimal=5)
+ assert_array_almost_equal(dot(a,v[:,i]),w[i]*v[:,i],decimal=5)
aw.real.sort()
w.real.sort()
assert_array_almost_equal(w,aw[-k:])
@@ -314,7 +307,7 @@
a,aw = self.get_a1(typ)
w,v = eigen(a,k,which='SM')
for i in range(k):
- assert_array_almost_equal(sb.dot(a,v[:,i]),w[i]*v[:,i])
+ assert_array_almost_equal(dot(a,v[:,i]),w[i]*v[:,i])
aw.real.sort()
w.real.sort()
assert_array_almost_equal(w,aw[:k])
@@ -323,7 +316,7 @@
a,aw = self.get_a1(typ)
w,v = eigen(a,k,which='LR')
for i in range(k):
- assert_array_almost_equal(sb.dot(a,v[:,i]),w[i]*v[:,i],decimal=5)
+ assert_array_almost_equal(dot(a,v[:,i]),w[i]*v[:,i],decimal=5)
aw.real.sort()
w.real.sort()
assert_array_almost_equal(w,aw[-k:],decimal=5)
@@ -333,18 +326,18 @@
a,aw = self.get_a1(typ)
w,v = eigen(a,k,which='SR')
for i in range(k):
- assert_array_almost_equal(sb.dot(a,v[:,i]),w[i]*v[:,i])
+ assert_array_almost_equal(dot(a,v[:,i]),w[i]*v[:,i])
aw.real.sort()
w.real.sort()
assert_array_almost_equal(w,aw[:k])
- def test_complex_symmetric(self):
- k=2
- for typ in 'FD':
- self.large_magnitude(typ,k)
- self.small_magnitude(typ,k)
- self.large_real(typ,k)
- self.small_real(typ,k)
+# def test_complex_symmetric(self):
+# k=2
+# for typ in 'FD':
+# self.large_magnitude(typ,k)
+# self.small_magnitude(typ,k)
+# self.large_real(typ,k)
+# self.small_real(typ,k)
Modified: trunk/scipy/splinalg/setup.py
===================================================================
--- trunk/scipy/splinalg/setup.py 2008-02-06 01:18:26 UTC (rev 3896)
+++ trunk/scipy/splinalg/setup.py 2008-02-06 17:56:43 UTC (rev 3897)
@@ -7,6 +7,7 @@
config.add_subpackage(('isolve'))
config.add_subpackage(('dsolve'))
+ config.add_subpackage(('eigen'))
config.add_data_dir('tests')
More information about the Scipy-svn
mailing list