[Scipy-svn] r2030 - in trunk/Lib/linalg: . tests
scipy-svn at scipy.org
scipy-svn at scipy.org
Fri Jun 30 07:38:14 EDT 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