[Scipy-svn] r5149 - in trunk/scipy/linalg: . tests

scipy-svn@scip... scipy-svn@scip...
Wed Nov 19 05:28:59 CST 2008


Author: tzito
Date: 2008-11-19 05:28:49 -0600 (Wed, 19 Nov 2008)
New Revision: 5149

Modified:
   trunk/scipy/linalg/decomp.py
   trunk/scipy/linalg/generic_flapack.pyf
   trunk/scipy/linalg/tests/test_decomp.py
Log:
Added new functionality and corresponding unit tests 
to linalg.eigh and linalg.eigvalsh.
This is the result of symeig integration in scipy.


Modified: trunk/scipy/linalg/decomp.py
===================================================================
--- trunk/scipy/linalg/decomp.py	2008-11-18 20:47:56 UTC (rev 5148)
+++ trunk/scipy/linalg/decomp.py	2008-11-19 11:28:49 UTC (rev 5149)
@@ -8,6 +8,7 @@
 # additions by Johannes Loehnert, June 2006
 # additions by Bart Vandereycken, June 2006
 # additions by Andrew D Straw, May 2007
+# additions by Tiziano Zito, November 2008
 
 __all__ = ['eig','eigh','eig_banded','eigvals','eigvalsh', 'eigvals_banded',
            'lu','svd','svdvals','diagsvd','cholesky','qr','qr_old','rq',
@@ -24,7 +25,7 @@
 from scipy.linalg import calc_lwork
 import numpy
 from numpy import array, asarray_chkfinite, asarray, diag, zeros, ones, \
-        single, isfinite, inexact, complexfloating
+        single, isfinite, inexact, complexfloating, nonzero, iscomplexobj
 
 cast = numpy.cast
 r_ = numpy.r_
@@ -203,97 +204,195 @@
         return w, vl
     return w, vr
 
-def eigh(a, lower=True, eigvals_only=False, overwrite_a=False):
-    """Solve the eigenvalue problem for a Hermitian or real symmetric matrix.
+def eigh(a, b=None, lower=True, eigvals_only=False, overwrite_a=False,
+         overwrite_b=False, turbo=True, eigvals=None, type=1):
+    """Solve an ordinary or generalized eigenvalue problem for a complex
+    Hermitian or real symmetric matrix. 
 
-    Find eigenvalues w and optionally right eigenvectors v of a::
+    Find eigenvalues w and optionally eigenvectors v of matrix a, where
+    b is positive definite::
 
-        a v[:,i] = w[i] v[:,i]
-        v.H v    = identity
+                      a v[:,i] = w[i] b v[:,i]
+        v[i,:].conj() a v[:,i] = w[i]
+        v[i,:].conj() b v[:,i] = 1
 
+
     Parameters
     ----------
     a : array, shape (M, M)
-        A complex Hermitian or symmetric real matrix whose eigenvalues
-        and eigenvectors will be computed.
+        A complex Hermitian or real symmetric matrix whose eigenvalues and
+        eigenvectors will be computed.
+    b : array, shape (M, M)
+        A complex Hermitian or real symmetric definite positive matrix in.
+        If omitted, identity matrix is assumed.
     lower : boolean
         Whether the pertinent array data is taken from the lower or upper
         triangle of a. (Default: lower)
     eigvals_only : boolean
         Whether to calculate only eigenvalues and no eigenvectors.
         (Default: both are calculated)
+    turbo : boolean
+        Use divide and conquer algorithm (faster but expensive in memory,
+        only for generalized eigenvalue problem and if eigvals=None)
+    eigvals : tuple (lo, hi)
+        Indexes of the smallest and largest (in ascending order) eigenvalues
+        and corresponding eigenvectors to be returned: 0 <= lo < hi <= M-1.
+        If omitted, all eigenvalues and eigenvectors are returned.
+    type: integer
+        Specifies the problem type to be solved:
+           type = 1: a   v[:,i] = w[i] b v[:,i]
+           type = 2: a b v[:,i] = w[i]   v[:,i]
+           type = 3: b a v[:,i] = w[i]   v[:,i]
     overwrite_a : boolean
-        Whether data in a is overwritten (may improve performance).
+        Whether to overwrite data in a (may improve performance)
+    overwrite_b : boolean
+        Whether to overwrite data in b (may improve performance)
 
     Returns
     -------
-    w : double array, shape (M,)
-        The eigenvalues, in ascending order, each repeated according to its
-        multiplicity.
+    w : real array, shape (N,)
+        The N (1<=N<=M) selected eigenvalues, in ascending order, each
+        repeated according to its multiplicity.
 
     (if eigvals_only == False)
-    v : double or complex double array, shape (M, M)
-        The normalized eigenvector corresponding to the eigenvalue w[i] is
-        the column v[:,i].
+    v : complex array, shape (M, N)
+        The normalized selected eigenvector corresponding to the
+        eigenvalue w[i] is the column v[:,i]. Normalization:
+        type 1 and 3:       v.conj() a      v  = w
+        type 2:        inv(v).conj() a  inv(v) = w
+        type = 1 or 2:      v.conj() b      v  = I
+        type = 3     :      v.conj() inv(b) v  = I
+        
+    Raises LinAlgError if eigenvalue computation does not converge, 
+    an error occurred, or b matrix is not definite positive. Note that
+    if input matrices are not symmetric or hermitian, no error is reported
+    but results will be wrong.
 
-    Raises LinAlgError if eigenvalue computation does not converge
-
     See Also
     --------
     eig : eigenvalues and right eigenvectors for non-symmetric arrays
 
     """
-    if eigvals_only or overwrite_a:
-        a1 = asarray_chkfinite(a)
-        overwrite_a = overwrite_a or (_datanotshared(a1,a))
-    else:
-        a1 = array(a)
-        if issubclass(a1.dtype.type, inexact) and not isfinite(a1).all():
-            raise ValueError, "array must not contain infs or NaNs"
-        overwrite_a = 1
-
+    a1 = asarray_chkfinite(a)
     if len(a1.shape) != 2 or a1.shape[0] != a1.shape[1]:
         raise ValueError, 'expected square matrix'
+    overwrite_a = overwrite_a or (_datanotshared(a1,a))
+    if iscomplexobj(a1):
+        cplx = True
+    else:
+        cplx = False
+    if b is not None:
+        b1 = asarray_chkfinite(b)
+        overwrite_b = overwrite_b or _datanotshared(b1,b)
+        if len(b1.shape) != 2 or b1.shape[0] != b1.shape[1]:
+            raise ValueError, 'expected square matrix'
 
-    if a1.dtype.char in 'FD':
-        heev, = get_lapack_funcs(('heev',),(a1,))
-        if heev.module_name[:7] == 'flapack':
-            lwork = calc_lwork.heev(heev.prefix,a1.shape[0],lower)
-            w,v,info = heev(a1,lwork = lwork,
-                            compute_v = not eigvals_only,
-                            lower = lower,
-                            overwrite_a = overwrite_a)
-        else: # 'clapack'
-            w,v,info = heev(a1,
-                            compute_v = not eigvals_only,
-                            lower = lower,
-                            overwrite_a = overwrite_a)
-        if info<0: raise ValueError,\
-           'illegal value in %-th argument of internal heev'%(-info)
-        if info>0: raise LinAlgError,"eig algorithm did not converge"
-    else: # a1.dtype.char in 'fd':
-        syev, = get_lapack_funcs(('syev',),(a1,))
-        if syev.module_name[:7] == 'flapack':
-            lwork = calc_lwork.syev(syev.prefix,a1.shape[0],lower)
-            w,v,info = syev(a1,lwork = lwork,
-                            compute_v = not eigvals_only,
-                            lower = lower,
-                            overwrite_a = overwrite_a)
-        else: # 'clapack'
-            w,v,info = syev(a1,
-                            compute_v = not eigvals_only,
-                            lower = lower,
-                            overwrite_a = overwrite_a)
-        if info<0: raise ValueError,\
-           'illegal value in %-th argument of internal syev'%(-info)
-        if info>0: raise LinAlgError,"eig algorithm did not converge"
+        if b1.shape != a1.shape:
+            raise ValueError("wrong b dimensions %s, should "
+                             "be %s" % (str(b1.shape), str(a1.shape)))
+        if iscomplexobj(b1):
+            cplx = True
+        else:
+            cplx = cplx or False
+    else:
+        b1 = None
+        
+    # Set job for fortran routines
+    _job = (eigvals_only and 'N') or 'V'
 
-    if eigvals_only:
-        return w
-    return w, v
+    # port eigenvalue range from python to fortran convention
+    if eigvals is not None:
+        lo, hi = eigvals
+        if lo < 0 or hi >= a1.shape[0]:
+            raise ValueError('The eigenvalue range specified is not valid.\n'
+                             'Valid range is [%s,%s]' % (0, a1.shape[0]-1))
+        lo += 1
+        hi += 1
+        eigvals = (lo, hi)
 
+    # set lower
+    if lower:
+        uplo = 'L'
+    else:
+        uplo = 'U'
 
+    # fix prefix for lapack routines
+    if cplx:
+        pfx = 'he'
+    else:
+        pfx = 'sy'
+        
+    #  Standard Eigenvalue Problem
+    #  Use '*evr' routines
+    # FIXME: implement calculation of optimal lwork
+    #        for all lapack routines
+    if b1 is None:
+        (evr,) = get_lapack_funcs((pfx+'evr',), (a1,))
+	if eigvals is None:
+            w, v, info = evr(a1, uplo=uplo, jobz=_job, range="A", il=1,
+                             iu=a1.shape[0], overwrite_a=overwrite_a)
+        else: 
+            (lo, hi)= eigvals
+            w_tot, v, info = evr(a1, uplo=uplo, jobz=_job, range="I",
+                                 il=lo, iu=hi, overwrite_a=overwrite_a)
+            w = w_tot[0:hi-lo+1]
 
+    # Generalized Eigenvalue Problem
+    else:
+        # Use '*gvx' routines if range is specified
+        if eigvals is not None:
+            (gvx,) = get_lapack_funcs((pfx+'gvx',), (a1,b1))
+            (lo, hi) = eigvals
+            w_tot, v, ifail, info = gvx(a1, b1, uplo=uplo, iu=hi,
+                                        itype=type,jobz=_job, il=lo,
+                                        overwrite_a=overwrite_a,
+                                        overwrite_b=overwrite_b)
+            w = w_tot[0:hi-lo+1]
+        # Use '*gvd' routine if turbo is on and no eigvals are specified
+        elif turbo:
+            (gvd,) = get_lapack_funcs((pfx+'gvd',), (a1,b1))
+            v, w, info = gvd(a1, b1, uplo=uplo, itype=type, jobz=_job,
+                             overwrite_a=overwrite_a,
+                             overwrite_b=overwrite_b)
+        # Use '*gv' routine if turbo is off and no eigvals are specified
+        else:
+            (gv,) = get_lapack_funcs((pfx+'gv',), (a1,b1))
+            v, w, info = gv(a1, b1, uplo=uplo, itype= type, jobz=_job,
+                            overwrite_a=overwrite_a,
+                            overwrite_b=overwrite_b)
+
+    # Check if we had a  successful exit
+    if info == 0:
+        if eigvals_only:
+            return w
+        else:
+            return w, v
+        
+    elif info < 0:
+        raise LinAlgError("illegal value in %i-th argument of internal"
+                          " fortran routine." % (-info))
+    elif info > 0 and b1 is None:
+        raise LinAlgError("unrecoverable internal error.")
+
+    # The algorithm failed to converge.
+    elif info > 0 and info <= b1.shape[0]:
+        if eigvals is not None:
+            raise LinAlgError("the eigenvectors %s failed to"
+                              " converge." % nonzero(ifail)-1)
+        else:
+            raise LinAlgError("internal fortran routine failed to converge: "
+                              "%i off-diagonal elements of an "
+                              "intermediate tridiagonal form did not converge"
+                              " to zero." % info)
+
+    # This occurs when b is not positive definite
+    else:
+        raise LinAlgError("the leading minor of order %i"
+                          " of 'b' is not positive definite. The"
+                          " factorization of 'b' could not be completed"
+                          " and no eigenvalues or eigenvectors were"
+                          " computed." % (info-b1.shape[0]))
+
 def eig_banded(a_band, lower=0, eigvals_only=0, overwrite_a_band=0,
                select='a', select_range=None, max_ev = 0):
     """Solve real symmetric or complex hermetian band matrix eigenvalue problem.
@@ -477,32 +576,56 @@
     """
     return eig(a,b=b,left=0,right=0,overwrite_a=overwrite_a)
 
-def eigvalsh(a,lower=1,overwrite_a=0):
-    """Solve the eigenvalue problem for a Hermitian or real symmetric matrix.
+def eigvalsh(a, b=None, lower=True, overwrite_a=False,
+             overwrite_b=False, turbo=True, eigvals=None, type=1):
+    """Solve an ordinary or generalized eigenvalue problem for a complex
+    Hermitian or real symmetric matrix. 
 
-    Find eigenvalues w of a::
+    Find eigenvalues w of matrix a, where b is positive definite::
 
-        a v[:,i] = w[i] v[:,i]
-        v.H v    = identity
+                      a v[:,i] = w[i] b v[:,i]
+        v[i,:].conj() a v[:,i] = w[i]
+        v[i,:].conj() b v[:,i] = 1
 
+
     Parameters
     ----------
     a : array, shape (M, M)
-        A complex Hermitian or symmetric real matrix whose eigenvalues
-        and eigenvectors will be computed.
+        A complex Hermitian or real symmetric matrix whose eigenvalues and
+        eigenvectors will be computed.
+    b : array, shape (M, M)
+        A complex Hermitian or real symmetric definite positive matrix in.
+        If omitted, identity matrix is assumed.
     lower : boolean
         Whether the pertinent array data is taken from the lower or upper
         triangle of a. (Default: lower)
+    turbo : boolean
+        Use divide and conquer algorithm (faster but expensive in memory,
+        only for generalized eigenvalue problem and if eigvals=None)
+    eigvals : tuple (lo, hi)
+        Indexes of the smallest and largest (in ascending order) eigenvalues
+        and corresponding eigenvectors to be returned: 0 <= lo < hi <= M-1.
+        If omitted, all eigenvalues and eigenvectors are returned.
+    type: integer
+        Specifies the problem type to be solved:
+           type = 1: a   v[:,i] = w[i] b v[:,i]
+           type = 2: a b v[:,i] = w[i]   v[:,i]
+           type = 3: b a v[:,i] = w[i]   v[:,i]
     overwrite_a : boolean
-        Whether data in a is overwritten (may improve performance).
+        Whether to overwrite data in a (may improve performance)
+    overwrite_b : boolean
+        Whether to overwrite data in b (may improve performance)
 
     Returns
     -------
-    w : double array, shape (M,)
-        The eigenvalues, in ascending order, each repeated according to its
-        multiplicity.
+    w : real array, shape (N,)
+        The N (1<=N<=M) selected eigenvalues, in ascending order, each
+        repeated according to its multiplicity.
 
-    Raises LinAlgError if eigenvalue computation does not converge
+    Raises LinAlgError if eigenvalue computation does not converge, 
+    an error occurred, or b matrix is not definite positive. Note that
+    if input matrices are not symmetric or hermitian, no error is reported
+    but results will be wrong.
 
     See Also
     --------
@@ -511,11 +634,13 @@
     eig : eigenvalues and right eigenvectors for non-symmetric arrays
 
     """
-    return eigh(a,lower=lower,eigvals_only=1,overwrite_a=overwrite_a)
+    return eigh(a, b=b, lower=lower, eigvals_only=True,
+                overwrite_a=overwrite_a, overwrite_b=overwrite_b,
+                turbo=turbo, eigvals=eigvals, type=type) 
 
 def eigvals_banded(a_band,lower=0,overwrite_a_band=0,
                    select='a', select_range=None):
-    """Solve real symmetric or complex hermetian band matrix eigenvalue problem.
+    """Solve real symmetric or complex hermitian band matrix eigenvalue problem.
 
     Find eigenvalues w of a::
 

Modified: trunk/scipy/linalg/generic_flapack.pyf
===================================================================
--- trunk/scipy/linalg/generic_flapack.pyf	2008-11-18 20:47:56 UTC (rev 5148)
+++ trunk/scipy/linalg/generic_flapack.pyf	2008-11-19 11:28:49 UTC (rev 5149)
@@ -6,7 +6,7 @@
 ! $Revision$ $Date$
 !
 ! Additions by Travis Oliphant
-! 
+! Additions by Tiziano Zito
 ! Usage:
 !   f2py -c flapack.pyf -L/usr/local/lib/atlas -llapack -lf77blas -lcblas -latlas -lg2c
 
@@ -1378,7 +1378,420 @@
     integer intent(out):: info
 end subroutine <tchar=s,d,c,z>gbtrs
 
+!
+! RRR routines for standard eigenvalue problem
+!
+subroutine ssyevr(jobz,range,uplo,n,a,lda,vl,vu,il,iu,abstol,m,w,z,ldz,isuppz,work,lwork,iwork,liwork,info)
+    ! Standard Eigenvalue Problem
+    ! simple/expert driver: all eigenvectors or optionally selected eigenvalues
+    ! algorithm: Relatively Robust Representation
+    ! matrix storage
+    ! Real - Single precision
+    character intent(in) :: jobz='V'
+    character intent(in) :: range='A'
+    character intent(in) :: uplo='L'
+    integer intent(hide) :: n=shape(a,0)
+    real intent(in,copy),dimension(n,n) :: a
+    integer intent(hide),depend(n,a) :: lda=n
+    real intent(hide) :: vl=0
+    real intent(hide) :: vu=1
+    integer optional,intent(in) :: il=1
+    integer optional,intent(in),depend(n) :: iu=n
+    real intent(hide) :: abstol=0.
+    integer  intent(hide),depend(iu) :: m=iu-il+1
+    real intent(out),dimension(n),depend(n) :: w
+    real intent(out),dimension(n,m),depend(n,m) :: z
+    integer intent(hide),check(shape(z,0)==ldz),depend(n,z) :: ldz=n
+    integer intent(hide),dimension(2*m) :: isuppz
+    integer intent(in),depend(n) :: lwork=26*n
+    real  intent(hide),dimension(lwork) :: work
+    integer intent(hide),depend(n):: liwork=10*n
+    integer intent(hide),dimension(liwork) :: iwork
+    integer  intent(out) :: info
+end subroutine ssyevr
+subroutine dsyevr(jobz,range,uplo,n,a,lda,vl,vu,il,iu,abstol,m,w,z,ldz,isuppz,work,lwork,iwork,liwork,info)
+    ! Standard Eigenvalue Problem
+    ! simple/expert driver: all eigenvectors or optionally selected eigenvalues
+    ! algorithm: Relatively Robust Representation
+    ! matrix storage
+    ! Real - Double precision
+    character intent(in) :: jobz='V'
+    character intent(in) :: range='A'
+    character intent(in) :: uplo='L'
+    integer intent(hide) :: n=shape(a,0)
+    double precision intent(in,copy),dimension(n,n) :: a
+    integer intent(hide),depend(n,a) :: lda=n
+    double precision intent(hide) :: vl=0
+    double precision intent(hide) :: vu=1
+    integer optional,intent(in) :: il=1
+    integer optional,intent(in),depend(n) :: iu=n
+    double precision intent(hide) :: abstol=0.
+    integer  intent(hide),depend(iu) :: m=iu-il+1
+    double precision intent(out),dimension(n),depend(n) :: w
+    double precision intent(out),dimension(n,m),depend(n,m) :: z
+    integer intent(hide),check(shape(z,0)==ldz),depend(n,z) :: ldz=n
+    integer intent(hide),dimension(2*m) :: isuppz
+    integer intent(in),depend(n) :: lwork=26*n
+    double precision intent(hide),dimension(lwork) :: work
+    integer intent(hide),depend(n):: liwork=10*n
+    integer intent(hide),dimension(liwork) :: iwork
+    integer  intent(out) :: info
+end subroutine dsyevr
+subroutine cheevr(jobz,range,uplo,n,a,lda,vl,vu,il,iu,abstol,m,w,z,ldz,isuppz,work,lwork,rwork,lrwork,iwork,liwork,info)
+    ! Standard Eigenvalue Problem
+    ! simple/expert driver: all eigenvectors or optionally selected eigenvalues
+    ! algorithm: Relatively Robust Representation
+    ! matrix storage
+    ! Complex - Single precision
+    character intent(in) :: jobz='V'
+    character intent(in) :: range='A'
+    character intent(in) :: uplo='L'
+    integer intent(hide) :: n=shape(a,0)
+    complex intent(in,copy),dimension(n,n) :: a
+    integer intent(hide),depend(n,a) :: lda=n
+    real intent(hide) :: vl=0
+    real intent(hide) :: vu=1
+    integer optional,intent(in) :: il=1
+    integer optional,intent(in),depend(n) :: iu=n
+    real intent(hide) :: abstol=0.
+    integer  intent(hide),depend(iu) :: m=iu-il+1
+    real intent(out),dimension(n),depend(n) :: w
+    complex intent(out),dimension(n,m),depend(n,m) :: z
+    integer intent(hide),check(shape(z,0)==ldz),depend(n,z) :: ldz=n
+    integer intent(hide),dimension(2*m) :: isuppz
+    integer intent(in),depend(n) :: lwork=2*n
+    complex  intent(hide),dimension(lwork) :: work
+    integer intent(hide),depend(n) :: lrwork=24*n
+    real  intent(hide),dimension(lrwork) :: rwork
+    integer intent(hide),depend(n):: liwork=10*n
+    integer intent(hide),dimension(liwork) :: iwork
+    integer  intent(out) :: info
+end subroutine cheevr
+subroutine zheevr(jobz,range,uplo,n,a,lda,vl,vu,il,iu,abstol,m,w,z,ldz,isuppz,work,lwork,rwork,lrwork,iwork,liwork,info)
+    ! Standard Eigenvalue Problem
+    ! simple/expert driver: all eigenvectors or optionally selected eigenvalues
+    ! algorithm: Relatively Robust Representation
+    ! matrix storage
+    ! Complex - Double precision
+    character intent(in) :: jobz='V'
+    character intent(in) :: range='A'
+    character intent(in) :: uplo='L'
+    integer intent(hide) :: n=shape(a,0)
+    complex*16 intent(in,copy),dimension(n,n) :: a
+    integer intent(hide),depend(n,a) :: lda=n
+    double precision intent(hide) :: vl=0
+    double precision intent(hide) :: vu=1
+    integer optional,intent(in) :: il=1
+    integer optional,intent(in),depend(n) :: iu=n
+    double precision intent(hide) :: abstol=0.
+    integer  intent(hide),depend(iu) :: m=iu-il+1
+    double precision  intent(out),dimension(n),depend(n) :: w
+    complex*16 intent(out),dimension(n,m),depend(n,m) :: z
+    integer intent(hide),check(shape(z,0)==ldz),depend(n,z) :: ldz=n
+    integer intent(hide),dimension(2*m) :: isuppz
+    integer intent(in),depend(n) :: lwork=2*n
+    complex*16  intent(hide),dimension(lwork) :: work
+    integer intent(hide),depend(n) :: lrwork=24*n
+    double precision  intent(hide),dimension(lrwork) :: rwork
+    integer intent(hide),depend(n):: liwork=10*n
+    integer intent(hide),dimension(liwork) :: iwork
+    integer  intent(out) :: info
+end subroutine zheevr
+subroutine ssygv(itype,jobz,uplo,n,a,lda,b,ldb,w,work,lwork,info)
+    ! Generalized Eigenvalue Problem
+    ! simple driver (all eigenvectors)
+    ! algorithm: standard
+    ! matrix storage
+    ! Real - Single precision
+    integer optional,intent(in) :: itype=1
+    character intent(in) :: jobz='V'
+    character intent(in) :: uplo='L'
+    integer intent(hide) :: n=shape(a,0)
+    real intent(in,copy,out),dimension(n,n) :: a
+    integer intent(hide),depend(n,a) :: lda=n
+    real intent(in,copy),dimension(n,n) :: b
+    integer intent(hide),depend(n,b) :: ldb=n
+    real intent(out),dimension(n),depend(n) :: w
+    integer intent(hide) :: lwork=3*n-1
+    real intent(hide),dimension(lwork) :: work
+    integer intent(out) :: info
+end subroutine ssygv
+subroutine dsygv(itype,jobz,uplo,n,a,lda,b,ldb,w,work,lwork,info)
+    ! Generalized Eigenvalue Problem
+    ! simple driver (all eigenvectors)
+    ! algorithm: standard
+    ! matrix storage
+    ! Real - Double precision
+    integer optional,intent(in) :: itype=1
+    character intent(in) :: jobz='V'
+    character intent(in) :: uplo='L'
+    integer intent(hide) :: n=shape(a,0)
+    double precision intent(in,copy,out),dimension(n,n) :: a
+    integer intent(hide),depend(n,a) :: lda=n
+    double precision intent(in,copy),dimension(n,n) :: b
+    integer intent(hide),depend(n,b) :: ldb=n
+    double precision intent(out),dimension(n),depend(n) :: w
+    integer intent(hide) :: lwork=3*n-1
+    double precision intent(hide),dimension(lwork) :: work
+    integer intent(out) :: info
+end subroutine dsygv
+subroutine chegv(itype,jobz,uplo,n,a,lda,b,ldb,w,work,lwork,rwork,info)
+    ! Generalized Eigenvalue Problem
+    ! simple driver (all eigenvectors)
+    ! algorithm: standard
+    ! matrix storage
+    ! Complex - Single precision
+    integer optional,intent(in) :: itype=1
+    character intent(in) :: jobz='V'
+    character intent(in) :: uplo='L'
+    integer intent(hide) :: n=shape(a,0)
+    complex intent(in,copy,out),dimension(n,n) :: a
+    integer intent(hide),depend(n,a) :: lda=n
+    complex intent(in,copy),dimension(n,n) :: b
+    integer intent(hide),depend(n,b) :: ldb=n
+    real intent(out),dimension(n),depend(n) :: w
+    integer intent(hide) :: lwork=2*n-1
+    complex intent(hide),dimension(lwork) :: work
+    real intent(hide),dimension(3*n-2) :: rwork
+    integer intent(out) :: info
+end subroutine chegv
+subroutine zhegv(itype,jobz,uplo,n,a,lda,b,ldb,w,work,lwork,rwork,info)
+    ! Generalized Eigenvalue Problem
+    ! simple driver (all eigenvectors)
+    ! algorithm: standard
+    ! matrix storage
+    ! Complex - Double precision
+    integer optional,intent(in) :: itype=1
+    character intent(in) :: jobz='V'
+    character intent(in) :: uplo='L'
+    integer intent(hide) :: n=shape(a,0)
+    complex*16 intent(in,copy,out),dimension(n,n) :: a
+    integer intent(hide),depend(n,a) :: lda=n
+    complex*16 intent(in,copy),dimension(n,n) :: b
+    integer intent(hide),depend(n,b) :: ldb=n
+    double precision intent(out),dimension(n),depend(n) :: w
+    integer intent(hide) :: lwork=2*n-1
+    complex*16 intent(hide),dimension(lwork) :: work
+    double precision intent(hide),dimension(3*n-2) :: rwork
+    integer intent(out) :: info
+end subroutine zhegv
+!
+! Divide and conquer routines for generalized eigenvalue problem
+!
+subroutine ssygvd(itype,jobz,uplo,n,a,lda,b,ldb,w,work,lwork,iwork,liwork,info)
+    ! Generalized Eigenvalue Problem
+    ! simple driver (all eigenvectors)
+    ! algorithm: divide and conquer
+    ! matrix storage
+    ! Real - Single precision
+    integer optional,intent(in) :: itype=1
+    character intent(in) :: jobz='V'
+    character intent(in) :: uplo='L'
+    integer intent(hide) :: n=shape(a,0)
+    real intent(in,copy,out),dimension(n,n) :: a
+    integer intent(hide),depend(n,a) :: lda=n
+    real intent(in,copy),dimension(n,n) :: b
+    integer intent(hide),depend(n,b) :: ldb=n
+    real intent(out),dimension(n),depend(n) :: w
+    integer intent(in),depend(n) :: lwork=1+6*n+2*n*n
+    real intent(hide),dimension(lwork) :: work
+    integer intent(hide),depend(n) :: liwork=3+5*n
+    integer intent(hide),dimension(liwork) :: iwork
+    integer intent(out) :: info
+end subroutine ssygvd
+subroutine dsygvd(itype,jobz,uplo,n,a,lda,b,ldb,w,work,lwork,iwork,liwork,info)
+    ! Generalized Eigenvalue Problem
+    ! simple driver (all eigenvectors)
+    ! algorithm: divide and conquer
+    ! matrix storage
+    ! Real - Double precision
+    integer optional,intent(in) :: itype=1
+    character intent(in) :: jobz='V'
+    character intent(in) :: uplo='L'
+    integer intent(hide) :: n=shape(a,0)
+    double precision intent(in,copy,out),dimension(n,n) :: a
+    integer intent(hide),depend(n,a) :: lda=n
+    double precision intent(in,copy),dimension(n,n) :: b
+    integer intent(hide),depend(n,b) :: ldb=n
+    double precision intent(out),dimension(n),depend(n) :: w
+    integer intent(in),depend(n) :: lwork=1+6*n+2*n*n
+    double precision intent(hide),dimension(lwork) :: work
+    integer intent(hide),depend(n) :: liwork=3+5*n
+    integer intent(hide),dimension(liwork) :: iwork
+    integer intent(out) :: info
+end subroutine dsygvd
+subroutine chegvd(itype,jobz,uplo,n,a,lda,b,ldb,w,work,lwork,rwork,lrwork,iwork,liwork,info)
+    ! Generalized Eigenvalue Problem
+    ! simple driver (all eigenvectors)
+    ! algorithm: divide and conquer
+    ! matrix storage
+    ! Complex - Single precision
+    integer optional,intent(in) :: itype=1
+    character intent(in) :: jobz='V'
+    character intent(in) :: uplo='L'
+    integer intent(hide) :: n=shape(a,0)
+    complex intent(in,copy,out),dimension(n,n) :: a
+    integer intent(hide),depend(n,a) :: lda=n
+    complex intent(in,copy),dimension(n,n) :: b
+    integer intent(hide),depend(n,b) :: ldb=n
+    real intent(out),dimension(n),depend(n) :: w
+    integer intent(in),depend(n) :: lwork=2*n+n*n
+    complex intent(hide),dimension(lwork) :: work
+    integer intent(hide),depend(n) :: lrwork=1+5*n+2*n*n
+    real intent(hide),dimension(lrwork) :: rwork
+    integer intent(hide),depend(n) :: liwork=3+5*n
+    integer intent(hide),dimension(liwork) :: iwork
+    integer intent(out) :: info
+end subroutine chegvd
+subroutine zhegvd(itype,jobz,uplo,n,a,lda,b,ldb,w,work,lwork,rwork,lrwork,iwork,liwork,info)
+    ! Generalized Eigenvalue Problem
+    ! simple driver (all eigenvectors)
+    ! algorithm: divide and conquer
+    ! matrix storage
+    ! Complex - Double precision
+    integer optional,intent(in) :: itype=1
+    character intent(in) :: jobz='V'
+    character intent(in) :: uplo='L'
+    integer intent(hide) :: n=shape(a,0)
+    complex*16 intent(in,copy,out),dimension(n,n) :: a
+    integer intent(hide),depend(n,a) :: lda=n
+    complex*16 intent(in,copy),dimension(n,n) :: b
+    integer intent(hide),depend(n,b) :: ldb=n
+    double precision intent(out),dimension(n),depend(n) :: w
+    integer intent(in),depend(n) :: lwork=2*n+n*n
+    complex*16 intent(hide),dimension(lwork) :: work
+    integer intent(hide),depend(n) :: lrwork=1+5*n+2*n*n
+    double precision intent(hide),dimension(lrwork) :: rwork
+    integer intent(hide),depend(n) :: liwork=3+5*n
+    integer intent(hide),dimension(liwork) :: iwork
+    integer intent(out) :: info
+end subroutine zhegvd
+! Expert routines for generalized eigenvalue problem
+!
+subroutine ssygvx(itype,jobz,range,uplo,n,a,lda,b,ldb,vl,vu,il,iu,abstol,m,w,z,ldz,work,lwork,iwork,ifail,info)
+    ! Generalized Eigenvalue Problem
+    ! expert driver (selected eigenvectors)
+    ! algorithm: standard
+    ! matrix storage
+    ! Real - Single precision
+    integer optional,intent(in) :: itype=1
+    character intent(in) :: jobz='V'
+    character intent(hide) :: range='I'
+    character intent(in) :: uplo='L'
+    integer intent(hide) :: n=shape(a,0)
+    real intent(in,copy),dimension(n,n) :: a
+    integer intent(hide),depend(n,a) :: lda=n
+    real intent(in,copy),dimension(n,n) :: b
+    integer intent(hide),depend(n,b) :: ldb=n
+    real intent(hide) :: vl=0.
+    real intent(hide) :: vu=0.
+    integer optional,intent(in) :: il=1
+    integer intent(in) :: iu
+    real intent(hide) :: abstol=0.
+    integer intent(hide),depend(iu) :: m=iu-il+1
+    real intent(out),dimension(n),depend(n) :: w
+    real intent(out),dimension(n,m),depend(n,m) :: z
+    integer intent(hide),check(shape(z,0)==ldz),depend(n,z) :: ldz=n
+    integer intent(in),depend(n) :: lwork=8*n
+    real intent(hide),dimension(lwork),depend(n,lwork) :: work
+    integer intent(hide),dimension(5*n) :: iwork
+    integer intent(out),dimension(n),depend(n) :: ifail
+    integer intent(out) :: info
+end subroutine ssygvx
+subroutine dsygvx(itype,jobz,range,uplo,n,a,lda,b,ldb,vl,vu,il,iu,abstol,m,w,z,ldz,work,lwork,iwork,ifail,info)
+    ! Generalized Eigenvalue Problem
+    ! expert driver (selected eigenvectors)
+    ! algorithm: standard
+    ! matrix storage
+    ! Real - Double precision
+    integer optional,intent(in) :: itype=1
+    character intent(in) :: jobz='V'
+    character intent(hide) :: range='I'
+    character intent(in) :: uplo='L'
+    integer intent(hide) :: n=shape(a,0)
+    double precision intent(in,copy),dimension(n,n) :: a
+    integer intent(hide),depend(n,a) :: lda=n
+    double precision intent(in,copy),dimension(n,n) :: b
+    integer intent(hide),depend(n,b) :: ldb=n
+    double precision intent(hide) :: vl=0.
+    double precision intent(hide) :: vu=0.
+    integer optional,intent(in) :: il=1
+    integer intent(in) :: iu
+    double precision intent(hide) :: abstol=0.
+    integer intent(hide),depend(iu) :: m=iu-il+1
+    double precision intent(out),dimension(n),depend(n) :: w
+    double precision intent(out),dimension(n,m),depend(n,m) :: z
+    integer intent(hide),check(shape(z,0)==ldz),depend(n,z) :: ldz=n
+    integer intent(in),depend(n) :: lwork=8*n
+    double precision intent(hide),dimension(lwork),depend(n,lwork) :: work
+    integer intent(hide),dimension(5*n) :: iwork
+    integer intent(out),dimension(n),depend(n) :: ifail
+    integer intent(out) :: info
+end subroutine dsygvx
+subroutine chegvx(itype,jobz,range,uplo,n,a,lda,b,ldb,vl,vu,il,iu,abstol,m,w,z,ldz,work,lwork,rwork,iwork,ifail,info)
+    ! Generalized Eigenvalue Problem
+    ! expert driver (selected eigenvectors)
+    ! algorithm: standard
+    ! matrix storage
+    ! Complex - Single precision
+    integer optional,intent(in) :: itype=1
+    character intent(in) :: jobz='V'
+    character intent(hide) :: range='I'
+    character intent(in) :: uplo='L'
+    integer intent(hide) :: n=shape(a,0)
+    complex intent(in,copy),dimension(n,n) :: a
+    integer intent(hide),depend(n,a) :: lda=n
+    complex intent(in,copy),dimension(n,n) :: b
+    integer intent(hide),depend(n,b) :: ldb=n
+    real intent(hide) :: vl=0.
+    real intent(hide) :: vu=0.
+    integer optional,intent(in) :: il=1
+    integer intent(in) :: iu
+    real intent(hide) :: abstol=0.
+    integer intent(hide),depend(iu) :: m=iu-il+1
+    real intent(out),dimension(n),depend(n) :: w
+    complex intent(out),dimension(n,m),depend(n,m) :: z
+    integer intent(hide),check(shape(z,0)==ldz),depend(n,z) :: ldz=n
+    integer intent(in),depend(n) :: lwork=2*n
+    complex intent(hide),dimension(lwork),depend(n,lwork) :: work
+    real intent(hide),dimension(7*n) :: rwork
+    integer intent(hide),dimension(5*n) :: iwork
+    integer intent(out),dimension(n),depend(n) :: ifail
+    integer intent(out) :: info
+end subroutine chegvx
+subroutine zhegvx(itype,jobz,range,uplo,n,a,lda,b,ldb,vl,vu,il,iu,abstol,m,w,z,ldz,work,lwork,rwork,iwork,ifail,info)
+    ! Generalized Eigenvalue Problem
+    ! expert driver (selected eigenvectors)
+    ! algorithm: standard
+    ! matrix storage
+    ! Complex - Double precision
+    integer optional,intent(in) :: itype=1
+    character intent(in) :: jobz='V'
+    character intent(hide) :: range='I'
+    character intent(in) :: uplo='L'
+    integer intent(hide) :: n=shape(a,0)
+    complex*16 intent(in,copy),dimension(n,n) :: a
+    integer intent(hide),depend(n,a) :: lda=n
+    complex*16 intent(in,copy),dimension(n,n) :: b
+    integer intent(hide),depend(n,b) :: ldb=n
+    double precision intent(hide) :: vl=0.
+    double precision intent(hide) :: vu=0.
+    integer optional,intent(in) :: il=1
+    integer intent(in) :: iu
+    double precision intent(hide) :: abstol=0.
+    integer intent(hide),depend(iu) :: m=iu-il+1
+    double precision intent(out),dimension(n),depend(n) :: w
+    complex*16 intent(out),dimension(n,m),depend(n,m) :: z
+    integer intent(hide),check(shape(z,0)==ldz),depend(n,z) :: ldz=n
+    integer intent(in),depend(n) :: lwork=2*n
+    complex*16 intent(hide),dimension(lwork),depend(n,lwork) :: work
+    double precision intent(hide),dimension(7*n) :: rwork
+    integer intent(hide),dimension(5*n) :: iwork
+    integer intent(out),dimension(n),depend(n) :: ifail
+    integer intent(out) :: info
+end subroutine zhegvx
 
+
 end interface
 
 end python module flapack

Modified: trunk/scipy/linalg/tests/test_decomp.py
===================================================================
--- trunk/scipy/linalg/tests/test_decomp.py	2008-11-18 20:47:56 UTC (rev 5148)
+++ trunk/scipy/linalg/tests/test_decomp.py	2008-11-19 11:28:49 UTC (rev 5149)
@@ -25,16 +25,14 @@
 from numpy import array, transpose, sometrue, diag, ones, linalg, \
      argsort, zeros, arange, float32, complex64, dot, conj, identity, \
      ravel, sqrt, iscomplex, shape, sort, conjugate, bmat, sign, \
-     asarray, matrix, isfinite, all, ndarray, outer, eye, dtype
+     asarray, matrix, isfinite, all, ndarray, outer, eye, dtype, empty,\
+     triu, tril
 
 from numpy.random import rand, normal
 
 # digit precision to use in asserts for different types
-DIGITS = {'d':12, 'D':12, 'f':6, 'F':6}
+DIGITS = {'d':11, 'D':11, 'f':4, 'F':4}
 
-# matrix dimension in tests
-DIM = 5
-
 # XXX: This function should be available through numpy.testing
 def assert_dtype_equal(act, des):
     if isinstance(act, ndarray):
@@ -51,24 +49,18 @@
 
 # XXX: This function should not be defined here, but somewhere in
 #      scipy.linalg namespace
-def hermitian(x):
-    """Return the Hermitian, i.e. conjugate transpose, of x."""
-    return x.T.conj()
-
-# XXX: This function should not be defined here, but somewhere in
-#      scipy.linalg namespace
-def symrand(dim_or_eigv, dtype="d"):
+def symrand(dim_or_eigv):
     """Return a random symmetric (Hermitian) matrix.
 
     If 'dim_or_eigv' is an integer N, return a NxN matrix, with eigenvalues
-        uniformly distributed on (0.1,1].
+        uniformly distributed on (-1,1).
 
     If 'dim_or_eigv' is  1-D real array 'a', return a matrix whose
-                      eigenvalues are sort(a).
+                      eigenvalues are 'a'.
     """
     if isinstance(dim_or_eigv, int):
         dim = dim_or_eigv
-        d = (rand(dim)*0.9)+0.1
+        d = (rand(dim)*2)-1
     elif (isinstance(dim_or_eigv, ndarray) and
           len(dim_or_eigv.shape) == 1):
         dim = dim_or_eigv.shape[0]
@@ -76,14 +68,15 @@
     else:
         raise TypeError("input type not supported.")
 
-    v = random_rot(dim, dtype=dtype)
-    h = dot(dot(hermitian(v), diag(d)), v)
+    v = random_rot(dim)
+    h = dot(dot(v.T.conj(), diag(d)), v)
     # to avoid roundoff errors, symmetrize the matrix (again)
-    return (0.5*(hermitian(h)+h)).astype(dtype)
+    h = 0.5*(h.T+h)
+    return h
 
 # XXX: This function should not be defined here, but somewhere in
 #      scipy.linalg namespace
-def random_rot(dim, dtype='d'):
+def random_rot(dim):
     """Return a random rotation matrix, drawn from the Haar distribution
     (the only uniform distribution on SO(n)).
     The algorithm is described in the paper
@@ -92,16 +85,16 @@
     on Numerical Analysis, 17(3), pp. 403-409, 1980.
     For more information see
     http://en.wikipedia.org/wiki/Orthogonal_matrix#Randomization"""
-    H = eye(dim, dtype=dtype)
-    D = ones((dim, ), dtype=dtype)
+    H = eye(dim)
+    D = ones((dim, ))
     for n in range(1, dim):
-        x = normal(size=(dim-n+1, )).astype(dtype)
+        x = normal(size=(dim-n+1, ))
         D[n-1] = sign(x[0])
         x[0] -= D[n-1]*sqrt((x*x).sum())
         # Householder transformation
 
-        Hx = eye(dim-n+1, dtype=dtype) - 2.*outer(x, x)/(x*x).sum()
-        mat = eye(dim, dtype=dtype)
+        Hx = eye(dim-n+1) - 2.*outer(x, x)/(x*x).sum()
+        mat = eye(dim)
         mat[n-1:,n-1:] = Hx
         H = dot(H, mat)
     # Fix the last sign such that the determinant is 1
@@ -493,37 +486,86 @@
         y_lin = linalg.solve(self.comp_mat, self.bc)
         assert_array_almost_equal(y, y_lin)
 
-class TestEigH(TestCase):
-    def eigenproblem_standard(self, dim, dtype, overwrite, lower):
-        """Solve a standard eigenvalue problem."""
-        a = symrand(dim, dtype)
-        if overwrite:
-            a_c = a.copy()
-        else:
-            a_c = a
-        w, z = eigh(a, lower=lower, overwrite_a = overwrite)
-        assert_dtype_equal(z.dtype, dtype)
-        w = w.astype(dtype)
-        diag_ = diag(dot(hermitian(z), dot(a_c, z))).real
-        assert_array_almost_equal(diag_, w, DIGITS[dtype])
+def test_eigh():
+    DIM = 6
+    v = {'dim': (DIM, ),
+         'dtype': ('f','d','F','D'),
+         'overwrite': (True, False),
+         'lower': (True, False),
+         'turbo': (True, False),
+         'eigvals': (None, (2, DIM-2))}
 
-    def test_eigh_real_standard(self):
-        self.eigenproblem_standard(DIM, 'd', False, False)
-        self.eigenproblem_standard(DIM, 'd', False, True)
-        self.eigenproblem_standard(DIM, 'd', True, True)
-        self.eigenproblem_standard(DIM, 'f', False, False)
-        self.eigenproblem_standard(DIM, 'f', False, True)
-        self.eigenproblem_standard(DIM, 'f', True, True)
+    for dim in v['dim']:
+        for typ in v['dtype']:
+            for overwrite in v['overwrite']:
+                for turbo in v['turbo']:
+                    for eigvals in v['eigvals']:
+                        for lower in v['lower']:
+                            yield (eigenhproblem_standard,
+                                   'ordinary',
+                                   dim, typ, overwrite, lower,
+                                   turbo, eigvals)
+                            yield (eigenhproblem_general,
+                                   'general ',
+                                   dim, typ, overwrite, lower,
+                                   turbo, eigvals)
 
-    def test_eigh_complex_standard(self):
-        self.eigenproblem_standard(DIM, 'D', False, False)
-        self.eigenproblem_standard(DIM, 'D', False, True)
-        self.eigenproblem_standard(DIM, 'D', True, True)
-        self.eigenproblem_standard(DIM, 'F', False, False)
-        self.eigenproblem_standard(DIM, 'F', False, True)
-        self.eigenproblem_standard(DIM, 'F', True, True)
+def _complex_symrand(dim, dtype):
+    a1, a2 = symrand(dim), symrand(dim)
+    # add antisymmetric matrix as imag part
+    a = a1 +1j*(triu(a2)-tril(a2))
+    return a.astype(dtype)
+  
+def eigenhproblem_standard(desc, dim, dtype,
+                           overwrite, lower, turbo,
+                           eigvals):
+    """Solve a standard eigenvalue problem."""
+    if iscomplex(empty(1, dtype=dtype)):
+        a = _complex_symrand(dim, dtype)
+    else:
+        a = symrand(dim).astype(dtype)
+    
+    if overwrite:
+        a_c = a.copy()
+    else:
+        a_c = a
+    w, z = eigh(a, overwrite_a=overwrite, lower=lower, eigvals=eigvals)
+    assert_dtype_equal(z.dtype, dtype)
+    w = w.astype(dtype)
+    diag_ = diag(dot(z.T.conj(), dot(a_c, z))).real
+    assert_array_almost_equal(diag_, w, DIGITS[dtype])
 
+def eigenhproblem_general(desc, dim, dtype,
+                          overwrite, lower, turbo,
+                          eigvals):
+    """Solve a generalized eigenvalue problem."""
+    if iscomplex(empty(1, dtype=dtype)):
+        a = _complex_symrand(dim, dtype)
+        b = _complex_symrand(dim, dtype)+diag([2.1]*dim).astype(dtype)
+    else:
+        a = symrand(dim).astype(dtype)
+        b = symrand(dim).astype(dtype)+diag([2.1]*dim).astype(dtype)
 
+    if overwrite:
+        a_c, b_c = a.copy(), b.copy()
+    else:
+        a_c, b_c = a, b
+
+    w, z = eigh(a, b, overwrite_a=overwrite, lower=lower,
+                overwrite_b=overwrite, turbo=turbo, eigvals=eigvals)
+    assert_dtype_equal(z.dtype, dtype)
+    w = w.astype(dtype)
+    diag1_ = diag(dot(z.T.conj(), dot(a_c, z))).real
+    assert_array_almost_equal(diag1_, w, DIGITS[dtype])
+    diag2_ = diag(dot(z.T.conj(), dot(b_c, z))).real
+    assert_array_almost_equal(diag2_, ones(diag2_.shape[0]), DIGITS[dtype])
+    
+def test_eigh_integer():
+    a = array([[1,2],[2,7]])
+    b = array([[3,1],[1,5]])
+    w,z = eigh(a)
+    w,z = eigh(a,b)
+
 class TestLU(TestCase):
 
     def __init__(self, *args, **kw):



More information about the Scipy-svn mailing list