dtbmv(3f) - [BLAS:DOUBLE_BLAS_LEVEL3]
subroutine dtbmv(uplo,trans,diag,n,k,a,lda,x,incx)
.. Scalar Arguments ..
integer,intent(in) :: incx,k,lda,n
character,intent(in) :: diag,trans,uplo
..
.. Array Arguments ..
double precision,intent(in) :: a(lda,*)
double precision,intent(inout) :: x(*)
..
DTBMV performs one of the matrix-vector operations
x := A*x, or x := A**T*x,
where x is an n element vector and A is an n by n unit, or non-unit, upper or lower triangular band matrix, with ( k + 1 ) diagonals.
UPLO
UPLO is CHARACTER*1
On entry, UPLO specifies whether the matrix is an upper or
lower triangular matrix as follows:
UPLO = 'U' or 'u' A is an upper triangular matrix.
UPLO = 'L' or 'l' A is a lower triangular matrix.
TRANS
TRANS is CHARACTER*1
On entry, TRANS specifies the operation to be performed as
follows:
TRANS = 'N' or 'n' x := A*x.
TRANS = 'T' or 't' x := A**T*x.
TRANS = 'C' or 'c' x := A**T*x.
DIAG
DIAG is CHARACTER*1
On entry, DIAG specifies whether or not A is unit
triangular as follows:
DIAG = 'U' or 'u' A is assumed to be unit triangular.
DIAG = 'N' or 'n' A is not assumed to be unit
triangular.
N
N is INTEGER
On entry, N specifies the order of the matrix A.
N must be at least zero.
K
K is INTEGER
On entry with UPLO = 'U' or 'u', K specifies the number of
super-diagonals of the matrix A.
On entry with UPLO = 'L' or 'l', K specifies the number of
sub-diagonals of the matrix A.
K must satisfy 0 .le. K.
A
A is DOUBLE PRECISION array, dimension ( LDA, N )
Before entry with UPLO = 'U' or 'u', the leading ( k + 1 )
by n part of the array A must contain the upper triangular
band part of the matrix of coefficients, supplied column by
column, with the leading diagonal of the matrix in row
( k + 1 ) of the array, the first super-diagonal starting at
position 2 in row k, and so on. The top left k by k triangle
of the array A is not referenced.
The following program segment will transfer an upper
triangular band matrix from conventional full matrix storage
to band storage:
DO 20, J = 1, N
M = K + 1 - J
DO 10, I = MAX( 1, J - K ), J
A( M + I, J ) = matrix( I, J )
10 CONTINUE
20 CONTINUE
Before entry with UPLO = 'L' or 'l', the leading ( k + 1 )
by n part of the array A must contain the lower triangular
band part of the matrix of coefficients, supplied column by
column, with the leading diagonal of the matrix in row 1 of
the array, the first sub-diagonal starting at position 1 in
row 2, and so on. The bottom right k by k triangle of the
array A is not referenced.
The following program segment will transfer a lower
triangular band matrix from conventional full matrix storage
to band storage:
DO 20, J = 1, N
M = 1 - J
DO 10, I = J, MIN( N, J + K )
A( M + I, J ) = matrix( I, J )
10 CONTINUE
20 CONTINUE
Note that when DIAG = 'U' or 'u' the elements of the array A
corresponding to the diagonal elements of the matrix are not
referenced, but are assumed to be unity.
LDA
LDA is INTEGER
On entry, LDA specifies the first dimension of A as declared
in the calling (sub) program. LDA must be at least
( k + 1 ).
X
X is DOUBLE PRECISION array, dimension at least
( 1 + ( n - 1 )*abs( INCX ) ).
Before entry, the incremented array X must contain the n
element vector x. On exit, X is overwritten with the
transformed vector x.
INCX
INCX is INTEGER
On entry, INCX specifies the increment for the elements of
X. INCX must not be zero.
date:December 2016
\ingroup double_blas_level2
FURTHER DETAILS
Level 2 Blas routine. The vector and matrix arguments are not referenced when N = 0, or M = 0
– Written on 22-October-1986. Jack Dongarra, Argonne National Lab. Jeremy Du Croz, Nag Central Office. Sven Hammarling, Nag Central Office. Richard Hanson, Sandia National Labs.
Online html documentation available at
http://www.netlib.org/lapack/explore-html/
Type | Intent | Optional | Attributes | Name | ||
---|---|---|---|---|---|---|
character(len=1), | intent(in) | :: | uplo | |||
character(len=1), | intent(in) | :: | trans | |||
character(len=1), | intent(in) | :: | diag | |||
integer, | intent(in) | :: | n | |||
integer, | intent(in) | :: | k | |||
double precision, | intent(in) | :: | a(lda,*) | |||
integer, | intent(in) | :: | lda | |||
double precision, | intent(inout) | :: | x(*) | |||
integer, | intent(in) | :: | incx |
Type | Visibility | Attributes | Name | Initial | |||
---|---|---|---|---|---|---|---|
integer, | public | :: | i | ||||
integer, | public | :: | info | ||||
integer, | public | :: | ix | ||||
integer, | public | :: | j | ||||
integer, | public | :: | jx | ||||
integer, | public | :: | kplus1 | ||||
integer, | public | :: | kx | ||||
integer, | public | :: | l | ||||
logical, | public | :: | nounit | ||||
double precision, | public | :: | temp | ||||
double precision, | public, | parameter | :: | zero | = | 0.0d+0 |
subroutine dtbmv(uplo,trans,diag,n,k,a,lda,x,incx)
implicit none
!
! -- Reference BLAS level2 routine (version 3.7.0) --
! -- Reference BLAS is a software package provided by Univ. of Tennessee, --
! -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! December 2016
!
! .. Scalar Arguments ..
integer,intent(in) :: incx,k,lda,n
character,intent(in) :: diag,trans,uplo
! ..
! .. Array Arguments ..
double precision,intent(in) :: a(lda,*)
double precision,intent(inout) :: x(*)
! ..
!
! =====================================================================
!
! .. Parameters ..
double precision zero
parameter (zero=0.0d+0)
! ..
! .. Local Scalars ..
double precision temp
integer i,info,ix,j,jx,kplus1,kx,l
logical nounit
! ..
! .. External Functions ..
! ..
! .. External Subroutines ..
! ..
! .. Intrinsic Functions ..
intrinsic max,min
! ..
!
! Test the input parameters.
!
info = 0
if (.not.lsame(uplo,'U') .and. .not.lsame(uplo,'L')) then
info = 1
elseif (.not.lsame(trans,'N') .and. .not.lsame(trans,'T') .and. .not.lsame(trans,'C')) then
info = 2
elseif (.not.lsame(diag,'U') .and. .not.lsame(diag,'N')) then
info = 3
elseif (n.lt.0) then
info = 4
elseif (k.lt.0) then
info = 5
elseif (lda.lt. (k+1)) then
info = 7
elseif (incx.eq.0) then
info = 9
endif
if (info.ne.0) then
call xerbla('DTBMV ',info)
return
endif
!
! Quick return if possible.
!
if (n.eq.0) return
!
nounit = lsame(diag,'N')
!
! Set up the start point in X if the increment is not unity. This
! will be ( N - 1 )*INCX too small for descending loops.
!
if (incx.le.0) then
kx = 1 - (n-1)*incx
elseif (incx.ne.1) then
kx = 1
endif
!
! Start the operations. In this version the elements of A are
! accessed sequentially with one pass through A.
!
if (lsame(trans,'N')) then
!
! Form x := A*x.
!
if (lsame(uplo,'U')) then
kplus1 = k + 1
if (incx.eq.1) then
do j = 1,n
if (x(j).ne.zero) then
temp = x(j)
l = kplus1 - j
do i = max(1,j-k),j - 1
x(i) = x(i) + temp*a(l+i,j)
enddo
if (nounit) x(j) = x(j)*a(kplus1,j)
endif
enddo
else
jx = kx
do j = 1,n
if (x(jx).ne.zero) then
temp = x(jx)
ix = kx
l = kplus1 - j
do i = max(1,j-k),j - 1
x(ix) = x(ix) + temp*a(l+i,j)
ix = ix + incx
enddo
if (nounit) x(jx) = x(jx)*a(kplus1,j)
endif
jx = jx + incx
if (j.gt.k) kx = kx + incx
enddo
endif
else
if (incx.eq.1) then
do j = n,1,-1
if (x(j).ne.zero) then
temp = x(j)
l = 1 - j
do i = min(n,j+k),j + 1,-1
x(i) = x(i) + temp*a(l+i,j)
enddo
if (nounit) x(j) = x(j)*a(1,j)
endif
enddo
else
kx = kx + (n-1)*incx
jx = kx
do j = n,1,-1
if (x(jx).ne.zero) then
temp = x(jx)
ix = kx
l = 1 - j
do i = min(n,j+k),j + 1,-1
x(ix) = x(ix) + temp*a(l+i,j)
ix = ix - incx
enddo
if (nounit) x(jx) = x(jx)*a(1,j)
endif
jx = jx - incx
if ((n-j).ge.k) kx = kx - incx
enddo
endif
endif
else
!
! Form x := A**T*x.
!
if (lsame(uplo,'U')) then
kplus1 = k + 1
if (incx.eq.1) then
do j = n,1,-1
temp = x(j)
l = kplus1 - j
if (nounit) temp = temp*a(kplus1,j)
do i = j - 1,max(1,j-k),-1
temp = temp + a(l+i,j)*x(i)
enddo
x(j) = temp
enddo
else
kx = kx + (n-1)*incx
jx = kx
do j = n,1,-1
temp = x(jx)
kx = kx - incx
ix = kx
l = kplus1 - j
if (nounit) temp = temp*a(kplus1,j)
do i = j - 1,max(1,j-k),-1
temp = temp + a(l+i,j)*x(ix)
ix = ix - incx
enddo
x(jx) = temp
jx = jx - incx
enddo
endif
else
if (incx.eq.1) then
do j = 1,n
temp = x(j)
l = 1 - j
if (nounit) temp = temp*a(1,j)
do i = j + 1,min(n,j+k)
temp = temp + a(l+i,j)*x(i)
enddo
x(j) = temp
enddo
else
jx = kx
do j = 1,n
temp = x(jx)
kx = kx + incx
ix = kx
l = 1 - j
if (nounit) temp = temp*a(1,j)
do i = j + 1,min(n,j+k)
temp = temp + a(l+i,j)*x(ix)
ix = ix + incx
enddo
x(jx) = temp
jx = jx + incx
enddo
endif
endif
endif
!
end subroutine dtbmv