[Scipy-svn] r2030 - in trunk/Lib/linalg: . tests

scipy-svn at scipy.org scipy-svn at scipy.org
Fri Jun 30 06:38:14 CDT 2006


Author: abaecker
Date: 2006-06-30 06:38:08 -0500 (Fri, 30 Jun 2006)
New Revision: 2030

Modified:
   trunk/Lib/linalg/decomp.py
   trunk/Lib/linalg/generic_flapack.pyf
   trunk/Lib/linalg/tests/test_decomp.py
Log:
band matrix wrappers added


Modified: trunk/Lib/linalg/decomp.py
===================================================================
--- trunk/Lib/linalg/decomp.py	2006-06-30 08:30:00 UTC (rev 2029)
+++ trunk/Lib/linalg/decomp.py	2006-06-30 11:38:08 UTC (rev 2030)
@@ -5,9 +5,11 @@
 #
 # additions by Travis Oliphant, March 2002
 # additions by Eric Jones,      June 2002
+# additions by Johannes Loehnert, June 2006
 
-__all__ = ['eig','eigh','eigvals','eigvalsh','lu','svd','svdvals','diagsvd',
-           'cholesky','qr',
+
+__all__ = ['eig','eigh','eig_banded','eigvals','eigvalsh', 'eigvals_banded',
+           'lu','svd','svdvals','diagsvd','cholesky','qr',
            'schur','rsf2csf','lu_factor','cho_factor','cho_solve','orth',
            'hessenberg']
 
@@ -21,7 +23,7 @@
 import calc_lwork
 import numpy
 from numpy import array, asarray_chkfinite, asarray, diag, zeros, ones, \
-        single, isfinite
+        single, isfinite, inexact
 
 cast = numpy.cast
 r_ = numpy.r_
@@ -245,6 +247,124 @@
     return w, v
 
 
+    
+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 problem.
+
+    Inputs:
+
+      a_band            -- A hermitian N x M matrix in 'packed storage'.
+                           Packed storage looks like this: ('upper form')
+                           [ ... (more off-diagonals) ...,
+                            [ *   *   a13 a24 a35 a46 ... a(n-2)(n)],
+                            [ *   a12 a23 a34 a45 a56 ... a(n-1)(n)],
+                            [ a11 a22 a33 a44 a55 a66 ... a(n)(n)  ]]
+                           The cells denoted with * may contain anything.
+      lower             -- a is in lower packed storage
+                           (default: upper packed form)
+      eigvals_only      -- if True, don't compute eigenvectors.
+      overwrite_a_band  -- content of a may be destroyed
+      select       -- 'a', 'all', 0   : return all eigenvalues/vectors
+                      'v', 'value', 1 : eigenvalues in the interval (min, max]
+                                        will be found
+                      'i', 'index', 2 : eigenvalues with index [min...max]
+                                        will be found
+      select_range -- select == 'v': eigenvalue limits as tuple (min, max)
+                      select == 'i': index limits as tuple (min, max)
+                      select == 'a': meaningless
+      max_ev       -- select == 'v': set to max. number of eigenvalues that is
+                                     expected. In doubt, leave untouched.
+                      select == 'i', 'a': meaningless
+
+    Outputs:
+
+      w,v     -- w: eigenvalues, v: eigenvectors [for eigvals_only == False]
+      w       -- eigenvalues [for eigvals_only == True].
+
+    Definitions:
+
+      a_full * v[:,i] = w[i] * v[:,i]  (with full matrix corresponding to a)
+      v.H * v = identity
+
+    """
+    if eigvals_only or overwrite_a_band:
+        a1 = asarray_chkfinite(a_band)
+        overwrite_a_band = overwrite_a_band or (_datanotshared(a1,a_band))
+    else:
+        a1 = array(a_band)
+        if issubclass(a1.dtype.type, inexact) and not isfinite(a1).all():
+            raise ValueError, "array must not contain infs or NaNs"
+        overwrite_a_band = 1
+
+    if len(a1.shape) != 2:
+        raise ValueError, 'expected two-dimensional array'
+    if select.lower() not in [0, 1, 2, 'a', 'v', 'i', 'all', 'value', 'index']:
+        raise ValueError, 'invalid argument for select'
+    if select.lower() in [0, 'a', 'all']:
+        if a1.dtype.char in 'GFD':
+            bevd, = get_lapack_funcs(('hbevd',),(a1,))
+            # FIXME: implement this somewhen, for now go with builtin values
+            # FIXME: calc optimal lwork by calling ?hbevd(lwork=-1)
+            #        or by using calc_lwork.f ???
+            # lwork = calc_lwork.hbevd(bevd.prefix, a1.shape[0], lower)
+            internal_name = 'hbevd'
+        else: # a1.dtype.char in 'fd':
+            bevd, = get_lapack_funcs(('sbevd',),(a1,))
+            # FIXME: implement this somewhen, for now go with builtin values
+            #         see above
+            # lwork = calc_lwork.sbevd(bevd.prefix, a1.shape[0], lower)
+            internal_name = 'sbevd'
+        w,v,info = bevd(a1, compute_v = not eigvals_only,
+                        lower = lower,
+                        overwrite_ab = overwrite_a_band)
+    if select.lower() in [1, 2, 'i', 'v', 'index', 'value']:
+        # calculate certain range only
+        if select.lower() in [2, 'i', 'index']:
+            select = 2
+            vl, vu, il, iu = 0.0, 0.0, min(select_range), max(select_range)
+            if min(il, iu) < 0 or max(il, iu) >= a1.shape[1]:
+                raise ValueError, 'select_range out of bounds'
+            max_ev = iu - il + 1
+        else:  # 1, 'v', 'value'
+            select = 1
+            vl, vu, il, iu = min(select_range), max(select_range), 0, 0
+            if max_ev == 0:
+                max_ev = a_band.shape[1]
+        if eigvals_only:
+            max_ev = 1
+        # calculate optimal abstol for dsbevx (see manpage)
+        if a1.dtype.char in 'fF':  # single precision
+            lamch, = get_lapack_funcs(('lamch',),(array(0, dtype='f'),))
+        else:
+            lamch, = get_lapack_funcs(('lamch',),(array(0, dtype='d'),))
+        abstol = 2 * lamch('s')
+        if a1.dtype.char in 'GFD':
+            bevx, = get_lapack_funcs(('hbevx',),(a1,))
+            internal_name = 'hbevx'
+        else: # a1.dtype.char in 'gfd'
+            bevx, = get_lapack_funcs(('sbevx',),(a1,))
+            internal_name = 'sbevx'
+        # il+1, iu+1: translate python indexing (0 ... N-1) into Fortran
+        # indexing (1 ... N)
+        w, v, m, ifail, info = bevx(a1, vl, vu, il+1, iu+1,
+                                    compute_v = not eigvals_only,
+                                    mmax = max_ev,
+                                    range = select, lower = lower,
+                                    overwrite_ab = overwrite_a_band,
+                                    abstol=abstol)
+        # crop off w and v
+        w = w[:m]
+        if not eigvals_only:
+            v = v[:, :m]
+    if info<0: raise ValueError,\
+    'illegal value in %-th argument of internal %s'%(-info, internal_name)
+    if info>0: raise LinAlgError,"eig algorithm did not converge"
+        
+    if eigvals_only:
+        return w
+    return w, v
+
 def eigvals(a,b=None,overwrite_a=0):
     """Return eigenvalues of square matrix."""
     return eig(a,b=b,left=0,right=0,overwrite_a=overwrite_a)
@@ -253,6 +373,13 @@
     """Return eigenvalues of hermitean or real symmetric matrix."""
     return eigh(a,lower=lower,eigvals_only=1,overwrite_a=overwrite_a)
 
+def eigvals_banded(a_band,lower=0,overwrite_a_band=0,
+                   select='a', select_range=None):
+    """Return eigenvalues of hermitean or real symmetric matrix."""
+    return eig_banded(a_band,lower=lower,eigvals_only=1,
+                      overwrite_a_band=overwrite_a_band, select=select,
+                      select_range=select_range)
+
 def lu_factor(a, overwrite_a=0):
     """Return raw LU decomposition of a matrix and pivots, for use in solving
     a system of linear equations.

Modified: trunk/Lib/linalg/generic_flapack.pyf
===================================================================
--- trunk/Lib/linalg/generic_flapack.pyf	2006-06-30 08:30:00 UTC (rev 2029)
+++ trunk/Lib/linalg/generic_flapack.pyf	2006-06-30 11:38:08 UTC (rev 2030)
@@ -856,8 +856,361 @@
 
 end subroutine <tchar>ggev
 
+! if anything is wrong with the following wrappers (until *gbtrs) 
+! blame Arnd Baecker and Johannes Loehnert and not Pearu
+subroutine <tchar=s,d>sbev(ab,compute_v,lower,n,ldab,kd,w,z,ldz,work,info) ! in :Band:dsbev.f
+  ! principally <s,d>sbevd does the same, and are recommended for use.
+  ! (see man dsbevd)
 
+  callstatement (*f2py_func)((compute_v?"V":"N"),(lower?"L":"U"),&n,&kd,ab,&ldab,w,z,&ldz,work,&info)
 
+  callprotoargument char*,char*,int*,int*,<type_in_c>*,int*,<type_in_c>*,<type_in_c>*,int*,<type_in_c>*,int*
+
+  ! Remark: if ab is fortran contigous on input 
+  !         and overwrite_ab=1  ab will be overwritten.
+  <type_in> dimension(ldab,*), intent(in,overwrite) :: ab
+
+  integer optional,intent(in):: compute_v = 1
+  check(compute_v==1||compute_v==0) compute_v
+  integer optional,intent(in),check(lower==0||lower==1) :: lower = 0   
+
+  integer optional,check(shape(ab,0)==ldab),depend(ab) :: ldab=shape(ab,0)
+  integer intent(hide),depend(ab) :: n=shape(ab,1)
+  integer intent(hide),depend(ab) :: kd=shape(ab,0)-1
+
+  <type_in> dimension(n),intent(out),depend(n) :: w
+  
+  ! For compute_v=1 z is used and contains the eigenvectors 
+  integer intent(hide),depend(n) :: ldz=(compute_v?n:1)
+  <type_in> dimension(ldz,ldz),intent(out),depend(ldz) :: z
+  
+  <type_in> dimension(MAX(1,3*n-1)),intent(hide),depend(n) :: work
+  integer intent(out)::info
+end subroutine <tchar=s,d>sbev
+
+
+
+subroutine <tchar=s,d>sbevd(ab,compute_v,lower,n,ldab,kd,w,z,ldz,work,lwork,iwork,liwork,info) ! in :Band:dsbevd.f
+
+  callstatement (*f2py_func)((compute_v?"V":"N"),(lower?"L":"U"),&n,&kd,ab,&ldab,w,z,&ldz,work,&lwork,iwork,&liwork,&info)
+
+  callprotoargument char*,char*,int*,int*,<type_in_c>*,int*,<type_in_c>*,<type_in_c>*,int*,<type_in_c>*,int*,int*,int*,int*
+
+  ! Remark: if ab is fortran contigous on input 
+  !         and overwrite_ab=1  ab will be overwritten.
+  <type_in> dimension(ldab,*), intent(in, overwrite) :: ab
+
+  integer optional,intent(in):: compute_v = 1
+  check( compute_v==1||compute_v==0) compute_v
+  integer optional,intent(in),check(lower==0||lower==1) :: lower = 0   
+
+  integer optional,check(shape(ab,0)==ldab),depend(ab) :: ldab=shape(ab,0)
+  integer intent(hide),depend(ab) :: n=shape(ab,1)
+  integer intent(hide),depend(ab) :: kd=shape(ab,0)-1
+
+  <type_in> dimension(n),intent(out),depend(n) :: w
+  <type_in> dimension(ldz,ldz),intent(out),depend(ldz) :: z
+  
+  ! For compute_v=1 z is used and contains the eigenvectors 
+  integer intent(hide),depend(n) :: ldz=(compute_v?n:1)
+  <type_in> dimension(ldz,ldz),depend(ldz) :: z
+  
+  integer intent(hide),depend(n) :: lwork=(compute_v?1+5*n+2*n*n:2*n)
+  <type_in> dimension(lwork),intent(hide),depend(lwork) :: work
+  integer intent(out)::info
+
+  integer optional,check(liwork>=(compute_v?3+5*n:1)),depend(n) :: liwork=(compute_v?3+5*n:1)
+  integer intent(hide),dimension(liwork),depend(liwork) :: iwork
+end subroutine <tchar=s,d>sbevd
+
+
+
+subroutine <tchar=s,d>sbevx(ab,ldab,compute_v,range,lower,n,kd,q,ldq,vl,vu,il,iu,abstol,w,z,m,mmax,ldz,work,iwork,ifail,info) ! in :Band:dsbevx.f
+
+  callstatement (*f2py_func)((compute_v?"V":"N"),(range>0?(range==1?"V":"I"):"A"),(lower?"L":"U"),&n,&kd,ab,&ldab,q,&ldq,&vl,&vu,&il,&iu,&abstol,&m,w,z,&ldz,work,iwork,ifail,&info)
+                             
+  callprotoargument char*,char*,char*,int*,int*,<type_in_c>*,int*,<type_in_c>*,int*,<type_in_c>*,<type_in_c>*,int*,int*,<type_in_c>*,int*,<type_in_c>*,<type_in_c>*,int*,<type_in_c>*, int*,int*,int*
+
+  integer optional,intent(in):: compute_v = 1
+  check(compute_v==1||compute_v==0) compute_v
+  integer optional,intent(in),check(lower==0||lower==1) :: lower = 0   
+
+  integer optional,check(shape(ab,0)==ldab),depend(ab) :: ldab=shape(ab,0)
+  integer intent(hide),depend(ab) :: n=shape(ab,1)
+  integer intent(hide),depend(ab) :: kd=shape(ab,0)-1
+
+  integer optional,intent(in):: range = 0
+  check(range==2||range==1||range==0) range
+
+
+  ! Remark: if ab is fortran contigous on input 
+  !         and overwrite_ab=1  ab will be overwritten.
+  <type_in> dimension(ldab,*),intent(in, overwrite) :: ab
+
+
+  ! FIXME: do we need to make q available for outside usage ???
+  !        If so: how to make this optional
+  !*  Q       (output) DOUBLE PRECISION array, dimension (LDQ, N)
+  !*          If JOBZ = 'V', the N-by-N orthogonal matrix used in the
+  !*                         reduction to tridiagonal form.
+  !*          If JOBZ = 'N', the array Q is not referenced.
+  integer intent(hide),depend(n) :: ldq=(compute_v?n:1)
+  <type_in> dimension(ldq,ldq),intent(hide),depend(ldq) :: q
+
+
+  <type_in> :: vl
+  <type_in> :: vu
+  integer,check((il>=1 && il<=n)),depend(n) :: il
+  integer,check((iu>=1 && iu<=n && iu>=il)),depend(n,il) :: iu
+
+  ! Remark, we don't use python indexing here, because
+  !  if someone uses ?sbevx directly,
+  !  he should expect Fortran style indexing.
+  !integer,check((il>=0 && il<n)),depend(n) :: il+1
+  !integer,check((iu>=0 && iu<n && iu>=il)),depend(n,il) :: iu+1
+
+  ! Remark:
+  ! Eigenvalues will be computed most accurately when ABSTOL is
+  ! set to twice the underflow threshold 2*DLAMCH('S'), not zero.
+  ! 
+  ! The easiest is to wrap DLAMCH (done below)
+  ! and let the user provide the value.
+  <type_in> optional,intent(in):: abstol=0.0
+
+  <type_in> dimension(n),intent(out),depend(n) :: w
+
+  <type_in> dimension(ldz,mmax),intent(out) :: z
+  integer intent(hide),depend(n) :: ldz=(compute_v?n:1)
+
+  ! We use the mmax parameter to fix the size of z
+  ! (only if eigenvalues are requested)
+  ! Otherwise we would allocate a (possibly) huge
+  ! region of memory for the eigenvectors, even
+  ! in cases where only a few are requested. 
+  ! If RANGE = 'V' (range=1) we a priori don't know the
+  ! number of eigenvalues in the interval in advance.
+  ! As default we use the maximum value
+  ! but the user should use an appropriate mmax.
+  integer intent(in),depend(n) :: mmax=(compute_v?(range==2?(iu-il+1):n):1)
+  integer intent(out) :: m
+
+  <type_in> dimension(7*n),intent(hide) :: work
+  integer dimension(5*n),intent(hide) :: iwork
+  integer dimension((compute_v?n:1)),intent(out) :: ifail
+  integer intent(out):: info
+end subroutine <tchar=s,d>sbevx
+
+
+subroutine <tchar=c,z>hbevd(ab,compute_v,lower,n,ldab,kd,w,z,ldz,work,lwork,rwork,lrwork,iwork,liwork,info) ! in :Band:zubevd.f
+
+  callstatement (*f2py_func)((compute_v?"V":"N"),(lower?"L":"U"),&n,&kd,ab,&ldab,w,z,&ldz,work,&lwork,rwork,&lrwork,iwork,&liwork,&info)
+
+  callprotoargument char*,char*,int*,int*,<type_in_c>*,int*,<rtype_in_c>*,<type_in_c>*,int*,<type_in_c>*,int*,<rtype_in_c>*,int*,int*,int*,int*
+
+  ! Remark: if ab is fortran contigous on input
+  !         and overwrite_ab=1  ab will be overwritten.
+  <type_in> dimension(ldab,*), intent(in, overwrite) :: ab
+
+  integer optional,intent(in):: compute_v = 1
+  check( compute_v==1||compute_v==0) compute_v
+  integer optional,intent(in),check(lower==0||lower==1) :: lower = 0
+
+  integer optional,check(shape(ab,0)==ldab),depend(ab) :: ldab=shape(ab,0)
+  integer intent(hide),depend(ab) :: n=shape(ab,1)
+  ! case n=0 is omitted in calculaton of lwork, lrwork, liwork
+  ! so we forbid it
+  check( n>0 ) n
+  integer intent(hide),depend(ab) :: kd=shape(ab,0)-1
+
+  <rtype_in> dimension(n),intent(out),depend(n) :: w
+
+  ! For compute_v=1 z is used and contains the eigenvectors
+  integer intent(hide),depend(n) :: ldz=(compute_v?n:1)
+  <type_in> dimension(ldz,ldz),intent(out),depend(ldz) :: z
+
+  integer intent(hide),depend(n) :: lwork=(compute_v?2*n*n:n)
+  <type_in> dimension(lwork),intent(hide),depend(lwork) :: work
+  integer intent(out)::info
+
+  integer optional, check(lrwork>=(compute_v?1+5*n+2*n*n:n)),depend(n) :: lrwork=(compute_v?1+5*n+2*n*n:n)
+
+  <rtype_in> intent(hide),dimension(lrwork),depend(lrwork) :: rwork
+
+  ! documentation says liwork >=2+5*n, but that crashes, +1 helps
+  integer optional, check(liwork>=(compute_v?3+5*n:1)),depend(n) :: liwork=(compute_v?3+5*n:1)
+  integer intent(hide),dimension(liwork),depend(liwork) :: iwork
+
+end subroutine <tchar=c,z>hbevd
+
+
+
+subroutine <tchar=c,z>hbevx(ab,ldab,compute_v,range,lower,n,kd,q,ldq,vl,vu,il,iu,abstol,w,z,m,mmax,ldz,work,rwork,iwork,ifail,info) ! in :Band:dsbevx.f
+
+  callstatement (*f2py_func)((compute_v?"V":"N"),(range>0?(range==1?"V":"I"):"A"),(lower?"L":"U"),&n,&kd,ab,&ldab,q,&ldq,&vl,&vu,&il,&iu,&abstol,&m,w,z,&ldz,work,rwork,iwork,ifail,&info)
+                             
+  callprotoargument
+  char*,char*,char*,int*,int*,<type_in_c>*,int*,<type_in_c>*,int*,<rtype_in_c>*,<rtype_in_c>*,int*,int*,<rtype_in_c>*,int*,<rtype_in_c>*,<type_in_c>*,int*,<type_in_c>*,<rtype_in_c>*,int*,int*,int*
+  
+  integer optional,intent(in):: compute_v = 1
+  check(compute_v==1||compute_v==0) compute_v
+  integer optional,intent(in),check(lower==0||lower==1) :: lower = 0   
+
+  integer optional,check(shape(ab,0)==ldab),depend(ab) :: ldab=shape(ab,0)
+  integer intent(hide),depend(ab) :: n=shape(ab,1)
+  integer intent(hide),depend(ab) :: kd=shape(ab,0)-1
+
+  integer optional,intent(in):: range = 0
+  check(range==2||range==1||range==0) range
+
+
+  ! Remark: if ab is fortran contigous on input 
+  !         and overwrite_ab=1  ab will be overwritten.
+  <type_in> dimension(ldab,*),intent(in, overwrite) :: ab
+
+
+  ! FIXME: do we need to make q available for outside usage ???
+  !        If so: how to make this optional
+  !*  Q       (output) DOUBLE PRECISION array, dimension (LDQ, N)
+  !*          If JOBZ = 'V', the N-by-N orthogonal matrix used in the
+  !*                         reduction to tridiagonal form.
+  !*          If JOBZ = 'N', the array Q is not referenced.
+  integer intent(hide),depend(n) :: ldq=(compute_v?n:1)
+  <type_in> dimension(ldq,ldq),intent(hide),depend(ldq) :: q
+
+
+  <rtype_in> :: vl
+  <rtype_in> :: vu
+  integer,check((il>=1 && il<=n)),depend(n) :: il
+  integer,check((iu>=1 && iu<=n && iu>=il)),depend(n,il) :: iu
+
+  ! Remark, we don't use python indexing here, because
+  !  if someone uses ?sbevx directly,
+  !  he should expect Fortran style indexing.
+  !integer,check((il>=0 && il<n)),depend(n) :: il+1
+  !integer,check((iu>=0 && iu<n && iu>=il)),depend(n,il) :: iu+1
+
+  ! Remark:
+  ! Eigenvalues will be computed most accurately when ABSTOL is
+  ! set to twice the underflow threshold 2*DLAMCH('S'), not zero.
+  ! 
+  ! The easiest is to wrap DLAMCH (done below)
+  ! and let the user provide the value.
+  <rtype_in> optional,intent(in):: abstol=0.0
+
+  <rtype_in> dimension(n),intent(out),depend(n) :: w
+
+  <type_in> dimension(ldz,mmax),intent(out) :: z
+  integer intent(hide),depend(n) :: ldz=(compute_v?n:1)
+
+  ! We use the mmax parameter to fix the size of z
+  ! (only if eigenvalues are requested)
+  ! Otherwise we would allocate a (possibly) huge
+  ! region of memory for the eigenvectors, even
+  ! in cases where only a few are requested. 
+  ! If RANGE = 'V' (range=1) we a priori don't know the
+  ! number of eigenvalues in the interval in advance.
+  ! As default we use the maximum value
+  ! but the user should use an appropriate mmax.
+  integer intent(in),depend(n) :: mmax=(compute_v?(range==2?(iu-il+1):n):1)
+  integer intent(out) :: m
+
+  <type_in> dimension(n),intent(hide) :: work
+  <rtype_in> dimension(7*n),intent(hide) :: rwork
+  integer dimension(5*n),intent(hide) :: iwork
+  integer dimension((compute_v?n:1)),intent(out) :: ifail
+  integer intent(out):: info
+end subroutine <tchar=c,z>hbevx
+
+
+! dlamch = dlamch(cmach)     
+!
+! determine double precision machine parameters
+!  CMACH   (input) CHARACTER*1
+!          Specifies the value to be returned by DLAMCH:
+!          = 'E' or 'e',   DLAMCH := eps
+!          = 'S' or 's ,   DLAMCH := sfmin
+!          = 'B' or 'b',   DLAMCH := base
+!          = 'P' or 'p',   DLAMCH := eps*base
+!          = 'N' or 'n',   DLAMCH := t
+!          = 'R' or 'r',   DLAMCH := rnd
+!          = 'M' or 'm',   DLAMCH := emin
+!          = 'U' or 'u',   DLAMCH := rmin
+!          = 'L' or 'l',   DLAMCH := emax
+!          = 'O' or 'o',   DLAMCH := rmax
+!
+!          where
+!
+!          eps   = relative machine precision
+!          sfmin = safe minimum, such that 1/sfmin does not overflow
+!          base  = base of the machine
+!          prec  = eps*base
+!          t     = number of (base) digits in the mantissa
+!          rnd   = 1.0 when rounding occurs in addition, 0.0 otherwise
+!          emin  = minimum exponent before (gradual) underflow
+!          rmin  = underflow threshold - base**(emin-1)
+!          emax  = largest exponent before overflow
+!          rmax  = overflow threshold  - (base**emax)*(1-eps)
+function <tchar=s,d>lamch(cmach)
+    character :: cmach
+    <type_in> intent(out):: dlamch
+end function <tchar=s,d>lamch
+
+
+
+! lu,ipiv,info = dgbtrf(ab,kl,ku,[m,n,ldab,overwrite_ab])
+! Compute  an  LU factorization of a real m-by-n band matrix 
+subroutine <tchar=s,d,c,z>gbtrf(m,n,ab,kl,ku,ldab,ipiv,info) ! in :Band:dgbtrf.f
+  ! threadsafe  ! FIXME: should this be added ?
+
+  callstatement {int i;(*f2py_func)(&m,&n,&kl,&ku,ab,&ldab,ipiv,&info); for(i=0,n=MIN(m,n);i<n;--ipiv[i++]);}
+                             
+  callprotoargument int*,int*,int*,int*,<type_in_c>*,int*,int*,int*
+  
+    ! let the default be a square matrix:
+    integer optional,depend(ab) :: m=shape(ab,1)
+    integer optional,depend(ab) :: n=shape(ab,1)
+    integer :: kl
+    integer :: ku
+
+    <type_in> dimension(ldab,*),intent(in,out,copy,out=lu) :: ab
+    integer optional,check(shape(ab,0)==ldab),depend(ab) :: ldab=shape(ab,0)
+    integer dimension(MIN(m,n)),depend(m,n),intent(out) :: ipiv
+    integer intent(out):: info
+end subroutine <tchar=s,d,c,z>gbtrf
+
+
+
+subroutine <tchar=s,d,c,z>gbtrs(ab,kl,ku,b,ipiv,trans,n,nrhs,ldab,ldb,info) ! in :Band:dgbtrs.f
+! x,info = dgbtrs(ab,kl,ku,b,ipiv,[trans,n,ldab,ldb,overwrite_b])
+! solve a system of linear equations A * X = B or A' * X = B
+! with a general band matrix A using the  LU  factorization  
+! computed by DGBTRF 
+!
+! TRANS   Specifies the form of the system of equations.  
+!  0  = 'N':  A * X =B  (No transpose)
+!  1  = 'T':  A'* X = B  (Transpose)
+!  2  = 'C':  A'* X = B  (Conjugate transpose = Transpose)
+
+callstatement {int i;for(i=0;i<n;++ipiv[i++]);(*f2py_func)((trans>0?(trans==1?"T":"C"):"N"),&n,&kl,&ku,&nrhs,ab,&ldab,ipiv,b,&ldb,&info);for(i=0;i<n;--ipiv[i++]);}
+lprotoargument char*,int*,int *,int*,int*,<type_in_c>*,int*,int*,<type_in_c>*,int*,int*
+    !character optional:: trans='N'
+    integer optional:: trans=0
+    integer optional,depend(ab) :: n=shape(ab,1)
+    integer :: kl
+    integer :: ku
+    integer intent(hide),depend(b):: nrhs=shape(b,1)
+
+    <type_in> dimension(ldab,*),intent(in) :: ab
+    integer optional,check(shape(ab,0)==ldab),depend(ab) :: ldab=shape(ab,0)
+
+    integer dimension(n),intent(in) :: ipiv
+    <type_in> dimension(ldb,*),intent(in,out,copy,out=x) :: b
+    integer optional,check(shape(b,0)==ldb),depend(b) :: ldb=shape(b,0)
+    integer optional,check(shape(b,0)==ldb),depend(b) :: ldb=shape(b,0)
+    integer intent(out):: info
+end subroutine <tchar=s,d,c,z>gbtrs
+
+
 end interface
 
 end python module generic_flapack

Modified: trunk/Lib/linalg/tests/test_decomp.py
===================================================================
--- trunk/Lib/linalg/tests/test_decomp.py	2006-06-30 08:30:00 UTC (rev 2029)
+++ trunk/Lib/linalg/tests/test_decomp.py	2006-06-30 11:38:08 UTC (rev 2030)
@@ -20,6 +20,10 @@
 set_package_path()
 from linalg import eig,eigvals,lu,svd,svdvals,cholesky,qr,schur,rsf2csf
 from linalg import lu_solve,lu_factor,solve,diagsvd,hessenberg
+from linalg import eig_banded,eigvals_banded
+from linalg.flapack import dgbtrf, dgbtrs, zgbtrf, zgbtrs
+from linalg.flapack import dsbev, dsbevd, dsbevx, zhbevd, zhbevx
+
 restore_path()
 
 from numpy import *
@@ -105,6 +109,293 @@
             assert_array_almost_equal(dot(conjugate(transpose(a)),vl[:,i]),
                                       conjugate(w[i])*vl[:,i])
 
+
+
+class test_eig_banded(ScipyTestCase):
+
+    def __init__(self, *args):
+        ScipyTestCase.__init__(self, *args)
+
+        self.create_bandmat()
+
+    def create_bandmat(self):
+        """Create the full matrix `self.fullmat` and 
+           the corresponding band matrix `self.bandmat`."""
+        N  = 10
+        self.KL = 2   # number of subdiagonals (below the diagonal)
+        self.KU = 2   # number of superdiagonals (above the diagonal)
+
+        # symmetric band matrix
+        self.sym_mat = ( diag(1.0*ones(N))
+                     +  diag(-1.0*ones(N-1), -1) + diag(-1.0*ones(N-1), 1) 
+                     + diag(-2.0*ones(N-2), -2) + diag(-2.0*ones(N-2), 2) )
+
+        # hermitian band matrix
+        self.herm_mat = ( diag(-1.0*ones(N))
+                     + 1j*diag(1.0*ones(N-1), -1) - 1j*diag(1.0*ones(N-1), 1)
+                     + diag(-2.0*ones(N-2), -2) + diag(-2.0*ones(N-2), 2) )
+
+        # general real band matrix
+        self.real_mat = ( diag(1.0*ones(N))
+                     +  diag(-1.0*ones(N-1), -1) + diag(-3.0*ones(N-1), 1) 
+                     + diag(2.0*ones(N-2), -2) + diag(-2.0*ones(N-2), 2) )
+
+        # general complex band matrix
+        self.comp_mat = ( 1j*diag(1.0*ones(N))
+                     +  diag(-1.0*ones(N-1), -1) + 1j*diag(-3.0*ones(N-1), 1) 
+                     + diag(2.0*ones(N-2), -2) + diag(-2.0*ones(N-2), 2) )
+
+
+        # Eigenvalues and -vectors from linalg.eig
+        ew, ev = linalg.eig(self.sym_mat)
+        ew = ew.real
+        args = argsort(ew)
+        self.w_sym_lin = ew[args]
+        self.evec_sym_lin = ev[:,args]
+
+        ew, ev = linalg.eig(self.herm_mat)
+        ew = ew.real
+        args = argsort(ew)
+        self.w_herm_lin = ew[args]
+        self.evec_herm_lin = ev[:,args]
+
+
+        # Extract upper bands from symmetric and hermitian band matrices
+        # (for use in dsbevd, dsbevx, zhbevd, zhbevx
+        #  and their single precision versions) 
+        LDAB = self.KU + 1
+        self.bandmat_sym  = zeros((LDAB, N), dtype=float)
+        self.bandmat_herm = zeros((LDAB, N), dtype=complex)
+        for i in xrange(LDAB):
+            self.bandmat_sym[LDAB-i-1,i:N]  = diag(self.sym_mat, i)
+            self.bandmat_herm[LDAB-i-1,i:N] = diag(self.herm_mat, i)
+
+
+        # Extract bands from general real and complex band matrix
+        # (for use in dgbtrf, dgbtrs and their single precision versions)
+        LDAB = 2*self.KL + self.KU + 1
+        self.bandmat_real = zeros((LDAB, N), dtype=float)
+        self.bandmat_real[2*self.KL,:] = diag(self.real_mat)     # diagonal
+        for i in xrange(self.KL):
+            # superdiagonals
+            self.bandmat_real[2*self.KL-1-i,i+1:N]   = diag(self.real_mat, i+1)
+            # subdiagonals
+            self.bandmat_real[2*self.KL+1+i,0:N-1-i] = diag(self.real_mat,-i-1)
+
+        self.bandmat_comp = zeros((LDAB, N), dtype=complex)
+        self.bandmat_comp[2*self.KL,:] = diag(self.comp_mat)     # diagonal
+        for i in xrange(self.KL):
+            # superdiagonals
+            self.bandmat_comp[2*self.KL-1-i,i+1:N]   = diag(self.comp_mat, i+1)
+            # subdiagonals
+            self.bandmat_comp[2*self.KL+1+i,0:N-1-i] = diag(self.comp_mat,-i-1)
+
+        # absolute value for linear equation system A*x = b
+        self.b = 1.0*arange(N)
+        self.bc = self.b *(1 + 1j) 
+        
+
+    #####################################################################
+
+        
+    def check_dsbev(self):
+        """Compare dsbev eigenvalues and eigenvectors with
+           the result of linalg.eig."""
+        w, evec, info  = dsbev(self.bandmat_sym, compute_v=1)
+        evec_ = evec[:,argsort(w)]
+        assert_array_almost_equal(sort(w), self.w_sym_lin)
+        assert_array_almost_equal(abs(evec_), abs(self.evec_sym_lin))
+
+
+    
+    def check_dsbevd(self):
+        """Compare dsbevd eigenvalues and eigenvectors with
+           the result of linalg.eig."""
+        w, evec, info = dsbevd(self.bandmat_sym, compute_v=1)
+        evec_ = evec[:,argsort(w)]
+        assert_array_almost_equal(sort(w), self.w_sym_lin)
+        assert_array_almost_equal(abs(evec_), abs(self.evec_sym_lin))
+
+
+
+    def check_dsbevx(self):
+        """Compare dsbevx eigenvalues and eigenvectors
+           with the result of linalg.eig."""
+        N,N = shape(self.sym_mat)
+        ## Achtung: Argumente 0.0,0.0,range?
+        w, evec, num, ifail, info = dsbevx(self.bandmat_sym, 0.0, 0.0, 1, N,
+                                       compute_v=1, range=2)
+        evec_ = evec[:,argsort(w)]
+        assert_array_almost_equal(sort(w), self.w_sym_lin)
+        assert_array_almost_equal(abs(evec_), abs(self.evec_sym_lin))
+
+
+    def check_zhbevd(self):
+        """Compare zhbevd eigenvalues and eigenvectors
+           with the result of linalg.eig."""
+        w, evec, info = zhbevd(self.bandmat_herm, compute_v=1)
+        evec_ = evec[:,argsort(w)]
+        assert_array_almost_equal(sort(w), self.w_herm_lin)
+        assert_array_almost_equal(abs(evec_), abs(self.evec_herm_lin))
+
+
+
+    def check_zhbevx(self):
+        """Compare zhbevx eigenvalues and eigenvectors
+           with the result of linalg.eig."""
+        N,N = shape(self.herm_mat)
+        ## Achtung: Argumente 0.0,0.0,range?
+        w, evec, num, ifail, info = zhbevx(self.bandmat_herm, 0.0, 0.0, 1, N,
+                                       compute_v=1, range=2)
+        evec_ = evec[:,argsort(w)]
+        assert_array_almost_equal(sort(w), self.w_herm_lin)
+        assert_array_almost_equal(abs(evec_), abs(self.evec_herm_lin))
+
+
+
+    def check_eigvals_banded(self):
+        """Compare eigenvalues of eigvals_banded with those of linalg.eig."""
+        w_sym = eigvals_banded(self.bandmat_sym)
+        w_sym = w_sym.real
+        assert_array_almost_equal(sort(w_sym), self.w_sym_lin)
+
+        w_herm = eigvals_banded(self.bandmat_herm)
+        w_herm = w_herm.real
+        assert_array_almost_equal(sort(w_herm), self.w_herm_lin)
+
+        # extracting eigenvalues with respect to an index range
+        ind1 = 2   
+        ind2 = 6
+        w_sym_ind = eigvals_banded(self.bandmat_sym,
+                                    select='i', select_range=(ind1, ind2) )
+        assert_array_almost_equal(sort(w_sym_ind),
+                                  self.w_sym_lin[ind1:ind2+1])
+        w_herm_ind = eigvals_banded(self.bandmat_herm,
+                                    select='i', select_range=(ind1, ind2) )
+        assert_array_almost_equal(sort(w_herm_ind),
+                                  self.w_herm_lin[ind1:ind2+1])
+
+        # extracting eigenvalues with respect to a value range
+        v_lower = self.w_sym_lin[ind1] - 1.0e-5
+        v_upper = self.w_sym_lin[ind2] + 1.0e-5
+        w_sym_val = eigvals_banded(self.bandmat_sym,
+                                select='v', select_range=(v_lower, v_upper) )
+        assert_array_almost_equal(sort(w_sym_val),
+                                  self.w_sym_lin[ind1:ind2+1])
+
+        v_lower = self.w_herm_lin[ind1] - 1.0e-5
+        v_upper = self.w_herm_lin[ind2] + 1.0e-5
+        w_herm_val = eigvals_banded(self.bandmat_herm,
+                                select='v', select_range=(v_lower, v_upper) )
+        assert_array_almost_equal(sort(w_herm_val),
+                                  self.w_herm_lin[ind1:ind2+1])
+
+
+
+    def check_eig_banded(self):
+        """Compare eigenvalues and eigenvectors of eig_banded
+           with those of linalg.eig. """
+        w_sym, evec_sym = eig_banded(self.bandmat_sym)
+        evec_sym_ = evec_sym[:,argsort(w_sym.real)]
+        assert_array_almost_equal(sort(w_sym), self.w_sym_lin)
+        assert_array_almost_equal(abs(evec_sym_), abs(self.evec_sym_lin))
+
+        w_herm, evec_herm = eig_banded(self.bandmat_herm)
+        evec_herm_ = evec_herm[:,argsort(w_herm.real)]
+        assert_array_almost_equal(sort(w_herm), self.w_herm_lin)
+        assert_array_almost_equal(abs(evec_herm_), abs(self.evec_herm_lin))
+        
+        # extracting eigenvalues with respect to an index range
+        ind1 = 2   
+        ind2 = 6
+        w_sym_ind, evec_sym_ind = eig_banded(self.bandmat_sym,
+                                    select='i', select_range=(ind1, ind2) )
+        assert_array_almost_equal(sort(w_sym_ind),
+                                  self.w_sym_lin[ind1:ind2+1])
+        assert_array_almost_equal(abs(evec_sym_ind),
+                                  abs(self.evec_sym_lin[:,ind1:ind2+1]) )
+
+        w_herm_ind, evec_herm_ind = eig_banded(self.bandmat_herm,
+                                    select='i', select_range=(ind1, ind2) )
+        assert_array_almost_equal(sort(w_herm_ind),
+                                  self.w_herm_lin[ind1:ind2+1])
+        assert_array_almost_equal(abs(evec_herm_ind),
+                                  abs(self.evec_herm_lin[:,ind1:ind2+1]) )
+
+        # extracting eigenvalues with respect to a value range
+        v_lower = self.w_sym_lin[ind1] - 1.0e-5
+        v_upper = self.w_sym_lin[ind2] + 1.0e-5
+        w_sym_val, evec_sym_val = eig_banded(self.bandmat_sym,
+                                select='v', select_range=(v_lower, v_upper) )
+        assert_array_almost_equal(sort(w_sym_val),
+                                  self.w_sym_lin[ind1:ind2+1])
+        assert_array_almost_equal(abs(evec_sym_val),
+                                  abs(self.evec_sym_lin[:,ind1:ind2+1]) )
+
+        v_lower = self.w_herm_lin[ind1] - 1.0e-5
+        v_upper = self.w_herm_lin[ind2] + 1.0e-5
+        w_herm_val, evec_herm_val = eig_banded(self.bandmat_herm,
+                                select='v', select_range=(v_lower, v_upper) )
+        assert_array_almost_equal(sort(w_herm_val),
+                                  self.w_herm_lin[ind1:ind2+1])
+        assert_array_almost_equal(abs(evec_herm_val),
+                                  abs(self.evec_herm_lin[:,ind1:ind2+1]) )
+
+
+    def check_dgbtrf(self):
+        """Compare dgbtrf  LU factorisation with the LU factorisation result
+           of linalg.lu."""
+        M,N = shape(self.real_mat)        
+        lu_symm_band, ipiv, info = dgbtrf(self.bandmat_real, self.KL, self.KU)
+
+        # extract matrix u from lu_symm_band
+        u = diag(lu_symm_band[2*self.KL,:])
+        for i in xrange(self.KL + self.KU):
+            u += diag(lu_symm_band[2*self.KL-1-i,i+1:N], i+1)
+
+        p_lin, l_lin, u_lin = lu(self.real_mat, permute_l=0)
+        assert_array_almost_equal(u, u_lin)
+
+
+    def check_zgbtrf(self):
+        """Compare zgbtrf  LU factorisation with the LU factorisation result
+           of linalg.lu."""
+        M,N = shape(self.comp_mat)        
+        lu_symm_band, ipiv, info = zgbtrf(self.bandmat_comp, self.KL, self.KU)
+
+        # extract matrix u from lu_symm_band
+        u = diag(lu_symm_band[2*self.KL,:])
+        for i in xrange(self.KL + self.KU):
+            u += diag(lu_symm_band[2*self.KL-1-i,i+1:N], i+1)
+
+        p_lin, l_lin, u_lin =lu(self.comp_mat, permute_l=0)
+        assert_array_almost_equal(u, u_lin)
+
+
+
+    def check_dgbtrs(self):
+        """Compare dgbtrs  solutions for linear equation system  A*x = b
+           with solutions of linalg.solve."""
+        
+        lu_symm_band, ipiv, info = dgbtrf(self.bandmat_real, self.KL, self.KU)
+        y, info = dgbtrs(lu_symm_band, self.KL, self.KU, self.b, ipiv)
+
+        y_lin = linalg.solve(self.real_mat, self.b)
+        assert_array_almost_equal(y, y_lin)
+
+    def check_zgbtrs(self):
+        """Compare zgbtrs  solutions for linear equation system  A*x = b
+           with solutions of linalg.solve."""
+        
+        lu_symm_band, ipiv, info = zgbtrf(self.bandmat_comp, self.KL, self.KU)
+        y, info = zgbtrs(lu_symm_band, self.KL, self.KU, self.bc, ipiv)
+
+        y_lin = linalg.solve(self.comp_mat, self.bc)
+        assert_array_almost_equal(y, y_lin)
+
+
+
+
 class test_lu(ScipyTestCase):
 
     def check_simple(self):



More information about the Scipy-svn mailing list