dtrsm(3f) - [BLAS:DOUBLE_BLAS_LEVEL3]
subroutine dtrsm(side,uplo,transa,diag,m,n,alpha,a,lda,b,ldb)
.. Scalar Arguments ..
double precision,intent(in) :: alpha
integer,intent(in) :: lda,ldb,m,n
character,intent(in) :: diag,side,transa,uplo
..
.. Array Arguments ..
double precision,intent(in) :: a(lda,*)
double precision,intent(inout) :: b(ldb,*)
..
DTRSM solves one of the matrix equations
op( A )*X = alpha*B, or X*op( A ) = alpha*B,
where alpha is a scalar, X and B are m by n matrices, A is a unit, or non-unit, upper or lower triangular matrix and op( A ) is one of
op( A ) = A or op( A ) = A**T.
The matrix X is overwritten on B.
SIDE
SIDE is CHARACTER*1
On entry, SIDE specifies whether op( A ) appears on the left
or right of X as follows:
SIDE = 'L' or 'l' op( A )*X = alpha*B.
SIDE = 'R' or 'r' X*op( A ) = alpha*B.
UPLO
UPLO is CHARACTER*1
On entry, UPLO specifies whether the matrix A 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.
TRANSA
TRANSA is CHARACTER*1
On entry, TRANSA specifies the form of op( A ) to be used in
the matrix multiplication as follows:
TRANSA = 'N' or 'n' op( A ) = A.
TRANSA = 'T' or 't' op( A ) = A**T.
TRANSA = 'C' or 'c' op( A ) = A**T.
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.
M
M is INTEGER
On entry, M specifies the number of rows of B. M must be at
least zero.
N
N is INTEGER
On entry, N specifies the number of columns of B. N must be
at least zero.
ALPHA
ALPHA is DOUBLE PRECISION.
On entry, ALPHA specifies the scalar alpha. When alpha is
zero then A is not referenced and B need not be set before
entry.
A
A is DOUBLE PRECISION array, dimension ( LDA, k ),
where k is m when SIDE = 'L' or 'l'
and k is n when SIDE = 'R' or 'r'.
Before entry with UPLO = 'U' or 'u', the leading k by k
upper triangular part of the array A must contain the upper
triangular matrix and the strictly lower triangular part of
A is not referenced.
Before entry with UPLO = 'L' or 'l', the leading k by k
lower triangular part of the array A must contain the lower
triangular matrix and the strictly upper triangular part of
A is not referenced.
Note that when DIAG = 'U' or 'u', the diagonal elements of
A are not referenced either, but are assumed to be unity.
LDA On entry, LDA specifies the first dimension of A as declared in the calling (sub) program. When SIDE = ‘L’ or ‘l’ then LDA must be at least max( 1, m ), when SIDE = ‘R’ or ‘r’ then LDA must be at least max( 1, n ).
B Before entry, the leading m by n part of the array B must contain the right-hand side matrix B, and on exit is overwritten by the solution matrix X.
LDB On entry, LDB specifies the first dimension of B as declared in the calling (sub) program. LDB must be at least max( 1, m ).
date:December 2016
FURTHER DETAILS
Level 3 Blas routine.
– Written on 8-February-1989. Jack Dongarra, Argonne National Laboratory. Iain Duff, AERE Harwell. Jeremy Du Croz, Numerical Algorithms Group Ltd. Sven Hammarling, Numerical Algorithms Group Ltd.
Online html documentation available at
http://www.netlib.org/lapack/explore-html/
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| character(len=1), | intent(in) | :: | side | |||
| character(len=1), | intent(in) | :: | uplo | |||
| character(len=1), | intent(in) | :: | transa | |||
| character(len=1), | intent(in) | :: | diag | |||
| integer, | intent(in) | :: | m | |||
| integer, | intent(in) | :: | n | |||
| double precision, | intent(in) | :: | alpha | |||
| double precision, | intent(in) | :: | a(lda,*) | |||
| integer, | intent(in) | :: | lda | |||
| double precision, | intent(inout) | :: | b(ldb,*) | |||
| integer, | intent(in) | :: | ldb |
| Type | Visibility | Attributes | Name | Initial | |||
|---|---|---|---|---|---|---|---|
| integer, | public | :: | i | ||||
| integer, | public | :: | info | ||||
| integer, | public | :: | j | ||||
| integer, | public | :: | k | ||||
| logical, | public | :: | lside | ||||
| logical, | public | :: | nounit | ||||
| integer, | public | :: | nrowa | ||||
| double precision, | public, | parameter | :: | one | = | 1.0d+0 | |
| double precision, | public | :: | temp | ||||
| logical, | public | :: | upper | ||||
| double precision, | public, | parameter | :: | zero | = | 0.0d+0 |
subroutine dtrsm(side,uplo,transa,diag,m,n,alpha,a,lda,b,ldb)
implicit none
!
! -- Reference BLAS level3 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 ..
double precision,intent(in) :: alpha
integer,intent(in) :: lda,ldb,m,n
character,intent(in) :: diag,side,transa,uplo
! ..
! .. Array Arguments ..
double precision,intent(in) :: a(lda,*)
double precision,intent(inout) :: b(ldb,*)
! ..
!
! =====================================================================
!
! .. External Functions ..
! ..
! .. External Subroutines ..
! ..
! .. Intrinsic Functions ..
intrinsic max
! ..
! .. Local Scalars ..
double precision temp
integer i,info,j,k,nrowa
logical lside,nounit,upper
! ..
! .. Parameters ..
double precision one,zero
parameter (one=1.0d+0,zero=0.0d+0)
! ..
!
! Test the input parameters.
!
lside = lsame(side,'L')
if (lside) then
nrowa = m
else
nrowa = n
endif
nounit = lsame(diag,'N')
upper = lsame(uplo,'U')
!
info = 0
if ((.not.lside) .and. (.not.lsame(side,'R'))) then
info = 1
elseif ((.not.upper) .and. (.not.lsame(uplo,'L'))) then
info = 2
elseif ((.not.lsame(transa,'N')) .and. (.not.lsame(transa,'T')) .and. (.not.lsame(transa,'C'))) then
info = 3
elseif ((.not.lsame(diag,'U')) .and. (.not.lsame(diag,'N'))) then
info = 4
elseif (m.lt.0) then
info = 5
elseif (n.lt.0) then
info = 6
elseif (lda.lt.max(1,nrowa)) then
info = 9
elseif (ldb.lt.max(1,m)) then
info = 11
endif
if (info.ne.0) then
call xerbla('DTRSM ',info)
return
endif
!
! Quick return if possible.
!
if (m.eq.0 .or. n.eq.0) return
!
! And when alpha.eq.zero.
!
if (alpha.eq.zero) then
b(1:m,1:n) = zero
return
endif
!
! Start the operations.
!
if (lside) then
if (lsame(transa,'N')) then
!
! Form B := alpha*inv( A )*B.
!
if (upper) then
do j = 1,n
if (alpha.ne.one) then
b(1:m,j) = alpha*b(1:m,j)
endif
do k = m,1,-1
if (b(k,j).ne.zero) then
if (nounit) b(k,j) = b(k,j)/a(k,k)
do i = 1,k - 1
b(i,j) = b(i,j) - b(k,j)*a(i,k)
enddo
endif
enddo
enddo
else
do j = 1,n
if (alpha.ne.one) then
b(1:m,j) = alpha*b(1:m,j)
endif
do k = 1,m
if (b(k,j).ne.zero) then
if (nounit) b(k,j) = b(k,j)/a(k,k)
do i = k + 1,m
b(i,j) = b(i,j) - b(k,j)*a(i,k)
enddo
endif
enddo
enddo
endif
else
!
! Form B := alpha*inv( A**T )*B.
!
if (upper) then
do j = 1,n
do i = 1,m
temp = alpha*b(i,j)
do k = 1,i - 1
temp = temp - a(k,i)*b(k,j)
enddo
if (nounit) temp = temp/a(i,i)
b(i,j) = temp
enddo
enddo
else
do j = 1,n
do i = m,1,-1
temp = alpha*b(i,j)
do k = i + 1,m
temp = temp - a(k,i)*b(k,j)
enddo
if (nounit) temp = temp/a(i,i)
b(i,j) = temp
enddo
enddo
endif
endif
else
if (lsame(transa,'N')) then
!
! Form B := alpha*B*inv( A ).
!
if (upper) then
do j = 1,n
if (alpha.ne.one) then
b(1:m,j) = alpha*b(1:m,j)
endif
do k = 1,j - 1
if (a(k,j).ne.zero) then
b(1:m,j) = b(1:m,j) - a(k,j)*b(1:m,k)
endif
enddo
if (nounit) then
temp = one/a(j,j)
b(1:m,j) = temp*b(1:m,j)
endif
enddo
else
do j = n,1,-1
if (alpha.ne.one) then
b(1:m,j) = alpha*b(1:m,j)
endif
do k = j + 1,n
if (a(k,j).ne.zero) then
b(1:m,j) = b(1:m,j) - a(k,j)*b(1:m,k)
endif
enddo
if (nounit) then
temp = one/a(j,j)
b(1:m,j) = temp*b(1:m,j)
endif
enddo
endif
else
!
! Form B := alpha*B*inv( A**T ).
!
if (upper) then
do k = n,1,-1
if (nounit) then
temp = one/a(k,k)
b(1:m,k) = temp*b(1:m,k)
endif
do j = 1,k - 1
if (a(j,k).ne.zero) then
temp = a(j,k)
b(1:m,j) = b(1:m,j) - temp*b(1:m,k)
endif
enddo
if (alpha.ne.one) then
b(1:m,k) = alpha*b(1:m,k)
endif
enddo
else
do k = 1,n
if (nounit) then
temp = one/a(k,k)
b(1:m,k) = temp*b(1:m,k)
endif
do j = k + 1,n
if (a(j,k).ne.zero) then
temp = a(j,k)
b(1:m,j) = b(1:m,j) - temp*b(1:m,k)
endif
enddo
if (alpha.ne.one) then
b(1:m,k) = alpha*b(1:m,k)
endif
enddo
endif
endif
endif
!
end subroutine dtrsm