!%f90 -*- f90 -*-
! Signatures for f2py wrappers of FORTRAN LAPACK functions.
!
! Author: Pearu Peterson
! Created: Jan-Feb 2002
! $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

python module flapack
interface

   subroutine <tchar=s,d>pbtrf(lower,n,kd,ab,ldab,info)
   
     ! Compute Cholesky decomposition of banded symmetric positive definite
     ! matrix:
     ! A = U^T * U, C = U if lower = 0
     ! A = L * L^T, C = L if lower = 1
     ! C is triangular matrix of the corresponding Cholesky decomposition.

     callstatement (*f2py_func)((lower?"L":"U"),&n,&kd,ab,&ldab,&info); 
     callprotoargument char*,int*,int*,<type_in_c>*,int*,int*

     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),check(lower==0||lower==1) :: lower = 0   

     <type_in> dimension(ldab,n),intent(in,out,copy,out=c) :: ab
     integer intent(out) :: info
     
   end subroutine <tchar=s,d>pbtrf

   subroutine <tchar=c,z>pbtrf(lower,n,kd,ab,ldab,info)
   

     ! Compute Cholesky decomposition of banded symmetric positive definite
     ! matrix:
     ! A = U^H * U, C = U if lower = 0
     ! A = L * L^H, C = L if lower = 1
     ! C is triangular matrix of the corresponding Cholesky decomposition.

     callstatement (*f2py_func)((lower?"L":"U"),&n,&kd,ab,&ldab,&info); 
     callprotoargument char*,int*,int*,<type_in_c>*,int*,int*

     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),check(lower==0||lower==1) :: lower = 0   

     <type_in> dimension(ldab,n),intent(in,out,copy,out=c) :: ab
     integer intent(out) :: info
     
   end subroutine <tchar=c,z>pbtrf

   subroutine <tchar=s,d>pbsv(lower,n,kd,nrhs,ab,ldab,b,ldb,info)
   
     !
     ! Computes the solution to a real system of linear equations
     ! A * X = B,
     !  where A is an N-by-N symmetric positive definite band matrix and X
     !  and B are N-by-NRHS matrices.
     !
     !  The Cholesky decomposition is used to factor A as
     !     A = U**T * U,  if lower=1, or
     !     A = L * L**T,  if lower=0
     !  where U is an upper triangular band matrix, and L is a lower
     !  triangular band matrix, with the same number of superdiagonals or
     !  subdiagonals as A.  The factored form of A is then used to solve the
     !  system of equations A * X = B.

     callstatement (*f2py_func)((lower?"L":"U"),&n,&kd,&nrhs,ab,&ldab,b,&ldb,&info); 
     callprotoargument char*,int*,int*,int*,<type_in_c>*,int*,<type_in_c>*,int*,int*

     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 intent(hide),depend(b) :: ldb=shape(b,1)
     integer intent(hide),depend(b) :: nrhs=shape(b,0)
     integer optional,intent(in),check(lower==0||lower==1) :: lower = 0   

     <type_in> dimension(nrhs,ldb),intent(in,out,copy,out=x) :: b
     <type_in> dimension(ldab,n),intent(in,out,copy,out=c) :: ab
     integer intent(out) :: info
     
   end subroutine <tchar=s,d>pbsv

   subroutine <tchar=c,z>pbsv(lower,n,kd,nrhs,ab,ldab,b,ldb,info)
   
     !
     ! Computes the solution to a real system of linear equations
     ! A * X = B,
     !  where A is an N-by-N Hermitian positive definite band matrix and X
     !  and B are N-by-NRHS matrices.
     !
     !  The Cholesky decomposition is used to factor A as
     !     A = U**H * U,  if lower=1, or
     !     A = L * L**H,  if lower=0
     !  where U is an upper triangular band matrix, and L is a lower
     !  triangular band matrix, with the same number of superdiagonals or
     !  subdiagonals as A.  The factored form of A is then used to solve the
     !  system of equations A * X = B.

     callstatement (*f2py_func)((lower?"L":"U"),&n,&kd,&nrhs,ab,&ldab,b,&ldb,&info); 
     callprotoargument char*,int*,int*,int*,<type_in_c>*,int*,<type_in_c>*,int*,int*

     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 intent(hide),depend(b) :: ldb=shape(b,1)
     integer intent(hide),depend(b) :: nrhs=shape(b,0)
     integer optional,intent(in),check(lower==0||lower==1) :: lower = 0   

     <type_in> dimension(nrhs,ldb),intent(in,out,copy,out=x) :: b
     <type_in> dimension(ldab,n),intent(in,out,copy,out=c) :: ab
     integer intent(out) :: info
     
   end subroutine <tchar=c,z>pbsv

   subroutine <tchar=s,d,c,z>gebal(scale,permute,n,a,m,lo,hi,pivscale,info)
   !
   ! ba,lo,hi,pivscale,info = gebal(a,scale=0,permute=0,overwrite_a=0)
   ! Balance general matrix a.
   ! hi,lo are such that ba[i][j]==0 if i>j and j=0...lo-1 or i=hi+1..n-1
   ! pivscale([0:lo], [lo:hi+1], [hi:n+1]) = (p1,d,p2) where (p1,p2)[j] is
   ! the index of the row and column interchanged with row and column j. 
   ! d[j] is the scaling factor applied to row and column j.
   ! The order in which the interchanges are made is n-1 to hi+1, then 0 to lo-1.
   !
   ! P * A * P = [[T1,X,Y],[0,B,Z],[0,0,T2]]
   ! BA = [[T1,X*D,Y],[0,inv(D)*B*D,ind(D)*Z],[0,0,T2]]
   ! where D = diag(d), T1,T2 are upper triangular matrices.
   ! lo,hi mark the starting and ending columns of submatrix B.

     callstatement { (*f2py_func)((permute?(scale?"B":"P"):(scale?"S":"N")),&n,a,&m,&lo,&hi,pivscale,&info); hi--; lo--; }
     callprotoargument char*,int*,<type_in_c>*,int*,int*,int*,<rtype_in_c>*,int*
     integer intent(in),optional :: permute = 0
     integer intent(in),optional :: scale = 0
     integer intent(hide),depend(a,n) :: m = shape(a,0)
     integer intent(hide),depend(a) :: n = shape(a,1)
     check(m>=n) m
     integer intent(out) :: hi,lo
     <rtype_in> dimension(n),intent(out),depend(n) :: pivscale
     <type_in> dimension(m,n),intent(in,out,copy,out=ba) :: a
     integer intent(out) :: info

   end subroutine <tchar=s,d,c,z>gebal

   subroutine <tchar=s,d,c,z>gehrd(n,lo,hi,a,tau,work,lwork,info)
   !
   ! hq,tau,info = gehrd(a,lo=0,hi=n-1,lwork=n,overwrite_a=0)
   ! Reduce general matrix A to upper Hessenberg form H by unitary similarity
   ! transform Q^H * A * Q = H
   !
   ! Q = H(lo) * H(lo+1) * ... * H(hi-1)
   ! H(i) = I - tau * v * v^H
   ! v[0:i+1] = 0, v[i+1]=1, v[hi+1:n] = 0
   ! v[i+2:hi+1] is stored in hq[i+2:hi+i,i]
   ! tau is tau[i]
   !
   ! hq for n=7,lo=1,hi=5:
   ! [a a h h h h a
   !    a h h h h a
   !    h h h h h h
   !    v2h h h h h
   !    v2v3h h h h
   !    v2v3v4h h h
   !              a]
   !
     callstatement { hi++; lo++; (*f2py_func)(&n,&lo,&hi,a,&n,tau,work,&lwork,&info); }
     callprotoargument int*,int*,int*,<type_in_c>*,int*,<type_in_c>*,<type_in_c>*,int*,int*
     integer intent(hide),depend(a) :: n = shape(a,0)
     <type_in> dimension(n,n),intent(in,out,copy,out=ht),check(shape(a,0)==shape(a,1)) :: a
     integer intent(in),optional :: lo = 0
     integer intent(in),optional,depend(n) :: hi = n-1
     <type_in> dimension(n-1),intent(out),depend(n) :: tau
     <type_in> dimension(lwork),intent(cahce,hide),depend(lwork) :: work
     integer intent(in),optional,depend(n),check(lwork>=MAX(n,1)) :: lwork = MAX(n,1)
     integer intent(out) :: info

   end subroutine <tchar=s,d,c,z>gehrd

   subroutine <tchar=s,d,c,z>gbsv(n,kl,ku,nrhs,ab,piv,b,info)
   ! 
   ! lub,piv,x,info = gbsv(kl,ku,ab,b,overwrite_ab=0,overwrite_b=0)
   ! Solve A * X = B
   ! A = P * L * U
   ! A is a band matrix of order n with kl subdiagonals and ku superdiagonals
   ! starting at kl-th row.
   ! X, B are n-by-nrhs matrices
   !
     callstatement {int i=2*kl+ku+1;(*f2py_func)(&n,&kl,&ku,&nrhs,ab,&i,piv,b,&n,&info);for(i=0;i<n;--piv[i++]);}
     callprotoargument int*,int*,int*,int*,<type_in_c>*,int*,int*,<type_in_c>*,int*,int*
     integer depend(ab),intent(hide):: n = shape(ab,1)
     integer intent(in) :: kl
     integer intent(in) :: ku
     integer depend(b),intent(hide) :: nrhs = shape(b,1)
     <type_in> dimension(2*kl+ku+1,n),depend(kl,ku), check(2*kl+ku+1==shape(ab,0)) :: ab
     integer dimension(n),depend(n),intent(out) :: piv
     <type_in> dimension(n,nrhs),depend(n),check(shape(ab,1)==shape(b,0)) :: b
     integer intent(out) :: info
     intent(in,out,copy,out=x) b
     intent(in,out,copy,out=lub) ab
   end subroutine <tchar=s,d,c,z>gbsv

   subroutine <tchar=s,d,c,z>gesv(n,nrhs,a,piv,b,info)

   ! lu,piv,x,info = gesv(a,b,overwrite_a=0,overwrite_b=0)
   ! Solve A * X = B.
   ! A = P * L * U
   ! U is upper diagonal triangular, L is unit lower triangular,
   ! piv pivots columns.

     callstatement {int i;(*f2py_func)(&n,&nrhs,a,&n,piv,b,&n,&info);for(i=0;i<n;--piv[i++]);}
     callprotoargument int*,int*,<type_in_c>*,int*,int*,<type_in_c>*,int*,int*

     integer depend(a),intent(hide):: n = shape(a,0)
     integer depend(b),intent(hide):: nrhs = shape(b,1)
     <type_in> dimension(n,n),check(shape(a,0)==shape(a,1)) :: a
     integer dimension(n),depend(n),intent(out) :: piv
     <type_in> dimension(n,nrhs),check(shape(a,0)==shape(b,0)),depend(n) :: b
     integer intent(out)::info
     intent(in,out,copy,out=x) b
     intent(in,out,copy,out=lu) a

   end subroutine <tchar=s,d,c,z>gesv

   subroutine <tchar=s,d,c,z>getrf(m,n,a,piv,info)

   ! lu,piv,info = getrf(a,overwrite_a=0)
   ! Compute an LU factorization of a  general  M-by-N  matrix  A.
   ! A = P * L * U
     threadsafe
     callstatement {int i;(*f2py_func)(&m,&n,a,&m,piv,&info);for(i=0,n=MIN(m,n);i<n;--piv[i++]);}
     callprotoargument int*,int*,<type_in_c>*,int*,int*,int*

     integer depend(a),intent(hide):: m = shape(a,0)
     integer depend(a),intent(hide):: n = shape(a,1)
     <type_in> dimension(m,n),intent(in,out,copy,out=lu) :: a
     integer dimension(MIN(m,n)),depend(m,n),intent(out) :: piv
     integer intent(out):: info

   end subroutine <tchar=s,d,c,z>getrf

   subroutine <tchar=s,d,c,z>getrs(n,nrhs,lu,piv,b,info,trans)

   ! x,info = getrs(lu,piv,b,trans=0,overwrite_b=0)
   ! Solve  A  * X = B if trans=0
   ! Solve A^T * X = B if trans=1
   ! Solve A^H * X = B if trans=2
   ! A = P * L * U
     threadsafe
     callstatement {int i;for(i=0;i<n;++piv[i++]);(*f2py_func)((trans?(trans==2?"C":"T"):"N"),&n,&nrhs,lu,&n,piv,b,&n,&info);for(i=0;i<n;--piv[i++]);}
     callprotoargument char*,int*,int*,<type_in_c>*,int*,int*,<type_in_c>*,int*,int*

     integer optional,intent(in),check(trans>=0 && trans <=2) :: trans = 0

     integer depend(lu),intent(hide):: n = shape(lu,0)
     integer depend(b),intent(hide):: nrhs = shape(b,1)
     <type_in> dimension(n,n),intent(in) :: lu
     check(shape(lu,0)==shape(lu,1)) :: lu
     integer dimension(n),intent(in),depend(n) :: piv
     <type_in> dimension(n,nrhs),intent(in,out,copy,out=x),depend(n),check(shape(lu,0)==shape(b,0)) :: b
     integer intent(out):: info
   end subroutine <tchar=s,d,c,z>getrs

   subroutine <tchar=s,d,c,z>getri(n,lu,piv,work,lwork,info)

   ! inv_a,info = getri(lu,piv,lwork=3*n,overwrite_lu=0)
   ! Find A inverse A^-1.
   ! A = P * L * U

     callstatement {int i;for(i=0;i<n;++piv[i++]);(*f2py_func)(&n,lu,&n,piv,work,&lwork,&info);for(i=0;i<n;--piv[i++]);}
     callprotoargument int*,<type_in_c>*,int*,int*,<type_in_c>*,int*,int*

     integer depend(lu),intent(hide):: n = shape(lu,0)
     <type_in> dimension(n,n),intent(in,out,copy,out=inv_a) :: lu
     check(shape(lu,0)==shape(lu,1)) :: lu
     integer dimension(n),intent(in),depend(n) :: piv
     integer intent(out):: info
     integer optional,intent(in),depend(n),check(lwork>=n) :: lwork=3*n
     <type_in> dimension(lwork),intent(hide,cache),depend(lwork) :: work

   end subroutine <tchar=s,d,c,z>getri

   subroutine <tchar=s,d>gesdd(m,n,minmn,du,dvt,a,compute_uv,u,s,vt,work,lwork,iwork,info)
   
   ! u,s,vt,info = gesdd(a,compute_uv=1,lwork=..,overwrite_a=0)
   ! Compute the singular value decomposition (SVD):
   !   A = U * SIGMA * transpose(V)
   ! A - M x N matrix
   ! U - M x M matrix
   ! SIGMA - M x N zero matrix with a main diagonal filled with min(M,N) 
   !               singular values  
   ! transpose(V) - N x N matrix

   callstatement (*f2py_func)((compute_uv?"A":"N"),&m,&n,a,&m,s,u,&du,vt,&dvt,work,&lwork,iwork,&info)
   callprotoargument char*,int*,int*,<type_in_c>*,int*,<type_in_c>*,<type_in_c>*,int*,<type_in_c>*,int*,<type_in_c>*,int*,int*,int*

   integer intent(in),optional,check(compute_uv==0||compute_uv==1):: compute_uv = 1
   integer intent(hide),depend(a):: m = shape(a,0)
   integer intent(hide),depend(a):: n = shape(a,1)
   integer intent(hide),depend(m,n):: minmn = MIN(m,n)
   integer intent(hide),depend(compute_uv,minmn) :: du = (compute_uv?m:1)
   integer intent(hide),depend(compute_uv,n) :: dvt = (compute_uv?n:1)
   <type_in> dimension(m,n),intent(in,copy) :: a
   <type_in> dimension(minmn),intent(out),depend(minmn) :: s
   <type_in> dimension(du,du),intent(out),depend(du) :: u
   <type_in> dimension(dvt,dvt),intent(out),depend(dvt) :: vt
   <type_in> dimension(lwork),intent(hide,cache),depend(lwork) :: work
   integer optional,intent(in),depend(minmn,compute_uv) &
        :: lwork = (compute_uv?4*minmn*minmn+MAX(m,n)+9*minmn:MAX(14*minmn+4,10*minmn+2+25*(25+8))+MAX(m,n))
   ! gesdd docs are mess: optimal turns out to be less than minimal in docs
   ! check(lwork>=(compute_uv?3*minmn*minmn+MAX(MAX(m,n),4*minmn*(minmn+1)):MAX(14*minmn+4,10*minmn+2+25*(25+8))+MAX(m,n))) :: lwork
   integer intent(hide,cache),dimension(8*minmn),depend(minmn) :: iwork
   integer intent(out)::info

   end subroutine <tchar=s,d>gesdd

   subroutine <tchar=c,z>gesdd(m,n,minmn,du,dvt,a,compute_uv,u,s,vt,work,rwork,lwork,iwork,info)
   
   ! u,s,vt,info = gesdd(a,compute_uv=1,lwork=..,overwrite_a=0)
   ! Compute the singular value decomposition (SVD):
   !   A = U * SIGMA * conjugate-transpose(V)
   ! A - M x N matrix
   ! U - M x M matrix
   ! SIGMA - M x N zero matrix with a main diagonal filled with min(M,N) 
   !               singular values  
   ! conjugate-transpose(V) - N x N matrix

   callstatement (*f2py_func)((compute_uv?"A":"N"),&m,&n,a,&m,s,u,&du,vt,&dvt,work,&lwork,rwork,iwork,&info)
   callprotoargument char*,int*,int*,<type_in_c>*,int*,<rtype_in_c>*,<type_in_c>*,int*,<type_in_c>*,int*,<type_in_c>*,int*,<rtype_in_c>*,int*,int*

   integer intent(in),optional,check(compute_uv==0||compute_uv==1):: compute_uv = 1
   integer intent(hide),depend(a):: m = shape(a,0)
   integer intent(hide),depend(a):: n = shape(a,1)
   integer intent(hide),depend(m,n):: minmn = MIN(m,n)
   integer intent(hide),depend(compute_uv,minmn) :: du = (compute_uv?m:1)
   integer intent(hide),depend(compute_uv,n) :: dvt = (compute_uv?n:1)
   <type_in> dimension(m,n),intent(in,copy) :: a
   <rtype_in> dimension(minmn),intent(out),depend(minmn) :: s
   <type_in> dimension(du,du),intent(out),depend(du) :: u
   <type_in> dimension(dvt,dvt),intent(out),depend(dvt) :: vt
   <type_in> dimension(lwork),intent(hide,cache),depend(lwork) :: work
   <rtype_in> dimension((compute_uv?5*minmn*minmn+7*minmn:5*minmn)),intent(hide,cache),depend(minmn,compute_uv) :: rwork
   integer optional,intent(in),depend(minmn,compute_uv) &
        :: lwork = (compute_uv?2*minmn*minmn+MAX(m,n)+2*minmn:2*minmn+MAX(m,n))
   check(lwork>=(compute_uv?2*minmn*minmn+MAX(m,n)+2*minmn:2*minmn+MAX(m,n))) :: lwork
   integer intent(hide,cache),dimension(8*minmn),depend(minmn) :: iwork
   integer intent(out)::info

   end subroutine <tchar=c,z>gesdd


   subroutine <tchar=s,d>gelss(m,n,minmn,maxmn,nrhs,a,b,s,cond,r,work,lwork,info)

   ! v,x,s,rank,info = gelss(a,b,cond=-1.0,overwrite_a=0,overwrite_b=0)
   ! Solve Minimize 2-norm(A * X - B).

     callstatement (*f2py_func)(&m,&n,&nrhs,a,&m,b,&maxmn,s,&cond,&r,work,&lwork,&info)
     callprotoargument int*,int*,int*,<type_in_c>*,int*,<type_in_c>*,int*,<type_in_c>*,<type_in_c>*,int*,<type_in_c>*,int*,int*

     integer intent(hide),depend(a):: m = shape(a,0)
     integer intent(hide),depend(a):: n = shape(a,1)
     integer intent(hide),depend(m,n):: minmn = MIN(m,n)
     integer intent(hide),depend(m,n):: maxmn = MAX(m,n)
     <type_in> dimension(m,n),intent(in,out,copy,out=v) :: a

     integer depend(b),intent(hide):: nrhs = shape(b,1)
     <type_in> dimension(maxmn,nrhs),check(maxmn==shape(b,0)),depend(maxmn) :: b
     intent(in,out,copy,out=x) b

     integer intent(out)::info

     integer optional,intent(in),depend(nrhs,minmn,maxmn), &
          check(lwork>=1) &
          :: lwork=3*minmn+MAX(2*minmn,MAX(maxmn,nrhs))
          !check(lwork>=3*minmn+MAX(2*minmn,MAX(maxmn,nrhs)))
     <type_in> dimension(lwork),intent(hide),depend(lwork) :: work

     <type_in> intent(in),optional :: cond = -1.0
     integer intent(out,out=rank) :: r
     <type_in> intent(out),dimension(minmn),depend(minmn) :: s

   end subroutine <tchar=s,d>gelss

   subroutine <tchar=c,z>gelss(m,n,minmn,maxmn,nrhs,a,b,s,cond,r,work,rwork,lwork,info)

   ! v,x,s,rank,info = gelss(a,b,cond=-1.0,overwrite_a=0,overwrite_b=0)
   ! Solve Minimize 2-norm(A * X - B).

     callstatement (*f2py_func)(&m,&n,&nrhs,a,&m,b,&maxmn,s,&cond,&r,work,&lwork,rwork,&info)
     callprotoargument int*,int*,int*,<type_in_c>*,int*,<type_in_c>*,int*,<rtype_in_c>*,<rtype_in_c>*,int*,<type_in_c>*,int*,<rtype_in_c>*,int*

     integer intent(hide),depend(a):: m = shape(a,0)
     integer intent(hide),depend(a):: n = shape(a,1)
     integer intent(hide),depend(m,n):: minmn = MIN(m,n)
     integer intent(hide),depend(m,n):: maxmn = MAX(m,n)
     <type_in> dimension(m,n),intent(in,out,copy,out=v) :: a

     integer depend(b),intent(hide):: nrhs = shape(b,1)
     <type_in> dimension(maxmn,nrhs),check(maxmn==shape(b,0)),depend(maxmn) :: b
     intent(in,out,copy,out=x) b

     integer intent(out)::info

     integer optional,intent(in),depend(nrhs,minmn,maxmn), &
          check(lwork>=1) &
          :: lwork=2*minmn+MAX(maxmn,nrhs)
          ! check(lwork>=2*minmn+MAX(maxmn,nrhs))
     <type_in> dimension(lwork),intent(hide),depend(lwork) :: work
     <rtype_in> dimension(5*minmn-1),intent(hide),depend(lwork) :: rwork

     <rtype_in> intent(in),optional :: cond = -1.0
     integer intent(out,out=rank) :: r
     <rtype_in> intent(out),dimension(minmn),depend(minmn) :: s

   end subroutine <tchar=c,z>gelss

   subroutine <tchar=s,d,c,z>geqrf(m,n,a,tau,work,lwork,info)

   ! qr_a,tau,work,info = geqrf(a,lwork=3*n,overwrite_a=0)
   ! Compute a QR factorization of a real M-by-N matrix A:
   !   A = Q * R.

     callstatement (*f2py_func)(&m,&n,a,&m,tau,work,&lwork,&info)
     callprotoargument int*,int*,<type_in_c>*,int*,<type_in_c>*,<type_in_c>*,int*,int*

     integer intent(hide),depend(a):: m = shape(a,0)
     integer intent(hide),depend(a):: n = shape(a,1)
     <type_in> dimension(m,n),intent(in,out,copy,out=qr) :: a
     <type_in> dimension(MIN(m,n)),intent(out) :: tau

     integer optional,intent(in),depend(n),check(lwork>=n||lwork==-1) :: lwork=3*n
     <type_in> dimension(MAX(lwork,1)),intent(out),depend(lwork) :: work
     integer intent(out) :: info
   end subroutine <tchar=s,d,c,z>geqrf

   subroutine <tchar=s,d,c,z>gerqf(m,n,a,tau,work,lwork,info)

   ! rq_a,tau,work,info = gerqf(a,lwork=3*n,overwrite_a=0)
   ! Compute an RQ factorization of a real M-by-N matrix A:
   !   A = R * Q.

     callstatement (*f2py_func)(&m,&n,a,&m,tau,work,&lwork,&info)
     callprotoargument int*,int*,<type_in_c>*,int*,<type_in_c>*,<type_in_c>*,int*,int*

     integer intent(hide),depend(a):: m = shape(a,0)
     integer intent(hide),depend(a):: n = shape(a,1)
     <type_in> dimension(m,n),intent(in,out,copy,out=qr) :: a
     <type_in> dimension(MIN(m,n)),intent(out) :: tau

     integer optional,intent(in),depend(n),check(lwork>=n||lwork==-1) :: lwork=3*n
     <type_in> dimension(MAX(lwork,1)),intent(out),depend(lwork) :: work
     integer intent(out) :: info
   end subroutine <tchar=s,d,c,z>gerqf

   subroutine <tchar=s,d>orgqr(m,n,k,a,tau,work,lwork,info)

   ! q,work,info = orgqr(a,lwork=3*n,overwrite_a=0)
   ! Generates an M-by-N real matrix Q with orthonormal columns,
   ! which is defined as the first N columns of a product of K elementary
   ! reflectors of order M (e.g. output of geqrf)

     callstatement (*f2py_func)(&m,&n,&k,a,&m,tau,work,&lwork,&info)
     callprotoargument int*,int*,int*,<type_in_c>*,int*,<type_in_c>*,<type_in_c>*,int*,int*

     integer intent(hide),depend(a):: m = shape(a,0)
     integer intent(hide),depend(a):: n = shape(a,1)
     integer intent(hide),depend(tau):: k = shape(tau,0)
     <type_in> dimension(m,n),intent(in,out,copy,out=q) :: a
     <type_in> dimension(k),intent(in) :: tau

     integer optional,intent(in),depend(n),check(lwork>=n||lwork==-1) :: lwork=3*n
     <type_in> dimension(MAX(lwork,1)),intent(out),depend(lwork) :: work
     integer intent(out) :: info
   end subroutine <tchar=s,d>orgqr

   subroutine <tchar=c,z>ungqr(m,n,k,a,tau,work,lwork,info)

   ! q,work,info = ungqr(a,lwork=3*n,overwrite_a=0)
   ! Generates an M-by-N complex matrix Q with unitary columns,
   ! which is defined as the first N columns of a product of K elementary
   ! reflectors of order M (e.g. output of geqrf)

     callstatement (*f2py_func)(&m,&n,&k,a,&m,tau,work,&lwork,&info)
     callprotoargument int*,int*,int*,<type_in_c>*,int*,<type_in_c>*,<type_in_c>*,int*,int*

     integer intent(hide),depend(a):: m = shape(a,0)
     integer intent(hide),depend(a):: n = shape(a,1)
     integer intent(hide),depend(tau):: k = shape(tau,0)
     <type_in> dimension(m,n),intent(in,out,copy,out=q) :: a
     <type_in> dimension(k),intent(in) :: tau

     integer optional,intent(in),depend(n),check(lwork>=n||lwork==-1) :: lwork=3*n
     <type_in> dimension(MAX(lwork,1)),intent(out),depend(lwork) :: work
     integer intent(out) :: info
   end subroutine <tchar=c,z>ungqr
+
   subroutine <tchar=s,d>geev(compute_vl,compute_vr,n,a,wr,wi,vl,ldvl,vr,ldvr,work,lwork,info)

     ! wr,wi,vl,vr,info = geev(a,compute_vl=1,compute_vr=1,lwork=4*n,overwrite_a=0)

     callstatement {(*f2py_func)((compute_vl?"V":"N"),(compute_vr?"V":"N"),&n,a,&n,wr,wi,vl,&ldvl,vr,&ldvr,work,&lwork,&info);}
     callprotoargument char*,char*,int*,<type_in_c>*,int*,<type_in_c>*,<type_in_c>*,<type_in_c>*,int*,<type_in_c>*,int*,<type_in_c>*,int*,int*

     integer optional,intent(in):: compute_vl = 1
     check(compute_vl==1||compute_vl==0) compute_vl
     integer optional,intent(in):: compute_vr = 1
     check(compute_vr==1||compute_vr==0) compute_vr

     integer intent(hide),depend(a) :: n = shape(a,0)
     <type_in>  dimension(n,n),intent(in,copy) :: a
     check(shape(a,0)==shape(a,1)) :: a

     <type_in>  dimension(n),intent(out),depend(n) :: wr
     <type_in>  dimension(n),intent(out),depend(n) :: wi

     <type_in>  dimension(ldvl,n),intent(out) :: vl
     integer intent(hide),depend(n,compute_vl) :: ldvl=(compute_vl?n:1)

     <type_in>  dimension(ldvr,n),intent(out) :: vr
     integer intent(hide),depend(n,compute_vr) :: ldvr=(compute_vr?n:1)

     integer optional,intent(in),depend(n,compute_vl,compute_vr) :: lwork=4*n
     check(lwork>=((compute_vl||compute_vr)?4*n:3*n)) :: lwork
     <type_in> dimension(lwork),intent(hide,cache),depend(lwork) :: work

     integer intent(out):: info
   end subroutine <tchar=s,d>geev

   subroutine <tchar=c,z>geev(compute_vl,compute_vr,n,a,w,vl,ldvl,vr,ldvr,work,lwork,rwork,info)

     ! w,vl,vr,info = geev(a,compute_vl=1,compute_vr=1,lwork=2*n,overwrite_a=0)

     callstatement (*f2py_func)((compute_vl?"V":"N"),(compute_vr?"V":"N"),&n,a,&n,w,vl,&ldvl,vr,&ldvr,work,&lwork,rwork,&info)
     callprotoargument char*,char*,int*,<type_in_c>*,int*,<type_in_c>*,<type_in_c>*,int*,<type_in_c>*,int*,<type_in_c>*,int*,<rtype_in_c>*,int*

     integer optional,intent(in):: compute_vl = 1
     check(compute_vl==1||compute_vl==0) compute_vl
     integer optional,intent(in):: compute_vr = 1
     check(compute_vr==1||compute_vr==0) compute_vr

     integer intent(hide),depend(a) :: n = shape(a,0)
     <type_in>  dimension(n,n),intent(in,copy) :: a
     check(shape(a,0)==shape(a,1)) :: a

     <type_in>  dimension(n),intent(out),depend(n) :: w

     <type_in>  dimension(ldvl,n),depend(ldvl),intent(out) :: vl
     integer intent(hide),depend(compute_vl,n) :: ldvl=(compute_vl?n:1)

     <type_in>  dimension(ldvr,n),depend(ldvr),intent(out) :: vr
     integer intent(hide),depend(compute_vr,n) :: ldvr=(compute_vr?n:1)
     
     integer optional,intent(in),depend(n) :: lwork=2*n
     check(lwork>=2*n) :: lwork
     <type_in> dimension(lwork),intent(hide),depend(lwork) :: work
     <rtype_in> dimension(2*n),intent(hide,cache),depend(n) :: rwork

     integer intent(out):: info
   end subroutine <tchar=c,z>geev

   subroutine <tchar=s,d>gegv(compute_vl,compute_vr,n,a,b,alphar,alphai,beta,vl,ldvl,vr,ldvr,work,lwork,info)
     ! Compute the generalized eigenvalues (alphar +/- alphai*i, beta)
     ! of the real nonsymmetric matrices A and B: det(A-w*B)=0 where w=alpha/beta.
     ! Optionally, compute the left and/or right generalized eigenvectors:
     ! (A - w B) r = 0, l^H * (A - w B) = 0
     !
     ! alphar,alphai,beta,vl,vr,info = gegv(a,b,compute_vl=1,compute_vr=1,lwork=8*n,overwrite_a=0,overwrite_b=0)

     callstatement (*f2py_func)((compute_vl?"V":"N"),(compute_vr?"V":"N"),&n,a,&n,b,&n,alphar,alphai,beta,vl,&ldvl,vr,&ldvr,work,&lwork,&info)
     callprotoargument char*,char*,int*,<type_in_c>*,int*,<type_in_c>*,int*,<type_in_c>*,<type_in_c>*,<type_in_c>*,<type_in_c>*,int*,<type_in_c>*,int*,<type_in_c>*,int*,int*

     integer optional,intent(in):: compute_vl = 1
     check(compute_vl==1||compute_vl==0) compute_vl
     integer optional,intent(in):: compute_vr = 1
     check(compute_vr==1||compute_vr==0) compute_vr

     integer intent(hide),depend(a) :: n = shape(a,0)
     <type_in>  dimension(n,n),intent(in,copy) :: a
     check(shape(a,0)==shape(a,1)) :: a

     <type_in>  dimension(n,n),depend(n),intent(in,copy) :: b
     check(shape(b,0)==shape(b,1) && shape(b,0)==n) :: b

     <type_in> dimension(n),depend(n),intent(out) :: alphar
     <type_in> dimension(n),depend(n),intent(out) :: alphai
     <type_in> dimension(n),depend(n),intent(out) :: beta

     <type_in>  dimension(ldvl,n),intent(out),depend(ldvl) :: vl
     integer intent(hide),depend(compute_vl,n) :: ldvl=(compute_vl?n:1)

     <type_in>  dimension(ldvr,n),intent(out),depend(ldvr) :: vr
     integer intent(hide),depend(compute_vr,n) :: ldvr=(compute_vr?n:1)

     integer optional,intent(in),depend(n) :: lwork=8*n
     check(lwork>=8*n) :: lwork
     <type_in> dimension(lwork),intent(hide),depend(lwork) :: work

     integer intent(out):: info

   end subroutine <tchar=s,d>gegv

   subroutine <tchar=c,z>gegv(compute_vl,compute_vr,n,a,b,alpha,beta,vl,ldvl,vr,ldvr,work,lwork,rwork,info)
     ! Compute the generalized eigenvalues (alpha, beta)
     ! of the comples nonsymmetric matrices A and B: det(A-w*B)=0 where w=alpha/beta.
     ! Optionally, compute the left and/or right generalized eigenvectors:
     ! (A - w B) r = 0, l^H * (A - w B) = 0
     !
     ! alpha,beta,vl,vr,info = gegv(a,b,compute_vl=1,compute_vr=1,lwork=2*n,overwrite_a=0,overwrite_b=0)

     callstatement (*f2py_func)((compute_vl?"V":"N"),(compute_vr?"V":"N"),&n,a,&n,b,&n,alpha,beta,vl,&ldvl,vr,&ldvr,work,&lwork,rwork,&info)
     callprotoargument char*,char*,int*,<type_in_c>*,int*,<type_in_c>*,int*,<type_in_c>*,<type_in_c>*,<type_in_c>*,int*,<type_in_c>*,int*,<type_in_c>*,int*,<rtype_in_c>*,int*

     integer optional,intent(in):: compute_vl = 1
     check(compute_vl==1||compute_vl==0) compute_vl
     integer optional,intent(in):: compute_vr = 1
     check(compute_vr==1||compute_vr==0) compute_vr

     integer intent(hide),depend(a) :: n = shape(a,0)
     <type_in>  dimension(n,n),intent(in,copy) :: a
     check(shape(a,0)==shape(a,1)) :: a

     <type_in>  dimension(n,n),depend(n),intent(in,copy) :: b
     check(shape(b,0)==shape(b,1) && shape(b,0)==n) :: b

     <type_in> dimension(n),depend(n),intent(out) :: alpha
     <type_in> dimension(n),depend(n),intent(out) :: beta

     <type_in>  dimension(ldvl,n),intent(out),depend(ldvl) :: vl
     integer intent(hide),depend(compute_vl,n) :: ldvl=(compute_vl?n:1)

     <type_in>  dimension(ldvr,n),intent(out),depend(ldvr) :: vr
     integer intent(hide),depend(compute_vr,n) :: ldvr=(compute_vr?n:1)

     integer optional,intent(in),depend(n) :: lwork=2*n
     check(lwork>=2*n) :: lwork
     <type_in> dimension(lwork),intent(hide),depend(lwork) :: work
     <rtype_in> dimension(8*n),intent(hide),depend(n) :: rwork

     integer intent(out):: info

   end subroutine <tchar=c,z>gegv


   subroutine <tchar=s,d>syev(compute_v,lower,n,w,a,work,lwork,info)

   ! w,v,info = syev(a,compute_v=1,lower=0,lwork=3*n-1,overwrite_a=0)
   ! Compute all eigenvalues and, optionally, eigenvectors of a
   ! real symmetric matrix A.
   !
   ! Performance tip:
   !   If compute_v=0 then set also overwrite_a=1.

     callstatement (*f2py_func)((compute_v?"V":"N"),(lower?"L":"U"),&n,a,&n,w,work,&lwork,&info)
     callprotoargument char*,char*,int*,<type_in_c>*,int*,<type_in_c>*,<type_in_c>*,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 intent(hide),depend(a):: n = shape(a,0)
     <type_in> dimension(n,n),check(shape(a,0)==shape(a,1)) :: a
     intent(in,copy,out,out=v) :: a

     <type_in> dimension(n),intent(out),depend(n) :: w

     integer optional,intent(in),depend(n) :: lwork=3*n-1
     check(lwork>=3*n-1) :: lwork
     <type_in> dimension(lwork),intent(hide),depend(lwork) :: work

     integer intent(out) :: info
   end subroutine <tchar=s,d>syev

   subroutine <tchar=c,z>heev(compute_v,lower,n,w,a,work,lwork,rwork,info)

   ! w,v,info = syev(a,compute_v=1,lower=0,lwork=3*n-1,overwrite_a=0)
   ! Compute all eigenvalues and, optionally, eigenvectors of a
   ! complex Hermitian matrix A.
   !
   ! Warning:
   !   If compute_v=0 and overwrite_a=1, the contents of a is destroyed.

     callstatement (*f2py_func)((compute_v?"V":"N"),(lower?"L":"U"),&n,a,&n,w,work,&lwork,rwork,&info)
     callprotoargument char*,char*,int*,<type_in_c>*,int*,<type_in_c>*,<type_in_c>*,int*,<rtype_in_c>*,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 intent(hide),depend(a):: n = shape(a,0)
     <type_in> dimension(n,n),check(shape(a,0)==shape(a,1)) :: a
     intent(in,copy,out,out=v) :: a

     <rtype_in> dimension(n),intent(out),depend(n) :: w

     integer optional,intent(in),depend(n) :: lwork=2*n-1
     check(lwork>=2*n-1) :: lwork
     <type_in> dimension(lwork),intent(hide),depend(lwork) :: work

     <rtype_in> dimension(3*n-1),intent(hide),depend(n) :: rwork

     integer intent(out) :: info
   end subroutine <tchar=c,z>heev

   subroutine <tchar=s,d,c,z>posv(n,nrhs,a,b,info,lower)

   ! c,x,info = posv(a,b,lower=0,overwrite_a=0,overwrite_b=0)
   ! Solve A * X = B.
   ! A is symmetric positive defined
   ! A = U^T * U, C = U if lower = 0
   ! A = L * L^T, C = L if lower = 1
   ! C is triangular matrix of the corresponding Cholesky decomposition.

     callstatement (*f2py_func)((lower?"L":"U"),&n,&nrhs,a,&n,b,&n,&info)
     callprotoargument char*,int*,int*,<type_in_c>*,int*,<type_in_c>*,int*,int*

     integer optional,intent(in),check(lower==0||lower==1) :: lower = 0

     integer depend(a),intent(hide):: n = shape(a,0)
     integer depend(b),intent(hide):: nrhs = shape(b,1)
     <type_in> dimension(n,n),intent(in,out,copy,out=c) :: a
     check(shape(a,0)==shape(a,1)) :: a
     <type_in> dimension(n,nrhs),intent(in,out,copy,out=x),depend(n):: b
     check(shape(a,0)==shape(b,0)) :: b
     integer intent(out) :: info

   end subroutine <tchar=s,d,c,z>posv
        
   subroutine <tchar=s,d>potrf(n,a,info,lower,clean)
   
     ! c,info = potrf(a,lower=0,clean=1,overwrite_a=0)
     ! Compute Cholesky decomposition of symmetric positive defined matrix:
     ! A = U^T * U, C = U if lower = 0
     ! A = L * L^T, C = L if lower = 1
     ! C is triangular matrix of the corresponding Cholesky decomposition.
     ! clean==1 zeros strictly lower or upper parts of U or L, respectively

     callstatement (*f2py_func)((lower?"L":"U"),&n,a,&n,&info); if(clean){int i,j;if(lower){for(i=0;i<n;++i) for(j=i+1;j<n;++j) *(a+j*n+i)=0.0;} else {for(i=0;i<n;++i) for(j=i+1;j<n;++j) *(a+i*n+j)=0.0;}}
     callprotoargument char*,int*,<type_in_c>*,int*,int*

     integer optional,intent(in),check(lower==0||lower==1) :: lower = 0
     integer optional,intent(in),check(clean==0||clean==1) :: clean = 1
     integer depend(a),intent(hide):: n = shape(a,0)
     <type_in> dimension(n,n),intent(in,out,copy,out=c) :: a
     check(shape(a,0)==shape(a,1)) :: a
     integer intent(out) :: info
     
   end subroutine <tchar=s,d>potrf

   subroutine <tchar=c,z>potrf(n,a,info,lower,clean)
   
     ! c,info = potrf(a,lower=0,clean=1,overwrite_a=0)
     ! Compute Cholesky decomposition of symmetric positive defined matrix:
     ! A = U^H * U, C = U if lower = 0
     ! A = L * L^H, C = L if lower = 1
     ! C is triangular matrix of the corresponding Cholesky decomposition.
     ! clean==1 zeros strictly lower or upper parts of U or L, respectively

     callstatement (*f2py_func)((lower?"L":"U"),&n,a,&n,&info); if(clean){int i,j,k;if(lower){for(i=0;i<n;++i) for(j=i+1;j<n;++j) {k=j*n+i;(a+k)->r=(a+k)->i=0.0;}} else {for(i=0;i<n;++i) for(j=i+1;j<n;++j) {k=i*n+j;(a+k)->r=(a+k)->i=0.0;}}}
     callprotoargument char*,int*,<type_in_c>*,int*,int*

     integer optional,intent(in),check(lower==0||lower==1) :: lower = 0
     integer optional,intent(in),check(clean==0||clean==1) :: clean = 1
     integer depend(a),intent(hide):: n = shape(a,0)
     <type_in> dimension(n,n),intent(in,out,copy,out=c) :: a
     check(shape(a,0)==shape(a,1)) :: a
     integer intent(out) :: info
     
   end subroutine <tchar=c,z>potrf

   subroutine <tchar=s,d,c,z>potrs(n,nrhs,c,b,info,lower)

   ! x,info = potrs(c,b,lower=0=1,overwrite_b=0)
   ! Solve A * X = B.
   ! A is symmetric positive defined
   ! A = U^T * U, C = U if lower = 0
   ! A = L * L^T, C = L if lower = 1
   ! C is triangular matrix of the corresponding Cholesky decomposition.

     callstatement (*f2py_func)((lower?"L":"U"),&n,&nrhs,c,&n,b,&n,&info)
     callprotoargument char*,int*,int*,<type_in_c>*,int*,<type_in_c>*,int*,int*

     integer optional,intent(in),check(lower==0||lower==1) :: lower = 0

     integer depend(c),intent(hide):: n = shape(c,0)
     integer depend(b),intent(hide):: nrhs = shape(b,1)
     <type_in> dimension(n,n),intent(in) :: c
     check(shape(c,0)==shape(c,1)) :: c
     <type_in> dimension(n,nrhs),intent(in,out,copy,out=x),depend(n):: b
     check(shape(c,0)==shape(b,0)) :: b
     integer intent(out) :: info

   end subroutine <tchar=s,d,c,z>potrs

   subroutine <tchar=s,d,c,z>potri(n,c,info,lower)
   
     ! inv_a,info = potri(c,lower=0,overwrite_c=0)
     ! Compute A inverse A^-1.
     ! A = U^T * U, C = U if lower = 0
     ! A = L * L^T, C = L if lower = 1
     ! C is triangular matrix of the corresponding Cholesky decomposition.

     callstatement (*f2py_func)((lower?"L":"U"),&n,c,&n,&info)
     callprotoargument char*,int*,<type_in_c>*,int*,int*

     integer optional,intent(in),check(lower==0||lower==1) :: lower = 0
     
     integer depend(c),intent(hide):: n = shape(c,0)
     <type_in> dimension(n,n),intent(c,in,out,copy,out=inv_a) :: c
     check(shape(c,0)==shape(c,1)) :: c
     integer intent(out) :: info
     
   end subroutine <tchar=s,d,c,z>potri


   subroutine <tchar=s,d,c,z>lauum(n,c,info,lower)
   
     ! a,info = lauum(c,lower=0,overwrite_c=0)
     ! Compute product
     ! U^T * U, C = U if lower = 0
     ! L * L^T, C = L if lower = 1
     ! C is triangular matrix of the corresponding Cholesky decomposition.

     callstatement (*f2py_func)((lower?"L":"U"),&n,c,&n,&info)
     callprotoargument char*,int*,<type_in_c>*,int*,int*

     integer optional,intent(in),check(lower==0||lower==1) :: lower = 0
     
     integer depend(c),intent(hide):: n = shape(c,0)
     <type_in> dimension(n,n),intent(in,out,copy,out=a) :: c
     check(shape(c,0)==shape(c,1)) :: c
     integer intent(out) :: info
     
   end subroutine <tchar=s,d,c,z>lauum

   subroutine <tchar=s,d,c,z>trtri(n,c,info,lower,unitdiag)
   
     ! inv_c,info = trtri(c,lower=0,unitdiag=1,overwrite_c=0)
     ! Compute C inverse C^-1 where
     ! C = U if lower = 0
     ! C = L if lower = 1
     ! C is non-unit triangular matrix if unitdiag = 0
     ! C is unit triangular matrix if unitdiag = 1

     callstatement (*f2py_func)((lower?"L":"U"),(unitdiag?"U":"N"),&n,c,&n,&info)
     callprotoargument char*,char*,int*,<type_in_c>*,int*,int*

     integer optional,intent(in),check(lower==0||lower==1) :: lower = 0
     integer optional,intent(in),check(unitdiag==0||unitdiag==1) :: unitdiag = 0
     
     integer depend(c),intent(hide):: n = shape(c,0)
     <type_in> dimension(n,n),intent(in,out,copy,out=inv_c) :: c
     check(shape(c,0)==shape(c,1)) :: c
     integer intent(out) :: info
     
   end subroutine <tchar=s,d,c,z>trtri


   subroutine <tchar=s,d,c,z>laswp(n,a,nrows,k1,k2,piv,off,inc,m)

   ! a = laswp(a,piv,k1=0,k2=len(piv)-1,off=0,inc=1,overwrite_a=0)
   ! Perform row interchanges on the matrix A for each of row k1 through k2
   ! 
   ! piv pivots rows.

     callstatement {int i;m=len(piv);for(i=0;i<m;++piv[i++]);++k1;++k2; (*f2py_func)(&n,a,&nrows,&k1,&k2,piv+off,&inc); for(i=0;i<m;--piv[i++]);}
     callprotoargument int*,<type_in_c>*,int*,int*,int*,int*,int*

     integer depend(a),intent(hide):: nrows = shape(a,0)
     integer depend(a),intent(hide):: n = shape(a,1)
     <type_in> dimension(nrows,n),intent(in,out,copy) :: a
     integer dimension(*),intent(in),depend(nrows) :: piv
     check(len(piv)<=nrows) :: piv
!XXX: how to check that all elements in piv are < n?

     integer optional,intent(in) :: k1 = 0
     check(0<=k1) :: k1
     integer optional,intent(in),depend(k1,piv,off) :: k2 = len(piv)-1
     check(k1<=k2 && k2<len(piv)-off) :: k2

     integer optional, intent(in),check(inc>0||inc<0) :: inc = 1
     integer optional,intent(in),depend(piv) :: off=0
     check(off>=0 && off<len(piv)) :: off

     integer intent(hide),depend(piv,inc,off) :: m = (len(piv)-off)/abs(inc)
     check(len(piv)-off>(m-1)*abs(inc)) :: m

   end subroutine <tchar=s,d,c,z>laswp

   subroutine <tchar=c,z>gees(compute_v,sort_t,<tchar>select,n,a,nrows,sdim,w,vs,ldvs,work,lwork,rwork,bwork,info)

     ! t,sdim,w,vs,work,info=gees(compute_v=1,sort_t=0,select,a,lwork=3*n)
     ! For an NxN matrix compute the eigenvalues, the schur form T, and optionally
     !  the matrix of Schur vectors Z.  This gives the Schur factorization 
     !  A = Z * T * Z^H  -- a complex matrix is in Schur form if it is upper 
     !  triangular

     callstatement (*f2py_func)((compute_v?"V":"N"),(sort_t?"S":"N"),cb_<tchar>select_in_<tchar>gees__user__routines,&n,a,&nrows,&sdim,w,vs,&ldvs,work,&lwork,rwork,bwork,&info,1,1)
     callprotoargument char*,char*,int(*)(<type_in_c>*),int*,<type_in_c>*,int*,int*,<type_in_c>*,<type_in_c>*,int*,<type_in_c>*,int*,<rtype_in_c>*,int*,int*,int,int

     use <tchar>gees__user__routines

     integer optional,intent(in),check(compute_v==0||compute_v==1) :: compute_v = 1
     integer optional,intent(in),check(sort_t==0||sort_t==1) :: sort_t = 0
     external <tchar>select
     integer intent(hide),depend(a) :: n = shape(a,1)
     <type_in> intent(in,out,copy,out=t),check(shape(a,0)==shape(a,1)),dimension(n,n) :: a
     integer intent(hide),depend(a) :: nrows=shape(a,0)
     integer intent(out) :: sdim=0
     <type_in> intent(out),dimension(n) :: w
     <type_in> intent(out),depend(ldvs,n),dimension(ldvs,n) :: vs
     integer intent(hide),depend(compute_v,n) :: ldvs=((compute_v==1)?n:1)
     <type_in> intent(out),depend(lwork),dimension(MAX(lwork,1)) :: work
     integer optional,intent(in),check((lwork==-1)||(lwork >= MAX(1,2*n))),depend(n) :: lwork = 3*n
     <rtype_in> optional,intent(hide),depend(n),dimension(n) :: rwork
     logical optional,intent(hide),depend(n),dimension(n) :: bwork
     integer intent(out) :: info
   end subroutine <tchar=c,z>gees

   subroutine <tchar=d,s>gees(compute_v,sort_t,<tchar>select,n,a,nrows,sdim,wr,wi,vs,ldvs,work,lwork,bwork,info)

     ! t,sdim,w,vs,work,info=gees(compute_v=1,sort_t=0,select,a,lwork=3*n)
     ! For an NxN matrix compute the eigenvalues, the schur form T, and optionally
     !  the matrix of Schur vectors Z.  This gives the Schur factorization 
     !  A = Z * T * Z^H  -- a real matrix is in Schur form if it is upper quasi-
     !  triangular with 1x1 and 2x2 blocks.

     callstatement (*f2py_func)((compute_v?"V":"N"),(sort_t?"S":"N"),cb_<tchar>select_in_<tchar>gees__user__routines,&n,a,&nrows,&sdim,wr,wi,vs,&ldvs,work,&lwork,bwork,&info,1,1)
     callprotoargument char*,char*,int(*)(<type_in_c>*,<type_in_c>*),int*,<type_in_c>*,int*,int*,<type_in_c>*,<type_in_c>*,<type_in_c>*,int*,<type_in_c>*,int*,int*,int*,int,int

    use <tchar>gees__user__routines

    integer optional,intent(in),check(compute_v==0||compute_v==1) :: compute_v = 1
    integer optional,intent(in),check(sort_t==0||sort_t==1) :: sort_t = 0
    external <tchar>select
    integer intent(hide),depend(a) :: n = shape(a,1)
    <type_in> intent(in,out,copy,out=t),check(shape(a,0)==shape(a,1)),dimension(n,n) :: a
    integer intent(hide),depend(a) :: nrows=shape(a,0)
    integer intent(out) :: sdim=0
    <type_in> intent(out),dimension(n) :: wr
    <type_in> intent(out),dimension(n) :: wi
    <type_in> intent(out),depend(ldvs,n),dimension(ldvs,n) :: vs
    integer intent(hide),depend(compute_v,n) :: ldvs=((compute_v==1)?n:1)
    <type_in> intent(out),depend(lwork),dimension(MAX(lwork,1)) :: work
    integer optional,intent(in),check((lwork==-1)||(lwork >= MAX(1,2*n))),depend(n) :: lwork = 3*n
    <rtype_in> optional,intent(hide),depend(n),dimension(n) :: rwork
    logical optional,intent(hide),depend(n),dimension(n) :: bwork
    integer intent(out) :: info
end subroutine <tchar=s,d>gees


subroutine <tchar=s,d>ggev(compute_vl,compute_vr,n,a,b,alphar,alphai,beta,vl,ldvl,vr,ldvr,work,lwork,info) 

     callstatement {(*f2py_func)((compute_vl?"V":"N"),(compute_vr?"V":"N"),&n,a,&n,b,&n,alphar,alphai,beta,vl,&ldvl,vr,&ldvr,work,&lwork,&info);}
     callprotoargument char*,char*,int*,<type_in_c>*,int*,<type_in_c>*,int*,<type_in_c>*,<type_in_c>*,<type_in_c>*,<type_in_c>*,int*,<type_in_c>*,int*,<type_in_c>*,int*,int*

    integer optional,intent(in):: compute_vl = 1
    check(compute_vl==1||compute_vl==0) compute_vl
    integer optional,intent(in):: compute_vr = 1
    check(compute_vr==1||compute_vr==0) compute_vr

    integer intent(hide),depend(a) :: n = shape(a,0)
    <type_in>  dimension(n,n),intent(in,copy) :: a
    check(shape(a,0)==shape(a,1)) :: a

    <type_in> intent(in,copy), dimension(n,n) :: b
    check(shape(b,0)==shape(b,1)) :: b

    <type_in> intent(out), dimension(n), depend(n) :: alphar
    <type_in> intent(out), dimension(n), depend(n) :: alphai
    <type_in> intent(out), dimension(n), depend(n) :: beta

    <type_in>  depend(ldvl,n), dimension(ldvl,n),intent(out) :: vl
    integer intent(hide),depend(n,compute_vl) :: ldvl=(compute_vl?n:1)
    
    <type_in>  depend(ldvr,n), dimension(ldvr,n),intent(out) :: vr
    integer intent(hide),depend(n,compute_vr) :: ldvr=(compute_vr?n:1)
    
    integer optional,intent(in),depend(n,compute_vl,compute_vr) :: lwork=8*n
    check((lwork==-1) || (lwork>=MAX(1,8*n))) :: lwork
    <type_in> intent(out), dimension(MAX(lwork,1)), depend(lwork) :: work

    integer intent(out):: info

end subroutine <tchar>ggev

subroutine <tchar=c,z>ggev(compute_vl,compute_vr,n,a,b,alpha,beta,vl,ldvl,vr,ldvr,work,lwork,rwork,info) 

     callstatement {(*f2py_func)((compute_vl?"V":"N"),(compute_vr?"V":"N"),&n,a,&n,b,&n,alpha,beta,vl,&ldvl,vr,&ldvr,work,&lwork,rwork,&info);}
     callprotoargument char*,char*,int*,<type_in_c>*,int*,<type_in_c>*,int*,<type_in_c>*,<type_in_c>*,<type_in_c>*,int*,<type_in_c>*,int*,<type_in_c>*,int*,<rtype_in_c>*,int*

    integer optional,intent(in):: compute_vl = 1
    check(compute_vl==1||compute_vl==0) compute_vl
    integer optional,intent(in):: compute_vr = 1
    check(compute_vr==1||compute_vr==0) compute_vr

    integer intent(hide),depend(a) :: n = shape(a,0)
    <type_in> dimension(n,n),intent(in,copy) :: a
    check(shape(a,0)==shape(a,1)) :: a

    <type_in> intent(in,copy), dimension(n,n) :: b
    check(shape(b,0)==shape(b,1)) :: b

    <type_in> intent(out), dimension(n), depend(n) :: alpha
    <type_in> intent(out), dimension(n), depend(n) :: beta

    <type_in>  depend(ldvl,n), dimension(ldvl,n),intent(out) :: vl
    integer intent(hide),depend(n,compute_vl) :: ldvl=(compute_vl?n:1)
    
    <type_in>  depend(ldvr,n), dimension(ldvr,n),intent(out) :: vr
    integer intent(hide),depend(n,compute_vr) :: ldvr=(compute_vr?n:1)
    
    integer optional,intent(in),depend(n,compute_vl,compute_vr) :: lwork=2*n
    check((lwork==-1) || (lwork>=MAX(1,2*n))) :: lwork
    <type_in> intent(out), dimension(MAX(lwork,1)), depend(lwork) :: work
    <rtype_in> intent(hide), dimension(8*n), depend(n) :: rwork

    integer intent(out):: info

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

!
! 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=18*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=18*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=18*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=18*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=18*n-1
    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=18*n-1
    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

! This file was auto-generated with f2py (version:2.10.173).
! See http://cens.ioc.ee/projects/f2py2e/