[Scipy-svn] r3897 - in trunk/scipy: sandbox splinalg splinalg/eigen splinalg/eigen/arpack splinalg/eigen/arpack/tests

scipy-svn@scip... scipy-svn@scip...
Wed Feb 6 11:57:19 CST 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