#:include "common.fypp"
submodule(stdlib_lapack_solve) stdlib_lapack_solve_tri_comp
implicit none
contains
#:for ik,it,ii in LINALG_INT_KINDS_TYPES
module subroutine stdlib${ii}$_strcon( norm, uplo, diag, n, a, lda, rcond, work,iwork, info )
!! STRCON estimates the reciprocal of the condition number of a
!! triangular matrix A, in either the 1-norm or the infinity-norm.
!! The norm of A is computed and an estimate is obtained for
!! norm(inv(A)), then the reciprocal of the condition number is
!! computed as
!! RCOND = 1 / ( norm(A) * norm(inv(A)) ).
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: diag, norm, uplo
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: lda, n
real(sp), intent(out) :: rcond
! Array Arguments
integer(${ik}$), intent(out) :: iwork(*)
real(sp), intent(in) :: a(lda,*)
real(sp), intent(out) :: work(*)
! =====================================================================
! Local Scalars
logical(lk) :: nounit, onenrm, upper
character :: normin
integer(${ik}$) :: ix, kase, kase1
real(sp) :: ainvnm, anorm, scale, smlnum, xnorm
! Local Arrays
integer(${ik}$) :: isave(3_${ik}$)
! Intrinsic Functions
! Executable Statements
! test the input parameters.
info = 0_${ik}$
upper = stdlib_lsame( uplo, 'U' )
onenrm = norm=='1' .or. stdlib_lsame( norm, 'O' )
nounit = stdlib_lsame( diag, 'N' )
if( .not.onenrm .and. .not.stdlib_lsame( norm, 'I' ) ) then
info = -1_${ik}$
else if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then
info = -2_${ik}$
else if( .not.nounit .and. .not.stdlib_lsame( diag, 'U' ) ) then
info = -3_${ik}$
else if( n<0_${ik}$ ) then
info = -4_${ik}$
else if( lda<max( 1_${ik}$, n ) ) then
info = -6_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'STRCON', -info )
return
end if
! quick return if possible
if( n==0_${ik}$ ) then
rcond = one
return
end if
rcond = zero
smlnum = stdlib${ii}$_slamch( 'SAFE MINIMUM' )*real( max( 1_${ik}$, n ),KIND=sp)
! compute the norm of the triangular matrix a.
anorm = stdlib${ii}$_slantr( norm, uplo, diag, n, n, a, lda, work )
! continue only if anorm > 0.
if( anorm>zero ) then
! estimate the norm of the inverse of a.
ainvnm = zero
normin = 'N'
if( onenrm ) then
kase1 = 1_${ik}$
else
kase1 = 2_${ik}$
end if
kase = 0_${ik}$
10 continue
call stdlib${ii}$_slacn2( n, work( n+1 ), work, iwork, ainvnm, kase, isave )
if( kase/=0_${ik}$ ) then
if( kase==kase1 ) then
! multiply by inv(a).
call stdlib${ii}$_slatrs( uplo, 'NO TRANSPOSE', diag, normin, n, a,lda, work, scale,&
work( 2_${ik}$*n+1 ), info )
else
! multiply by inv(a**t).
call stdlib${ii}$_slatrs( uplo, 'TRANSPOSE', diag, normin, n, a, lda,work, scale, &
work( 2_${ik}$*n+1 ), info )
end if
normin = 'Y'
! multiply by 1/scale if doing so will not cause overflow.
if( scale/=one ) then
ix = stdlib${ii}$_isamax( n, work, 1_${ik}$ )
xnorm = abs( work( ix ) )
if( scale<xnorm*smlnum .or. scale==zero )go to 20
call stdlib${ii}$_srscl( n, scale, work, 1_${ik}$ )
end if
go to 10
end if
! compute the estimate of the reciprocal condition number.
if( ainvnm/=zero )rcond = ( one / anorm ) / ainvnm
end if
20 continue
return
end subroutine stdlib${ii}$_strcon
module subroutine stdlib${ii}$_dtrcon( norm, uplo, diag, n, a, lda, rcond, work,iwork, info )
!! DTRCON estimates the reciprocal of the condition number of a
!! triangular matrix A, in either the 1-norm or the infinity-norm.
!! The norm of A is computed and an estimate is obtained for
!! norm(inv(A)), then the reciprocal of the condition number is
!! computed as
!! RCOND = 1 / ( norm(A) * norm(inv(A)) ).
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: diag, norm, uplo
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: lda, n
real(dp), intent(out) :: rcond
! Array Arguments
integer(${ik}$), intent(out) :: iwork(*)
real(dp), intent(in) :: a(lda,*)
real(dp), intent(out) :: work(*)
! =====================================================================
! Local Scalars
logical(lk) :: nounit, onenrm, upper
character :: normin
integer(${ik}$) :: ix, kase, kase1
real(dp) :: ainvnm, anorm, scale, smlnum, xnorm
! Local Arrays
integer(${ik}$) :: isave(3_${ik}$)
! Intrinsic Functions
! Executable Statements
! test the input parameters.
info = 0_${ik}$
upper = stdlib_lsame( uplo, 'U' )
onenrm = norm=='1' .or. stdlib_lsame( norm, 'O' )
nounit = stdlib_lsame( diag, 'N' )
if( .not.onenrm .and. .not.stdlib_lsame( norm, 'I' ) ) then
info = -1_${ik}$
else if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then
info = -2_${ik}$
else if( .not.nounit .and. .not.stdlib_lsame( diag, 'U' ) ) then
info = -3_${ik}$
else if( n<0_${ik}$ ) then
info = -4_${ik}$
else if( lda<max( 1_${ik}$, n ) ) then
info = -6_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'DTRCON', -info )
return
end if
! quick return if possible
if( n==0_${ik}$ ) then
rcond = one
return
end if
rcond = zero
smlnum = stdlib${ii}$_dlamch( 'SAFE MINIMUM' )*real( max( 1_${ik}$, n ),KIND=dp)
! compute the norm of the triangular matrix a.
anorm = stdlib${ii}$_dlantr( norm, uplo, diag, n, n, a, lda, work )
! continue only if anorm > 0.
if( anorm>zero ) then
! estimate the norm of the inverse of a.
ainvnm = zero
normin = 'N'
if( onenrm ) then
kase1 = 1_${ik}$
else
kase1 = 2_${ik}$
end if
kase = 0_${ik}$
10 continue
call stdlib${ii}$_dlacn2( n, work( n+1 ), work, iwork, ainvnm, kase, isave )
if( kase/=0_${ik}$ ) then
if( kase==kase1 ) then
! multiply by inv(a).
call stdlib${ii}$_dlatrs( uplo, 'NO TRANSPOSE', diag, normin, n, a,lda, work, scale,&
work( 2_${ik}$*n+1 ), info )
else
! multiply by inv(a**t).
call stdlib${ii}$_dlatrs( uplo, 'TRANSPOSE', diag, normin, n, a, lda,work, scale, &
work( 2_${ik}$*n+1 ), info )
end if
normin = 'Y'
! multiply by 1/scale if doing so will not cause overflow.
if( scale/=one ) then
ix = stdlib${ii}$_idamax( n, work, 1_${ik}$ )
xnorm = abs( work( ix ) )
if( scale<xnorm*smlnum .or. scale==zero )go to 20
call stdlib${ii}$_drscl( n, scale, work, 1_${ik}$ )
end if
go to 10
end if
! compute the estimate of the reciprocal condition number.
if( ainvnm/=zero )rcond = ( one / anorm ) / ainvnm
end if
20 continue
return
end subroutine stdlib${ii}$_dtrcon
#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
module subroutine stdlib${ii}$_${ri}$trcon( norm, uplo, diag, n, a, lda, rcond, work,iwork, info )
!! DTRCON: estimates the reciprocal of the condition number of a
!! triangular matrix A, in either the 1-norm or the infinity-norm.
!! The norm of A is computed and an estimate is obtained for
!! norm(inv(A)), then the reciprocal of the condition number is
!! computed as
!! RCOND = 1 / ( norm(A) * norm(inv(A)) ).
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_${rk}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: diag, norm, uplo
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: lda, n
real(${rk}$), intent(out) :: rcond
! Array Arguments
integer(${ik}$), intent(out) :: iwork(*)
real(${rk}$), intent(in) :: a(lda,*)
real(${rk}$), intent(out) :: work(*)
! =====================================================================
! Local Scalars
logical(lk) :: nounit, onenrm, upper
character :: normin
integer(${ik}$) :: ix, kase, kase1
real(${rk}$) :: ainvnm, anorm, scale, smlnum, xnorm
! Local Arrays
integer(${ik}$) :: isave(3_${ik}$)
! Intrinsic Functions
! Executable Statements
! test the input parameters.
info = 0_${ik}$
upper = stdlib_lsame( uplo, 'U' )
onenrm = norm=='1' .or. stdlib_lsame( norm, 'O' )
nounit = stdlib_lsame( diag, 'N' )
if( .not.onenrm .and. .not.stdlib_lsame( norm, 'I' ) ) then
info = -1_${ik}$
else if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then
info = -2_${ik}$
else if( .not.nounit .and. .not.stdlib_lsame( diag, 'U' ) ) then
info = -3_${ik}$
else if( n<0_${ik}$ ) then
info = -4_${ik}$
else if( lda<max( 1_${ik}$, n ) ) then
info = -6_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'DTRCON', -info )
return
end if
! quick return if possible
if( n==0_${ik}$ ) then
rcond = one
return
end if
rcond = zero
smlnum = stdlib${ii}$_${ri}$lamch( 'SAFE MINIMUM' )*real( max( 1_${ik}$, n ),KIND=${rk}$)
! compute the norm of the triangular matrix a.
anorm = stdlib${ii}$_${ri}$lantr( norm, uplo, diag, n, n, a, lda, work )
! continue only if anorm > 0.
if( anorm>zero ) then
! estimate the norm of the inverse of a.
ainvnm = zero
normin = 'N'
if( onenrm ) then
kase1 = 1_${ik}$
else
kase1 = 2_${ik}$
end if
kase = 0_${ik}$
10 continue
call stdlib${ii}$_${ri}$lacn2( n, work( n+1 ), work, iwork, ainvnm, kase, isave )
if( kase/=0_${ik}$ ) then
if( kase==kase1 ) then
! multiply by inv(a).
call stdlib${ii}$_${ri}$latrs( uplo, 'NO TRANSPOSE', diag, normin, n, a,lda, work, scale,&
work( 2_${ik}$*n+1 ), info )
else
! multiply by inv(a**t).
call stdlib${ii}$_${ri}$latrs( uplo, 'TRANSPOSE', diag, normin, n, a, lda,work, scale, &
work( 2_${ik}$*n+1 ), info )
end if
normin = 'Y'
! multiply by 1/scale if doing so will not cause overflow.
if( scale/=one ) then
ix = stdlib${ii}$_i${ri}$amax( n, work, 1_${ik}$ )
xnorm = abs( work( ix ) )
if( scale<xnorm*smlnum .or. scale==zero )go to 20
call stdlib${ii}$_${ri}$rscl( n, scale, work, 1_${ik}$ )
end if
go to 10
end if
! compute the estimate of the reciprocal condition number.
if( ainvnm/=zero )rcond = ( one / anorm ) / ainvnm
end if
20 continue
return
end subroutine stdlib${ii}$_${ri}$trcon
#:endif
#:endfor
module subroutine stdlib${ii}$_ctrcon( norm, uplo, diag, n, a, lda, rcond, work,rwork, info )
!! CTRCON estimates the reciprocal of the condition number of a
!! triangular matrix A, in either the 1-norm or the infinity-norm.
!! The norm of A is computed and an estimate is obtained for
!! norm(inv(A)), then the reciprocal of the condition number is
!! computed as
!! RCOND = 1 / ( norm(A) * norm(inv(A)) ).
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: diag, norm, uplo
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: lda, n
real(sp), intent(out) :: rcond
! Array Arguments
real(sp), intent(out) :: rwork(*)
complex(sp), intent(in) :: a(lda,*)
complex(sp), intent(out) :: work(*)
! =====================================================================
! Local Scalars
logical(lk) :: nounit, onenrm, upper
character :: normin
integer(${ik}$) :: ix, kase, kase1
real(sp) :: ainvnm, anorm, scale, smlnum, xnorm
complex(sp) :: zdum
! Local Arrays
integer(${ik}$) :: isave(3_${ik}$)
! Intrinsic Functions
! Statement Functions
real(sp) :: cabs1
! Statement Function Definitions
cabs1( zdum ) = abs( real( zdum,KIND=sp) ) + abs( aimag( zdum ) )
! Executable Statements
! test the input parameters.
info = 0_${ik}$
upper = stdlib_lsame( uplo, 'U' )
onenrm = norm=='1' .or. stdlib_lsame( norm, 'O' )
nounit = stdlib_lsame( diag, 'N' )
if( .not.onenrm .and. .not.stdlib_lsame( norm, 'I' ) ) then
info = -1_${ik}$
else if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then
info = -2_${ik}$
else if( .not.nounit .and. .not.stdlib_lsame( diag, 'U' ) ) then
info = -3_${ik}$
else if( n<0_${ik}$ ) then
info = -4_${ik}$
else if( lda<max( 1_${ik}$, n ) ) then
info = -6_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'CTRCON', -info )
return
end if
! quick return if possible
if( n==0_${ik}$ ) then
rcond = one
return
end if
rcond = zero
smlnum = stdlib${ii}$_slamch( 'SAFE MINIMUM' )*real( max( 1_${ik}$, n ),KIND=sp)
! compute the norm of the triangular matrix a.
anorm = stdlib${ii}$_clantr( norm, uplo, diag, n, n, a, lda, rwork )
! continue only if anorm > 0.
if( anorm>zero ) then
! estimate the norm of the inverse of a.
ainvnm = zero
normin = 'N'
if( onenrm ) then
kase1 = 1_${ik}$
else
kase1 = 2_${ik}$
end if
kase = 0_${ik}$
10 continue
call stdlib${ii}$_clacn2( n, work( n+1 ), work, ainvnm, kase, isave )
if( kase/=0_${ik}$ ) then
if( kase==kase1 ) then
! multiply by inv(a).
call stdlib${ii}$_clatrs( uplo, 'NO TRANSPOSE', diag, normin, n, a,lda, work, scale,&
rwork, info )
else
! multiply by inv(a**h).
call stdlib${ii}$_clatrs( uplo, 'CONJUGATE TRANSPOSE', diag, normin,n, a, lda, work,&
scale, rwork, info )
end if
normin = 'Y'
! multiply by 1/scale if doing so will not cause overflow.
if( scale/=one ) then
ix = stdlib${ii}$_icamax( n, work, 1_${ik}$ )
xnorm = cabs1( work( ix ) )
if( scale<xnorm*smlnum .or. scale==zero )go to 20
call stdlib${ii}$_csrscl( n, scale, work, 1_${ik}$ )
end if
go to 10
end if
! compute the estimate of the reciprocal condition number.
if( ainvnm/=zero )rcond = ( one / anorm ) / ainvnm
end if
20 continue
return
end subroutine stdlib${ii}$_ctrcon
module subroutine stdlib${ii}$_ztrcon( norm, uplo, diag, n, a, lda, rcond, work,rwork, info )
!! ZTRCON estimates the reciprocal of the condition number of a
!! triangular matrix A, in either the 1-norm or the infinity-norm.
!! The norm of A is computed and an estimate is obtained for
!! norm(inv(A)), then the reciprocal of the condition number is
!! computed as
!! RCOND = 1 / ( norm(A) * norm(inv(A)) ).
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: diag, norm, uplo
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: lda, n
real(dp), intent(out) :: rcond
! Array Arguments
real(dp), intent(out) :: rwork(*)
complex(dp), intent(in) :: a(lda,*)
complex(dp), intent(out) :: work(*)
! =====================================================================
! Local Scalars
logical(lk) :: nounit, onenrm, upper
character :: normin
integer(${ik}$) :: ix, kase, kase1
real(dp) :: ainvnm, anorm, scale, smlnum, xnorm
complex(dp) :: zdum
! Local Arrays
integer(${ik}$) :: isave(3_${ik}$)
! Intrinsic Functions
! Statement Functions
real(dp) :: cabs1
! Statement Function Definitions
cabs1( zdum ) = abs( real( zdum,KIND=dp) ) + abs( aimag( zdum ) )
! Executable Statements
! test the input parameters.
info = 0_${ik}$
upper = stdlib_lsame( uplo, 'U' )
onenrm = norm=='1' .or. stdlib_lsame( norm, 'O' )
nounit = stdlib_lsame( diag, 'N' )
if( .not.onenrm .and. .not.stdlib_lsame( norm, 'I' ) ) then
info = -1_${ik}$
else if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then
info = -2_${ik}$
else if( .not.nounit .and. .not.stdlib_lsame( diag, 'U' ) ) then
info = -3_${ik}$
else if( n<0_${ik}$ ) then
info = -4_${ik}$
else if( lda<max( 1_${ik}$, n ) ) then
info = -6_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'ZTRCON', -info )
return
end if
! quick return if possible
if( n==0_${ik}$ ) then
rcond = one
return
end if
rcond = zero
smlnum = stdlib${ii}$_dlamch( 'SAFE MINIMUM' )*real( max( 1_${ik}$, n ),KIND=dp)
! compute the norm of the triangular matrix a.
anorm = stdlib${ii}$_zlantr( norm, uplo, diag, n, n, a, lda, rwork )
! continue only if anorm > 0.
if( anorm>zero ) then
! estimate the norm of the inverse of a.
ainvnm = zero
normin = 'N'
if( onenrm ) then
kase1 = 1_${ik}$
else
kase1 = 2_${ik}$
end if
kase = 0_${ik}$
10 continue
call stdlib${ii}$_zlacn2( n, work( n+1 ), work, ainvnm, kase, isave )
if( kase/=0_${ik}$ ) then
if( kase==kase1 ) then
! multiply by inv(a).
call stdlib${ii}$_zlatrs( uplo, 'NO TRANSPOSE', diag, normin, n, a,lda, work, scale,&
rwork, info )
else
! multiply by inv(a**h).
call stdlib${ii}$_zlatrs( uplo, 'CONJUGATE TRANSPOSE', diag, normin,n, a, lda, work,&
scale, rwork, info )
end if
normin = 'Y'
! multiply by 1/scale if doing so will not cause overflow.
if( scale/=one ) then
ix = stdlib${ii}$_izamax( n, work, 1_${ik}$ )
xnorm = cabs1( work( ix ) )
if( scale<xnorm*smlnum .or. scale==zero )go to 20
call stdlib${ii}$_zdrscl( n, scale, work, 1_${ik}$ )
end if
go to 10
end if
! compute the estimate of the reciprocal condition number.
if( ainvnm/=zero )rcond = ( one / anorm ) / ainvnm
end if
20 continue
return
end subroutine stdlib${ii}$_ztrcon
#:for ck,ct,ci in CMPLX_KINDS_TYPES
#:if not ck in ["sp","dp"]
module subroutine stdlib${ii}$_${ci}$trcon( norm, uplo, diag, n, a, lda, rcond, work,rwork, info )
!! ZTRCON: estimates the reciprocal of the condition number of a
!! triangular matrix A, in either the 1-norm or the infinity-norm.
!! The norm of A is computed and an estimate is obtained for
!! norm(inv(A)), then the reciprocal of the condition number is
!! computed as
!! RCOND = 1 / ( norm(A) * norm(inv(A)) ).
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_${ck}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: diag, norm, uplo
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: lda, n
real(${ck}$), intent(out) :: rcond
! Array Arguments
real(${ck}$), intent(out) :: rwork(*)
complex(${ck}$), intent(in) :: a(lda,*)
complex(${ck}$), intent(out) :: work(*)
! =====================================================================
! Local Scalars
logical(lk) :: nounit, onenrm, upper
character :: normin
integer(${ik}$) :: ix, kase, kase1
real(${ck}$) :: ainvnm, anorm, scale, smlnum, xnorm
complex(${ck}$) :: zdum
! Local Arrays
integer(${ik}$) :: isave(3_${ik}$)
! Intrinsic Functions
! Statement Functions
real(${ck}$) :: cabs1
! Statement Function Definitions
cabs1( zdum ) = abs( real( zdum,KIND=${ck}$) ) + abs( aimag( zdum ) )
! Executable Statements
! test the input parameters.
info = 0_${ik}$
upper = stdlib_lsame( uplo, 'U' )
onenrm = norm=='1' .or. stdlib_lsame( norm, 'O' )
nounit = stdlib_lsame( diag, 'N' )
if( .not.onenrm .and. .not.stdlib_lsame( norm, 'I' ) ) then
info = -1_${ik}$
else if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then
info = -2_${ik}$
else if( .not.nounit .and. .not.stdlib_lsame( diag, 'U' ) ) then
info = -3_${ik}$
else if( n<0_${ik}$ ) then
info = -4_${ik}$
else if( lda<max( 1_${ik}$, n ) ) then
info = -6_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'ZTRCON', -info )
return
end if
! quick return if possible
if( n==0_${ik}$ ) then
rcond = one
return
end if
rcond = zero
smlnum = stdlib${ii}$_${c2ri(ci)}$lamch( 'SAFE MINIMUM' )*real( max( 1_${ik}$, n ),KIND=${ck}$)
! compute the norm of the triangular matrix a.
anorm = stdlib${ii}$_${ci}$lantr( norm, uplo, diag, n, n, a, lda, rwork )
! continue only if anorm > 0.
if( anorm>zero ) then
! estimate the norm of the inverse of a.
ainvnm = zero
normin = 'N'
if( onenrm ) then
kase1 = 1_${ik}$
else
kase1 = 2_${ik}$
end if
kase = 0_${ik}$
10 continue
call stdlib${ii}$_${ci}$lacn2( n, work( n+1 ), work, ainvnm, kase, isave )
if( kase/=0_${ik}$ ) then
if( kase==kase1 ) then
! multiply by inv(a).
call stdlib${ii}$_${ci}$latrs( uplo, 'NO TRANSPOSE', diag, normin, n, a,lda, work, scale,&
rwork, info )
else
! multiply by inv(a**h).
call stdlib${ii}$_${ci}$latrs( uplo, 'CONJUGATE TRANSPOSE', diag, normin,n, a, lda, work,&
scale, rwork, info )
end if
normin = 'Y'
! multiply by 1/scale if doing so will not cause overflow.
if( scale/=one ) then
ix = stdlib${ii}$_i${ci}$amax( n, work, 1_${ik}$ )
xnorm = cabs1( work( ix ) )
if( scale<xnorm*smlnum .or. scale==zero )go to 20
call stdlib${ii}$_${ci}$drscl( n, scale, work, 1_${ik}$ )
end if
go to 10
end if
! compute the estimate of the reciprocal condition number.
if( ainvnm/=zero )rcond = ( one / anorm ) / ainvnm
end if
20 continue
return
end subroutine stdlib${ii}$_${ci}$trcon
#:endif
#:endfor
pure module subroutine stdlib${ii}$_strtrs( uplo, trans, diag, n, nrhs, a, lda, b, ldb,info )
!! STRTRS solves a triangular system of the form
!! A * X = B or A**T * X = B,
!! where A is a triangular matrix of order N, and B is an N-by-NRHS
!! matrix. A check is made to verify that A is nonsingular.
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: diag, trans, uplo
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: lda, ldb, n, nrhs
! Array Arguments
real(sp), intent(in) :: a(lda,*)
real(sp), intent(inout) :: b(ldb,*)
! =====================================================================
! Local Scalars
logical(lk) :: nounit
! Intrinsic Functions
! Executable Statements
! test the input parameters.
info = 0_${ik}$
nounit = stdlib_lsame( diag, 'N' )
if( .not.stdlib_lsame( uplo, 'U' ) .and. .not.stdlib_lsame( uplo, 'L' ) ) then
info = -1_${ik}$
else if( .not.stdlib_lsame( trans, 'N' ) .and. .not.stdlib_lsame( trans, 'T' ) .and. &
.not.stdlib_lsame( trans, 'C' ) ) then
info = -2_${ik}$
else if( .not.nounit .and. .not.stdlib_lsame( diag, 'U' ) ) then
info = -3_${ik}$
else if( n<0_${ik}$ ) then
info = -4_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -5_${ik}$
else if( lda<max( 1_${ik}$, n ) ) then
info = -7_${ik}$
else if( ldb<max( 1_${ik}$, n ) ) then
info = -9_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'STRTRS', -info )
return
end if
! quick return if possible
if( n==0 )return
! check for singularity.
if( nounit ) then
do info = 1, n
if( a( info, info )==zero )return
end do
end if
info = 0_${ik}$
! solve a * x = b or a**t * x = b.
call stdlib${ii}$_strsm( 'LEFT', uplo, trans, diag, n, nrhs, one, a, lda, b,ldb )
return
end subroutine stdlib${ii}$_strtrs
pure module subroutine stdlib${ii}$_dtrtrs( uplo, trans, diag, n, nrhs, a, lda, b, ldb,info )
!! DTRTRS solves a triangular system of the form
!! A * X = B or A**T * X = B,
!! where A is a triangular matrix of order N, and B is an N-by-NRHS
!! matrix. A check is made to verify that A is nonsingular.
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: diag, trans, uplo
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: lda, ldb, n, nrhs
! Array Arguments
real(dp), intent(in) :: a(lda,*)
real(dp), intent(inout) :: b(ldb,*)
! =====================================================================
! Local Scalars
logical(lk) :: nounit
! Intrinsic Functions
! Executable Statements
! test the input parameters.
info = 0_${ik}$
nounit = stdlib_lsame( diag, 'N' )
if( .not.stdlib_lsame( uplo, 'U' ) .and. .not.stdlib_lsame( uplo, 'L' ) ) then
info = -1_${ik}$
else if( .not.stdlib_lsame( trans, 'N' ) .and. .not.stdlib_lsame( trans, 'T' ) .and. &
.not.stdlib_lsame( trans, 'C' ) ) then
info = -2_${ik}$
else if( .not.nounit .and. .not.stdlib_lsame( diag, 'U' ) ) then
info = -3_${ik}$
else if( n<0_${ik}$ ) then
info = -4_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -5_${ik}$
else if( lda<max( 1_${ik}$, n ) ) then
info = -7_${ik}$
else if( ldb<max( 1_${ik}$, n ) ) then
info = -9_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'DTRTRS', -info )
return
end if
! quick return if possible
if( n==0 )return
! check for singularity.
if( nounit ) then
do info = 1, n
if( a( info, info )==zero )return
end do
end if
info = 0_${ik}$
! solve a * x = b or a**t * x = b.
call stdlib${ii}$_dtrsm( 'LEFT', uplo, trans, diag, n, nrhs, one, a, lda, b,ldb )
return
end subroutine stdlib${ii}$_dtrtrs
#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
pure module subroutine stdlib${ii}$_${ri}$trtrs( uplo, trans, diag, n, nrhs, a, lda, b, ldb,info )
!! DTRTRS: solves a triangular system of the form
!! A * X = B or A**T * X = B,
!! where A is a triangular matrix of order N, and B is an N-by-NRHS
!! matrix. A check is made to verify that A is nonsingular.
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_${rk}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: diag, trans, uplo
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: lda, ldb, n, nrhs
! Array Arguments
real(${rk}$), intent(in) :: a(lda,*)
real(${rk}$), intent(inout) :: b(ldb,*)
! =====================================================================
! Local Scalars
logical(lk) :: nounit
! Intrinsic Functions
! Executable Statements
! test the input parameters.
info = 0_${ik}$
nounit = stdlib_lsame( diag, 'N' )
if( .not.stdlib_lsame( uplo, 'U' ) .and. .not.stdlib_lsame( uplo, 'L' ) ) then
info = -1_${ik}$
else if( .not.stdlib_lsame( trans, 'N' ) .and. .not.stdlib_lsame( trans, 'T' ) .and. &
.not.stdlib_lsame( trans, 'C' ) ) then
info = -2_${ik}$
else if( .not.nounit .and. .not.stdlib_lsame( diag, 'U' ) ) then
info = -3_${ik}$
else if( n<0_${ik}$ ) then
info = -4_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -5_${ik}$
else if( lda<max( 1_${ik}$, n ) ) then
info = -7_${ik}$
else if( ldb<max( 1_${ik}$, n ) ) then
info = -9_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'DTRTRS', -info )
return
end if
! quick return if possible
if( n==0 )return
! check for singularity.
if( nounit ) then
do info = 1, n
if( a( info, info )==zero )return
end do
end if
info = 0_${ik}$
! solve a * x = b or a**t * x = b.
call stdlib${ii}$_${ri}$trsm( 'LEFT', uplo, trans, diag, n, nrhs, one, a, lda, b,ldb )
return
end subroutine stdlib${ii}$_${ri}$trtrs
#:endif
#:endfor
pure module subroutine stdlib${ii}$_ctrtrs( uplo, trans, diag, n, nrhs, a, lda, b, ldb,info )
!! CTRTRS solves a triangular system of the form
!! A * X = B, A**T * X = B, or A**H * X = B,
!! where A is a triangular matrix of order N, and B is an N-by-NRHS
!! matrix. A check is made to verify that A is nonsingular.
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: diag, trans, uplo
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: lda, ldb, n, nrhs
! Array Arguments
complex(sp), intent(in) :: a(lda,*)
complex(sp), intent(inout) :: b(ldb,*)
! =====================================================================
! Local Scalars
logical(lk) :: nounit
! Intrinsic Functions
! Executable Statements
! test the input parameters.
info = 0_${ik}$
nounit = stdlib_lsame( diag, 'N' )
if( .not.stdlib_lsame( uplo, 'U' ) .and. .not.stdlib_lsame( uplo, 'L' ) ) then
info = -1_${ik}$
else if( .not.stdlib_lsame( trans, 'N' ) .and. .not.stdlib_lsame( trans, 'T' ) .and. &
.not.stdlib_lsame( trans, 'C' ) ) then
info = -2_${ik}$
else if( .not.nounit .and. .not.stdlib_lsame( diag, 'U' ) ) then
info = -3_${ik}$
else if( n<0_${ik}$ ) then
info = -4_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -5_${ik}$
else if( lda<max( 1_${ik}$, n ) ) then
info = -7_${ik}$
else if( ldb<max( 1_${ik}$, n ) ) then
info = -9_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'CTRTRS', -info )
return
end if
! quick return if possible
if( n==0 )return
! check for singularity.
if( nounit ) then
do info = 1, n
if( a( info, info )==czero )return
end do
end if
info = 0_${ik}$
! solve a * x = b, a**t * x = b, or a**h * x = b.
call stdlib${ii}$_ctrsm( 'LEFT', uplo, trans, diag, n, nrhs, cone, a, lda, b,ldb )
return
end subroutine stdlib${ii}$_ctrtrs
pure module subroutine stdlib${ii}$_ztrtrs( uplo, trans, diag, n, nrhs, a, lda, b, ldb,info )
!! ZTRTRS solves a triangular system of the form
!! A * X = B, A**T * X = B, or A**H * X = B,
!! where A is a triangular matrix of order N, and B is an N-by-NRHS
!! matrix. A check is made to verify that A is nonsingular.
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: diag, trans, uplo
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: lda, ldb, n, nrhs
! Array Arguments
complex(dp), intent(in) :: a(lda,*)
complex(dp), intent(inout) :: b(ldb,*)
! =====================================================================
! Local Scalars
logical(lk) :: nounit
! Intrinsic Functions
! Executable Statements
! test the input parameters.
info = 0_${ik}$
nounit = stdlib_lsame( diag, 'N' )
if( .not.stdlib_lsame( uplo, 'U' ) .and. .not.stdlib_lsame( uplo, 'L' ) ) then
info = -1_${ik}$
else if( .not.stdlib_lsame( trans, 'N' ) .and. .not.stdlib_lsame( trans, 'T' ) .and. &
.not.stdlib_lsame( trans, 'C' ) ) then
info = -2_${ik}$
else if( .not.nounit .and. .not.stdlib_lsame( diag, 'U' ) ) then
info = -3_${ik}$
else if( n<0_${ik}$ ) then
info = -4_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -5_${ik}$
else if( lda<max( 1_${ik}$, n ) ) then
info = -7_${ik}$
else if( ldb<max( 1_${ik}$, n ) ) then
info = -9_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'ZTRTRS', -info )
return
end if
! quick return if possible
if( n==0 )return
! check for singularity.
if( nounit ) then
do info = 1, n
if( a( info, info )==czero )return
end do
end if
info = 0_${ik}$
! solve a * x = b, a**t * x = b, or a**h * x = b.
call stdlib${ii}$_ztrsm( 'LEFT', uplo, trans, diag, n, nrhs, cone, a, lda, b,ldb )
return
end subroutine stdlib${ii}$_ztrtrs
#:for ck,ct,ci in CMPLX_KINDS_TYPES
#:if not ck in ["sp","dp"]
pure module subroutine stdlib${ii}$_${ci}$trtrs( uplo, trans, diag, n, nrhs, a, lda, b, ldb,info )
!! ZTRTRS: solves a triangular system of the form
!! A * X = B, A**T * X = B, or A**H * X = B,
!! where A is a triangular matrix of order N, and B is an N-by-NRHS
!! matrix. A check is made to verify that A is nonsingular.
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_${ck}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: diag, trans, uplo
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: lda, ldb, n, nrhs
! Array Arguments
complex(${ck}$), intent(in) :: a(lda,*)
complex(${ck}$), intent(inout) :: b(ldb,*)
! =====================================================================
! Local Scalars
logical(lk) :: nounit
! Intrinsic Functions
! Executable Statements
! test the input parameters.
info = 0_${ik}$
nounit = stdlib_lsame( diag, 'N' )
if( .not.stdlib_lsame( uplo, 'U' ) .and. .not.stdlib_lsame( uplo, 'L' ) ) then
info = -1_${ik}$
else if( .not.stdlib_lsame( trans, 'N' ) .and. .not.stdlib_lsame( trans, 'T' ) .and. &
.not.stdlib_lsame( trans, 'C' ) ) then
info = -2_${ik}$
else if( .not.nounit .and. .not.stdlib_lsame( diag, 'U' ) ) then
info = -3_${ik}$
else if( n<0_${ik}$ ) then
info = -4_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -5_${ik}$
else if( lda<max( 1_${ik}$, n ) ) then
info = -7_${ik}$
else if( ldb<max( 1_${ik}$, n ) ) then
info = -9_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'ZTRTRS', -info )
return
end if
! quick return if possible
if( n==0 )return
! check for singularity.
if( nounit ) then
do info = 1, n
if( a( info, info )==czero )return
end do
end if
info = 0_${ik}$
! solve a * x = b, a**t * x = b, or a**h * x = b.
call stdlib${ii}$_${ci}$trsm( 'LEFT', uplo, trans, diag, n, nrhs, cone, a, lda, b,ldb )
return
end subroutine stdlib${ii}$_${ci}$trtrs
#:endif
#:endfor
pure module subroutine stdlib${ii}$_slatrs( uplo, trans, diag, normin, n, a, lda, x, scale,cnorm, info )
!! SLATRS solves one of the triangular systems
!! A *x = s*b or A**T*x = s*b
!! with scaling to prevent overflow. Here A is an upper or lower
!! triangular matrix, A**T denotes the transpose of A, x and b are
!! n-element vectors, and s is a scaling factor, usually less than
!! or equal to 1, chosen so that the components of x will be less than
!! the overflow threshold. If the unscaled problem will not cause
!! overflow, the Level 2 BLAS routine STRSV is called. If the matrix A
!! is singular (A(j,j) = 0 for some j), then s is set to 0 and a
!! non-trivial solution to A*x = 0 is returned.
! -- lapack auxiliary routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: diag, normin, trans, uplo
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: lda, n
real(sp), intent(out) :: scale
! Array Arguments
real(sp), intent(in) :: a(lda,*)
real(sp), intent(inout) :: cnorm(*), x(*)
! =====================================================================
! Local Scalars
logical(lk) :: notran, nounit, upper
integer(${ik}$) :: i, imax, j, jfirst, jinc, jlast
real(sp) :: bignum, grow, rec, smlnum, sumj, tjj, tjjs, tmax, tscal, uscal, xbnd, xj, &
xmax
! Intrinsic Functions
! Executable Statements
info = 0_${ik}$
upper = stdlib_lsame( uplo, 'U' )
notran = stdlib_lsame( trans, 'N' )
nounit = stdlib_lsame( diag, 'N' )
! test the input parameters.
if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then
info = -1_${ik}$
else if( .not.notran .and. .not.stdlib_lsame( trans, 'T' ) .and. .not.stdlib_lsame( &
trans, 'C' ) ) then
info = -2_${ik}$
else if( .not.nounit .and. .not.stdlib_lsame( diag, 'U' ) ) then
info = -3_${ik}$
else if( .not.stdlib_lsame( normin, 'Y' ) .and. .not.stdlib_lsame( normin, 'N' ) ) &
then
info = -4_${ik}$
else if( n<0_${ik}$ ) then
info = -5_${ik}$
else if( lda<max( 1_${ik}$, n ) ) then
info = -7_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'SLATRS', -info )
return
end if
! quick return if possible
if( n==0 )return
! determine machine dependent parameters to control overflow.
smlnum = stdlib${ii}$_slamch( 'SAFE MINIMUM' ) / stdlib${ii}$_slamch( 'PRECISION' )
bignum = one / smlnum
scale = one
if( stdlib_lsame( normin, 'N' ) ) then
! compute the 1-norm of each column, not including the diagonal.
if( upper ) then
! a is upper triangular.
do j = 1, n
cnorm( j ) = stdlib${ii}$_sasum( j-1, a( 1_${ik}$, j ), 1_${ik}$ )
end do
else
! a is lower triangular.
do j = 1, n - 1
cnorm( j ) = stdlib${ii}$_sasum( n-j, a( j+1, j ), 1_${ik}$ )
end do
cnorm( n ) = zero
end if
end if
! scale the column norms by tscal if the maximum element in cnorm is
! greater than bignum.
imax = stdlib${ii}$_isamax( n, cnorm, 1_${ik}$ )
tmax = cnorm( imax )
if( tmax<=bignum ) then
tscal = one
else
tscal = one / ( smlnum*tmax )
call stdlib${ii}$_sscal( n, tscal, cnorm, 1_${ik}$ )
end if
! compute a bound on the computed solution vector to see if the
! level 2 blas routine stdlib${ii}$_strsv can be used.
j = stdlib${ii}$_isamax( n, x, 1_${ik}$ )
xmax = abs( x( j ) )
xbnd = xmax
if( notran ) then
! compute the growth in a * x = b.
if( upper ) then
jfirst = n
jlast = 1_${ik}$
jinc = -1_${ik}$
else
jfirst = 1_${ik}$
jlast = n
jinc = 1_${ik}$
end if
if( tscal/=one ) then
grow = zero
go to 50
end if
if( nounit ) then
! a is non-unit triangular.
! compute grow = 1/g(j) and xbnd = 1/m(j).
! initially, g(0) = max{x(i), i=1,...,n}.
grow = one / max( xbnd, smlnum )
xbnd = grow
do j = jfirst, jlast, jinc
! exit the loop if the growth factor is too small.
if( grow<=smlnum )go to 50
! m(j) = g(j-1) / abs(a(j,j))
tjj = abs( a( j, j ) )
xbnd = min( xbnd, min( one, tjj )*grow )
if( tjj+cnorm( j )>=smlnum ) then
! g(j) = g(j-1)*( 1 + cnorm(j) / abs(a(j,j)) )
grow = grow*( tjj / ( tjj+cnorm( j ) ) )
else
! g(j) could overflow, set grow to 0.
grow = zero
end if
end do
grow = xbnd
else
! a is unit triangular.
! compute grow = 1/g(j), where g(0) = max{x(i), i=1,...,n}.
grow = min( one, one / max( xbnd, smlnum ) )
do j = jfirst, jlast, jinc
! exit the loop if the growth factor is too small.
if( grow<=smlnum )go to 50
! g(j) = g(j-1)*( 1 + cnorm(j) )
grow = grow*( one / ( one+cnorm( j ) ) )
end do
end if
50 continue
else
! compute the growth in a**t * x = b.
if( upper ) then
jfirst = 1_${ik}$
jlast = n
jinc = 1_${ik}$
else
jfirst = n
jlast = 1_${ik}$
jinc = -1_${ik}$
end if
if( tscal/=one ) then
grow = zero
go to 80
end if
if( nounit ) then
! a is non-unit triangular.
! compute grow = 1/g(j) and xbnd = 1/m(j).
! initially, m(0) = max{x(i), i=1,...,n}.
grow = one / max( xbnd, smlnum )
xbnd = grow
do j = jfirst, jlast, jinc
! exit the loop if the growth factor is too small.
if( grow<=smlnum )go to 80
! g(j) = max( g(j-1), m(j-1)*( 1 + cnorm(j) ) )
xj = one + cnorm( j )
grow = min( grow, xbnd / xj )
! m(j) = m(j-1)*( 1 + cnorm(j) ) / abs(a(j,j))
tjj = abs( a( j, j ) )
if( xj>tjj )xbnd = xbnd*( tjj / xj )
end do
grow = min( grow, xbnd )
else
! a is unit triangular.
! compute grow = 1/g(j), where g(0) = max{x(i), i=1,...,n}.
grow = min( one, one / max( xbnd, smlnum ) )
do j = jfirst, jlast, jinc
! exit the loop if the growth factor is too small.
if( grow<=smlnum )go to 80
! g(j) = ( 1 + cnorm(j) )*g(j-1)
xj = one + cnorm( j )
grow = grow / xj
end do
end if
80 continue
end if
if( ( grow*tscal )>smlnum ) then
! use the level 2 blas solve if the reciprocal of the bound on
! elements of x is not too small.
call stdlib${ii}$_strsv( uplo, trans, diag, n, a, lda, x, 1_${ik}$ )
else
! use a level 1 blas solve, scaling intermediate results.
if( xmax>bignum ) then
! scale x so that its components are less than or equal to
! bignum in absolute value.
scale = bignum / xmax
call stdlib${ii}$_sscal( n, scale, x, 1_${ik}$ )
xmax = bignum
end if
if( notran ) then
! solve a * x = b
loop_100: do j = jfirst, jlast, jinc
! compute x(j) = b(j) / a(j,j), scaling x if necessary.
xj = abs( x( j ) )
if( nounit ) then
tjjs = a( j, j )*tscal
else
tjjs = tscal
if( tscal==one )go to 95
end if
tjj = abs( tjjs )
if( tjj>smlnum ) then
! abs(a(j,j)) > smlnum:
if( tjj<one ) then
if( xj>tjj*bignum ) then
! scale x by 1/b(j).
rec = one / xj
call stdlib${ii}$_sscal( n, rec, x, 1_${ik}$ )
scale = scale*rec
xmax = xmax*rec
end if
end if
x( j ) = x( j ) / tjjs
xj = abs( x( j ) )
else if( tjj>zero ) then
! 0 < abs(a(j,j)) <= smlnum:
if( xj>tjj*bignum ) then
! scale x by (1/abs(x(j)))*abs(a(j,j))*bignum
! to avoid overflow when dividing by a(j,j).
rec = ( tjj*bignum ) / xj
if( cnorm( j )>one ) then
! scale by 1/cnorm(j) to avoid overflow when
! multiplying x(j) times column j.
rec = rec / cnorm( j )
end if
call stdlib${ii}$_sscal( n, rec, x, 1_${ik}$ )
scale = scale*rec
xmax = xmax*rec
end if
x( j ) = x( j ) / tjjs
xj = abs( x( j ) )
else
! a(j,j) = 0: set x(1:n) = 0, x(j) = 1, and
! scale = 0, and compute a solution to a*x = 0.
do i = 1, n
x( i ) = zero
end do
x( j ) = one
xj = one
scale = zero
xmax = zero
end if
95 continue
! scale x if necessary to avoid overflow when adding a
! multiple of column j of a.
if( xj>one ) then
rec = one / xj
if( cnorm( j )>( bignum-xmax )*rec ) then
! scale x by 1/(2*abs(x(j))).
rec = rec*half
call stdlib${ii}$_sscal( n, rec, x, 1_${ik}$ )
scale = scale*rec
end if
else if( xj*cnorm( j )>( bignum-xmax ) ) then
! scale x by 1/2.
call stdlib${ii}$_sscal( n, half, x, 1_${ik}$ )
scale = scale*half
end if
if( upper ) then
if( j>1_${ik}$ ) then
! compute the update
! x(1:j-1) := x(1:j-1) - x(j) * a(1:j-1,j)
call stdlib${ii}$_saxpy( j-1, -x( j )*tscal, a( 1_${ik}$, j ), 1_${ik}$, x,1_${ik}$ )
i = stdlib${ii}$_isamax( j-1, x, 1_${ik}$ )
xmax = abs( x( i ) )
end if
else
if( j<n ) then
! compute the update
! x(j+1:n) := x(j+1:n) - x(j) * a(j+1:n,j)
call stdlib${ii}$_saxpy( n-j, -x( j )*tscal, a( j+1, j ), 1_${ik}$,x( j+1 ), 1_${ik}$ )
i = j + stdlib${ii}$_isamax( n-j, x( j+1 ), 1_${ik}$ )
xmax = abs( x( i ) )
end if
end if
end do loop_100
else
! solve a**t * x = b
loop_140: do j = jfirst, jlast, jinc
! compute x(j) = b(j) - sum a(k,j)*x(k).
! k<>j
xj = abs( x( j ) )
uscal = tscal
rec = one / max( xmax, one )
if( cnorm( j )>( bignum-xj )*rec ) then
! if x(j) could overflow, scale x by 1/(2*xmax).
rec = rec*half
if( nounit ) then
tjjs = a( j, j )*tscal
else
tjjs = tscal
end if
tjj = abs( tjjs )
if( tjj>one ) then
! divide by a(j,j) when scaling x if a(j,j) > 1.
rec = min( one, rec*tjj )
uscal = uscal / tjjs
end if
if( rec<one ) then
call stdlib${ii}$_sscal( n, rec, x, 1_${ik}$ )
scale = scale*rec
xmax = xmax*rec
end if
end if
sumj = zero
if( uscal==one ) then
! if the scaling needed for a in the dot product is 1,
! call stdlib${ii}$_sdot to perform the dot product.
if( upper ) then
sumj = stdlib${ii}$_sdot( j-1, a( 1_${ik}$, j ), 1_${ik}$, x, 1_${ik}$ )
else if( j<n ) then
sumj = stdlib${ii}$_sdot( n-j, a( j+1, j ), 1_${ik}$, x( j+1 ), 1_${ik}$ )
end if
else
! otherwise, use in-line code for the dot product.
if( upper ) then
do i = 1, j - 1
sumj = sumj + ( a( i, j )*uscal )*x( i )
end do
else if( j<n ) then
do i = j + 1, n
sumj = sumj + ( a( i, j )*uscal )*x( i )
end do
end if
end if
if( uscal==tscal ) then
! compute x(j) := ( x(j) - sumj ) / a(j,j) if 1/a(j,j)
! was not used to scale the dotproduct.
x( j ) = x( j ) - sumj
xj = abs( x( j ) )
if( nounit ) then
tjjs = a( j, j )*tscal
else
tjjs = tscal
if( tscal==one )go to 135
end if
! compute x(j) = x(j) / a(j,j), scaling if necessary.
tjj = abs( tjjs )
if( tjj>smlnum ) then
! abs(a(j,j)) > smlnum:
if( tjj<one ) then
if( xj>tjj*bignum ) then
! scale x by 1/abs(x(j)).
rec = one / xj
call stdlib${ii}$_sscal( n, rec, x, 1_${ik}$ )
scale = scale*rec
xmax = xmax*rec
end if
end if
x( j ) = x( j ) / tjjs
else if( tjj>zero ) then
! 0 < abs(a(j,j)) <= smlnum:
if( xj>tjj*bignum ) then
! scale x by (1/abs(x(j)))*abs(a(j,j))*bignum.
rec = ( tjj*bignum ) / xj
call stdlib${ii}$_sscal( n, rec, x, 1_${ik}$ )
scale = scale*rec
xmax = xmax*rec
end if
x( j ) = x( j ) / tjjs
else
! a(j,j) = 0: set x(1:n) = 0, x(j) = 1, and
! scale = 0, and compute a solution to a**t*x = 0.
do i = 1, n
x( i ) = zero
end do
x( j ) = one
scale = zero
xmax = zero
end if
135 continue
else
! compute x(j) := x(j) / a(j,j) - sumj if the dot
! product has already been divided by 1/a(j,j).
x( j ) = x( j ) / tjjs - sumj
end if
xmax = max( xmax, abs( x( j ) ) )
end do loop_140
end if
scale = scale / tscal
end if
! scale the column norms by 1/tscal for return.
if( tscal/=one ) then
call stdlib${ii}$_sscal( n, one / tscal, cnorm, 1_${ik}$ )
end if
return
end subroutine stdlib${ii}$_slatrs
pure module subroutine stdlib${ii}$_dlatrs( uplo, trans, diag, normin, n, a, lda, x, scale,cnorm, info )
!! DLATRS solves one of the triangular systems
!! A *x = s*b or A**T *x = s*b
!! with scaling to prevent overflow. Here A is an upper or lower
!! triangular matrix, A**T denotes the transpose of A, x and b are
!! n-element vectors, and s is a scaling factor, usually less than
!! or equal to 1, chosen so that the components of x will be less than
!! the overflow threshold. If the unscaled problem will not cause
!! overflow, the Level 2 BLAS routine DTRSV is called. If the matrix A
!! is singular (A(j,j) = 0 for some j), then s is set to 0 and a
!! non-trivial solution to A*x = 0 is returned.
! -- lapack auxiliary routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: diag, normin, trans, uplo
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: lda, n
real(dp), intent(out) :: scale
! Array Arguments
real(dp), intent(in) :: a(lda,*)
real(dp), intent(inout) :: cnorm(*), x(*)
! =====================================================================
! Local Scalars
logical(lk) :: notran, nounit, upper
integer(${ik}$) :: i, imax, j, jfirst, jinc, jlast
real(dp) :: bignum, grow, rec, smlnum, sumj, tjj, tjjs, tmax, tscal, uscal, xbnd, xj, &
xmax
! Intrinsic Functions
! Executable Statements
info = 0_${ik}$
upper = stdlib_lsame( uplo, 'U' )
notran = stdlib_lsame( trans, 'N' )
nounit = stdlib_lsame( diag, 'N' )
! test the input parameters.
if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then
info = -1_${ik}$
else if( .not.notran .and. .not.stdlib_lsame( trans, 'T' ) .and. .not.stdlib_lsame( &
trans, 'C' ) ) then
info = -2_${ik}$
else if( .not.nounit .and. .not.stdlib_lsame( diag, 'U' ) ) then
info = -3_${ik}$
else if( .not.stdlib_lsame( normin, 'Y' ) .and. .not.stdlib_lsame( normin, 'N' ) ) &
then
info = -4_${ik}$
else if( n<0_${ik}$ ) then
info = -5_${ik}$
else if( lda<max( 1_${ik}$, n ) ) then
info = -7_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'DLATRS', -info )
return
end if
! quick return if possible
if( n==0 )return
! determine machine dependent parameters to control overflow.
smlnum = stdlib${ii}$_dlamch( 'SAFE MINIMUM' ) / stdlib${ii}$_dlamch( 'PRECISION' )
bignum = one / smlnum
scale = one
if( stdlib_lsame( normin, 'N' ) ) then
! compute the 1-norm of each column, not including the diagonal.
if( upper ) then
! a is upper triangular.
do j = 1, n
cnorm( j ) = stdlib${ii}$_dasum( j-1, a( 1_${ik}$, j ), 1_${ik}$ )
end do
else
! a is lower triangular.
do j = 1, n - 1
cnorm( j ) = stdlib${ii}$_dasum( n-j, a( j+1, j ), 1_${ik}$ )
end do
cnorm( n ) = zero
end if
end if
! scale the column norms by tscal if the maximum element in cnorm is
! greater than bignum.
imax = stdlib${ii}$_idamax( n, cnorm, 1_${ik}$ )
tmax = cnorm( imax )
if( tmax<=bignum ) then
tscal = one
else
tscal = one / ( smlnum*tmax )
call stdlib${ii}$_dscal( n, tscal, cnorm, 1_${ik}$ )
end if
! compute a bound on the computed solution vector to see if the
! level 2 blas routine stdlib${ii}$_dtrsv can be used.
j = stdlib${ii}$_idamax( n, x, 1_${ik}$ )
xmax = abs( x( j ) )
xbnd = xmax
if( notran ) then
! compute the growth in a * x = b.
if( upper ) then
jfirst = n
jlast = 1_${ik}$
jinc = -1_${ik}$
else
jfirst = 1_${ik}$
jlast = n
jinc = 1_${ik}$
end if
if( tscal/=one ) then
grow = zero
go to 50
end if
if( nounit ) then
! a is non-unit triangular.
! compute grow = 1/g(j) and xbnd = 1/m(j).
! initially, g(0) = max{x(i), i=1,...,n}.
grow = one / max( xbnd, smlnum )
xbnd = grow
do j = jfirst, jlast, jinc
! exit the loop if the growth factor is too small.
if( grow<=smlnum )go to 50
! m(j) = g(j-1) / abs(a(j,j))
tjj = abs( a( j, j ) )
xbnd = min( xbnd, min( one, tjj )*grow )
if( tjj+cnorm( j )>=smlnum ) then
! g(j) = g(j-1)*( 1 + cnorm(j) / abs(a(j,j)) )
grow = grow*( tjj / ( tjj+cnorm( j ) ) )
else
! g(j) could overflow, set grow to 0.
grow = zero
end if
end do
grow = xbnd
else
! a is unit triangular.
! compute grow = 1/g(j), where g(0) = max{x(i), i=1,...,n}.
grow = min( one, one / max( xbnd, smlnum ) )
do j = jfirst, jlast, jinc
! exit the loop if the growth factor is too small.
if( grow<=smlnum )go to 50
! g(j) = g(j-1)*( 1 + cnorm(j) )
grow = grow*( one / ( one+cnorm( j ) ) )
end do
end if
50 continue
else
! compute the growth in a**t * x = b.
if( upper ) then
jfirst = 1_${ik}$
jlast = n
jinc = 1_${ik}$
else
jfirst = n
jlast = 1_${ik}$
jinc = -1_${ik}$
end if
if( tscal/=one ) then
grow = zero
go to 80
end if
if( nounit ) then
! a is non-unit triangular.
! compute grow = 1/g(j) and xbnd = 1/m(j).
! initially, m(0) = max{x(i), i=1,...,n}.
grow = one / max( xbnd, smlnum )
xbnd = grow
do j = jfirst, jlast, jinc
! exit the loop if the growth factor is too small.
if( grow<=smlnum )go to 80
! g(j) = max( g(j-1), m(j-1)*( 1 + cnorm(j) ) )
xj = one + cnorm( j )
grow = min( grow, xbnd / xj )
! m(j) = m(j-1)*( 1 + cnorm(j) ) / abs(a(j,j))
tjj = abs( a( j, j ) )
if( xj>tjj )xbnd = xbnd*( tjj / xj )
end do
grow = min( grow, xbnd )
else
! a is unit triangular.
! compute grow = 1/g(j), where g(0) = max{x(i), i=1,...,n}.
grow = min( one, one / max( xbnd, smlnum ) )
do j = jfirst, jlast, jinc
! exit the loop if the growth factor is too small.
if( grow<=smlnum )go to 80
! g(j) = ( 1 + cnorm(j) )*g(j-1)
xj = one + cnorm( j )
grow = grow / xj
end do
end if
80 continue
end if
if( ( grow*tscal )>smlnum ) then
! use the level 2 blas solve if the reciprocal of the bound on
! elements of x is not too small.
call stdlib${ii}$_dtrsv( uplo, trans, diag, n, a, lda, x, 1_${ik}$ )
else
! use a level 1 blas solve, scaling intermediate results.
if( xmax>bignum ) then
! scale x so that its components are less than or equal to
! bignum in absolute value.
scale = bignum / xmax
call stdlib${ii}$_dscal( n, scale, x, 1_${ik}$ )
xmax = bignum
end if
if( notran ) then
! solve a * x = b
loop_110: do j = jfirst, jlast, jinc
! compute x(j) = b(j) / a(j,j), scaling x if necessary.
xj = abs( x( j ) )
if( nounit ) then
tjjs = a( j, j )*tscal
else
tjjs = tscal
if( tscal==one )go to 100
end if
tjj = abs( tjjs )
if( tjj>smlnum ) then
! abs(a(j,j)) > smlnum:
if( tjj<one ) then
if( xj>tjj*bignum ) then
! scale x by 1/b(j).
rec = one / xj
call stdlib${ii}$_dscal( n, rec, x, 1_${ik}$ )
scale = scale*rec
xmax = xmax*rec
end if
end if
x( j ) = x( j ) / tjjs
xj = abs( x( j ) )
else if( tjj>zero ) then
! 0 < abs(a(j,j)) <= smlnum:
if( xj>tjj*bignum ) then
! scale x by (1/abs(x(j)))*abs(a(j,j))*bignum
! to avoid overflow when dividing by a(j,j).
rec = ( tjj*bignum ) / xj
if( cnorm( j )>one ) then
! scale by 1/cnorm(j) to avoid overflow when
! multiplying x(j) times column j.
rec = rec / cnorm( j )
end if
call stdlib${ii}$_dscal( n, rec, x, 1_${ik}$ )
scale = scale*rec
xmax = xmax*rec
end if
x( j ) = x( j ) / tjjs
xj = abs( x( j ) )
else
! a(j,j) = 0: set x(1:n) = 0, x(j) = 1, and
! scale = 0, and compute a solution to a*x = 0.
do i = 1, n
x( i ) = zero
end do
x( j ) = one
xj = one
scale = zero
xmax = zero
end if
100 continue
! scale x if necessary to avoid overflow when adding a
! multiple of column j of a.
if( xj>one ) then
rec = one / xj
if( cnorm( j )>( bignum-xmax )*rec ) then
! scale x by 1/(2*abs(x(j))).
rec = rec*half
call stdlib${ii}$_dscal( n, rec, x, 1_${ik}$ )
scale = scale*rec
end if
else if( xj*cnorm( j )>( bignum-xmax ) ) then
! scale x by 1/2.
call stdlib${ii}$_dscal( n, half, x, 1_${ik}$ )
scale = scale*half
end if
if( upper ) then
if( j>1_${ik}$ ) then
! compute the update
! x(1:j-1) := x(1:j-1) - x(j) * a(1:j-1,j)
call stdlib${ii}$_daxpy( j-1, -x( j )*tscal, a( 1_${ik}$, j ), 1_${ik}$, x,1_${ik}$ )
i = stdlib${ii}$_idamax( j-1, x, 1_${ik}$ )
xmax = abs( x( i ) )
end if
else
if( j<n ) then
! compute the update
! x(j+1:n) := x(j+1:n) - x(j) * a(j+1:n,j)
call stdlib${ii}$_daxpy( n-j, -x( j )*tscal, a( j+1, j ), 1_${ik}$,x( j+1 ), 1_${ik}$ )
i = j + stdlib${ii}$_idamax( n-j, x( j+1 ), 1_${ik}$ )
xmax = abs( x( i ) )
end if
end if
end do loop_110
else
! solve a**t * x = b
loop_160: do j = jfirst, jlast, jinc
! compute x(j) = b(j) - sum a(k,j)*x(k).
! k<>j
xj = abs( x( j ) )
uscal = tscal
rec = one / max( xmax, one )
if( cnorm( j )>( bignum-xj )*rec ) then
! if x(j) could overflow, scale x by 1/(2*xmax).
rec = rec*half
if( nounit ) then
tjjs = a( j, j )*tscal
else
tjjs = tscal
end if
tjj = abs( tjjs )
if( tjj>one ) then
! divide by a(j,j) when scaling x if a(j,j) > 1.
rec = min( one, rec*tjj )
uscal = uscal / tjjs
end if
if( rec<one ) then
call stdlib${ii}$_dscal( n, rec, x, 1_${ik}$ )
scale = scale*rec
xmax = xmax*rec
end if
end if
sumj = zero
if( uscal==one ) then
! if the scaling needed for a in the dot product is 1,
! call stdlib${ii}$_ddot to perform the dot product.
if( upper ) then
sumj = stdlib${ii}$_ddot( j-1, a( 1_${ik}$, j ), 1_${ik}$, x, 1_${ik}$ )
else if( j<n ) then
sumj = stdlib${ii}$_ddot( n-j, a( j+1, j ), 1_${ik}$, x( j+1 ), 1_${ik}$ )
end if
else
! otherwise, use in-line code for the dot product.
if( upper ) then
do i = 1, j - 1
sumj = sumj + ( a( i, j )*uscal )*x( i )
end do
else if( j<n ) then
do i = j + 1, n
sumj = sumj + ( a( i, j )*uscal )*x( i )
end do
end if
end if
if( uscal==tscal ) then
! compute x(j) := ( x(j) - sumj ) / a(j,j) if 1/a(j,j)
! was not used to scale the dotproduct.
x( j ) = x( j ) - sumj
xj = abs( x( j ) )
if( nounit ) then
tjjs = a( j, j )*tscal
else
tjjs = tscal
if( tscal==one )go to 150
end if
! compute x(j) = x(j) / a(j,j), scaling if necessary.
tjj = abs( tjjs )
if( tjj>smlnum ) then
! abs(a(j,j)) > smlnum:
if( tjj<one ) then
if( xj>tjj*bignum ) then
! scale x by 1/abs(x(j)).
rec = one / xj
call stdlib${ii}$_dscal( n, rec, x, 1_${ik}$ )
scale = scale*rec
xmax = xmax*rec
end if
end if
x( j ) = x( j ) / tjjs
else if( tjj>zero ) then
! 0 < abs(a(j,j)) <= smlnum:
if( xj>tjj*bignum ) then
! scale x by (1/abs(x(j)))*abs(a(j,j))*bignum.
rec = ( tjj*bignum ) / xj
call stdlib${ii}$_dscal( n, rec, x, 1_${ik}$ )
scale = scale*rec
xmax = xmax*rec
end if
x( j ) = x( j ) / tjjs
else
! a(j,j) = 0: set x(1:n) = 0, x(j) = 1, and
! scale = 0, and compute a solution to a**t*x = 0.
do i = 1, n
x( i ) = zero
end do
x( j ) = one
scale = zero
xmax = zero
end if
150 continue
else
! compute x(j) := x(j) / a(j,j) - sumj if the dot
! product has already been divided by 1/a(j,j).
x( j ) = x( j ) / tjjs - sumj
end if
xmax = max( xmax, abs( x( j ) ) )
end do loop_160
end if
scale = scale / tscal
end if
! scale the column norms by 1/tscal for return.
if( tscal/=one ) then
call stdlib${ii}$_dscal( n, one / tscal, cnorm, 1_${ik}$ )
end if
return
end subroutine stdlib${ii}$_dlatrs
#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
pure module subroutine stdlib${ii}$_${ri}$latrs( uplo, trans, diag, normin, n, a, lda, x, scale,cnorm, info )
!! DLATRS: solves one of the triangular systems
!! A *x = s*b or A**T *x = s*b
!! with scaling to prevent overflow. Here A is an upper or lower
!! triangular matrix, A**T denotes the transpose of A, x and b are
!! n-element vectors, and s is a scaling factor, usually less than
!! or equal to 1, chosen so that the components of x will be less than
!! the overflow threshold. If the unscaled problem will not cause
!! overflow, the Level 2 BLAS routine DTRSV is called. If the matrix A
!! is singular (A(j,j) = 0 for some j), then s is set to 0 and a
!! non-trivial solution to A*x = 0 is returned.
! -- lapack auxiliary routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_${rk}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: diag, normin, trans, uplo
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: lda, n
real(${rk}$), intent(out) :: scale
! Array Arguments
real(${rk}$), intent(in) :: a(lda,*)
real(${rk}$), intent(inout) :: cnorm(*), x(*)
! =====================================================================
! Local Scalars
logical(lk) :: notran, nounit, upper
integer(${ik}$) :: i, imax, j, jfirst, jinc, jlast
real(${rk}$) :: bignum, grow, rec, smlnum, sumj, tjj, tjjs, tmax, tscal, uscal, xbnd, xj, &
xmax
! Intrinsic Functions
! Executable Statements
info = 0_${ik}$
upper = stdlib_lsame( uplo, 'U' )
notran = stdlib_lsame( trans, 'N' )
nounit = stdlib_lsame( diag, 'N' )
! test the input parameters.
if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then
info = -1_${ik}$
else if( .not.notran .and. .not.stdlib_lsame( trans, 'T' ) .and. .not.stdlib_lsame( &
trans, 'C' ) ) then
info = -2_${ik}$
else if( .not.nounit .and. .not.stdlib_lsame( diag, 'U' ) ) then
info = -3_${ik}$
else if( .not.stdlib_lsame( normin, 'Y' ) .and. .not.stdlib_lsame( normin, 'N' ) ) &
then
info = -4_${ik}$
else if( n<0_${ik}$ ) then
info = -5_${ik}$
else if( lda<max( 1_${ik}$, n ) ) then
info = -7_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'DLATRS', -info )
return
end if
! quick return if possible
if( n==0 )return
! determine machine dependent parameters to control overflow.
smlnum = stdlib${ii}$_${ri}$lamch( 'SAFE MINIMUM' ) / stdlib${ii}$_${ri}$lamch( 'PRECISION' )
bignum = one / smlnum
scale = one
if( stdlib_lsame( normin, 'N' ) ) then
! compute the 1-norm of each column, not including the diagonal.
if( upper ) then
! a is upper triangular.
do j = 1, n
cnorm( j ) = stdlib${ii}$_${ri}$asum( j-1, a( 1_${ik}$, j ), 1_${ik}$ )
end do
else
! a is lower triangular.
do j = 1, n - 1
cnorm( j ) = stdlib${ii}$_${ri}$asum( n-j, a( j+1, j ), 1_${ik}$ )
end do
cnorm( n ) = zero
end if
end if
! scale the column norms by tscal if the maximum element in cnorm is
! greater than bignum.
imax = stdlib${ii}$_i${ri}$amax( n, cnorm, 1_${ik}$ )
tmax = cnorm( imax )
if( tmax<=bignum ) then
tscal = one
else
tscal = one / ( smlnum*tmax )
call stdlib${ii}$_${ri}$scal( n, tscal, cnorm, 1_${ik}$ )
end if
! compute a bound on the computed solution vector to see if the
! level 2 blas routine stdlib${ii}$_${ri}$trsv can be used.
j = stdlib${ii}$_i${ri}$amax( n, x, 1_${ik}$ )
xmax = abs( x( j ) )
xbnd = xmax
if( notran ) then
! compute the growth in a * x = b.
if( upper ) then
jfirst = n
jlast = 1_${ik}$
jinc = -1_${ik}$
else
jfirst = 1_${ik}$
jlast = n
jinc = 1_${ik}$
end if
if( tscal/=one ) then
grow = zero
go to 50
end if
if( nounit ) then
! a is non-unit triangular.
! compute grow = 1/g(j) and xbnd = 1/m(j).
! initially, g(0) = max{x(i), i=1,...,n}.
grow = one / max( xbnd, smlnum )
xbnd = grow
do j = jfirst, jlast, jinc
! exit the loop if the growth factor is too small.
if( grow<=smlnum )go to 50
! m(j) = g(j-1) / abs(a(j,j))
tjj = abs( a( j, j ) )
xbnd = min( xbnd, min( one, tjj )*grow )
if( tjj+cnorm( j )>=smlnum ) then
! g(j) = g(j-1)*( 1 + cnorm(j) / abs(a(j,j)) )
grow = grow*( tjj / ( tjj+cnorm( j ) ) )
else
! g(j) could overflow, set grow to 0.
grow = zero
end if
end do
grow = xbnd
else
! a is unit triangular.
! compute grow = 1/g(j), where g(0) = max{x(i), i=1,...,n}.
grow = min( one, one / max( xbnd, smlnum ) )
do j = jfirst, jlast, jinc
! exit the loop if the growth factor is too small.
if( grow<=smlnum )go to 50
! g(j) = g(j-1)*( 1 + cnorm(j) )
grow = grow*( one / ( one+cnorm( j ) ) )
end do
end if
50 continue
else
! compute the growth in a**t * x = b.
if( upper ) then
jfirst = 1_${ik}$
jlast = n
jinc = 1_${ik}$
else
jfirst = n
jlast = 1_${ik}$
jinc = -1_${ik}$
end if
if( tscal/=one ) then
grow = zero
go to 80
end if
if( nounit ) then
! a is non-unit triangular.
! compute grow = 1/g(j) and xbnd = 1/m(j).
! initially, m(0) = max{x(i), i=1,...,n}.
grow = one / max( xbnd, smlnum )
xbnd = grow
do j = jfirst, jlast, jinc
! exit the loop if the growth factor is too small.
if( grow<=smlnum )go to 80
! g(j) = max( g(j-1), m(j-1)*( 1 + cnorm(j) ) )
xj = one + cnorm( j )
grow = min( grow, xbnd / xj )
! m(j) = m(j-1)*( 1 + cnorm(j) ) / abs(a(j,j))
tjj = abs( a( j, j ) )
if( xj>tjj )xbnd = xbnd*( tjj / xj )
end do
grow = min( grow, xbnd )
else
! a is unit triangular.
! compute grow = 1/g(j), where g(0) = max{x(i), i=1,...,n}.
grow = min( one, one / max( xbnd, smlnum ) )
do j = jfirst, jlast, jinc
! exit the loop if the growth factor is too small.
if( grow<=smlnum )go to 80
! g(j) = ( 1 + cnorm(j) )*g(j-1)
xj = one + cnorm( j )
grow = grow / xj
end do
end if
80 continue
end if
if( ( grow*tscal )>smlnum ) then
! use the level 2 blas solve if the reciprocal of the bound on
! elements of x is not too small.
call stdlib${ii}$_${ri}$trsv( uplo, trans, diag, n, a, lda, x, 1_${ik}$ )
else
! use a level 1 blas solve, scaling intermediate results.
if( xmax>bignum ) then
! scale x so that its components are less than or equal to
! bignum in absolute value.
scale = bignum / xmax
call stdlib${ii}$_${ri}$scal( n, scale, x, 1_${ik}$ )
xmax = bignum
end if
if( notran ) then
! solve a * x = b
loop_110: do j = jfirst, jlast, jinc
! compute x(j) = b(j) / a(j,j), scaling x if necessary.
xj = abs( x( j ) )
if( nounit ) then
tjjs = a( j, j )*tscal
else
tjjs = tscal
if( tscal==one )go to 100
end if
tjj = abs( tjjs )
if( tjj>smlnum ) then
! abs(a(j,j)) > smlnum:
if( tjj<one ) then
if( xj>tjj*bignum ) then
! scale x by 1/b(j).
rec = one / xj
call stdlib${ii}$_${ri}$scal( n, rec, x, 1_${ik}$ )
scale = scale*rec
xmax = xmax*rec
end if
end if
x( j ) = x( j ) / tjjs
xj = abs( x( j ) )
else if( tjj>zero ) then
! 0 < abs(a(j,j)) <= smlnum:
if( xj>tjj*bignum ) then
! scale x by (1/abs(x(j)))*abs(a(j,j))*bignum
! to avoid overflow when dividing by a(j,j).
rec = ( tjj*bignum ) / xj
if( cnorm( j )>one ) then
! scale by 1/cnorm(j) to avoid overflow when
! multiplying x(j) times column j.
rec = rec / cnorm( j )
end if
call stdlib${ii}$_${ri}$scal( n, rec, x, 1_${ik}$ )
scale = scale*rec
xmax = xmax*rec
end if
x( j ) = x( j ) / tjjs
xj = abs( x( j ) )
else
! a(j,j) = 0: set x(1:n) = 0, x(j) = 1, and
! scale = 0, and compute a solution to a*x = 0.
do i = 1, n
x( i ) = zero
end do
x( j ) = one
xj = one
scale = zero
xmax = zero
end if
100 continue
! scale x if necessary to avoid overflow when adding a
! multiple of column j of a.
if( xj>one ) then
rec = one / xj
if( cnorm( j )>( bignum-xmax )*rec ) then
! scale x by 1/(2*abs(x(j))).
rec = rec*half
call stdlib${ii}$_${ri}$scal( n, rec, x, 1_${ik}$ )
scale = scale*rec
end if
else if( xj*cnorm( j )>( bignum-xmax ) ) then
! scale x by 1/2.
call stdlib${ii}$_${ri}$scal( n, half, x, 1_${ik}$ )
scale = scale*half
end if
if( upper ) then
if( j>1_${ik}$ ) then
! compute the update
! x(1:j-1) := x(1:j-1) - x(j) * a(1:j-1,j)
call stdlib${ii}$_${ri}$axpy( j-1, -x( j )*tscal, a( 1_${ik}$, j ), 1_${ik}$, x,1_${ik}$ )
i = stdlib${ii}$_i${ri}$amax( j-1, x, 1_${ik}$ )
xmax = abs( x( i ) )
end if
else
if( j<n ) then
! compute the update
! x(j+1:n) := x(j+1:n) - x(j) * a(j+1:n,j)
call stdlib${ii}$_${ri}$axpy( n-j, -x( j )*tscal, a( j+1, j ), 1_${ik}$,x( j+1 ), 1_${ik}$ )
i = j + stdlib${ii}$_i${ri}$amax( n-j, x( j+1 ), 1_${ik}$ )
xmax = abs( x( i ) )
end if
end if
end do loop_110
else
! solve a**t * x = b
loop_160: do j = jfirst, jlast, jinc
! compute x(j) = b(j) - sum a(k,j)*x(k).
! k<>j
xj = abs( x( j ) )
uscal = tscal
rec = one / max( xmax, one )
if( cnorm( j )>( bignum-xj )*rec ) then
! if x(j) could overflow, scale x by 1/(2*xmax).
rec = rec*half
if( nounit ) then
tjjs = a( j, j )*tscal
else
tjjs = tscal
end if
tjj = abs( tjjs )
if( tjj>one ) then
! divide by a(j,j) when scaling x if a(j,j) > 1.
rec = min( one, rec*tjj )
uscal = uscal / tjjs
end if
if( rec<one ) then
call stdlib${ii}$_${ri}$scal( n, rec, x, 1_${ik}$ )
scale = scale*rec
xmax = xmax*rec
end if
end if
sumj = zero
if( uscal==one ) then
! if the scaling needed for a in the dot product is 1,
! call stdlib${ii}$_${ri}$dot to perform the dot product.
if( upper ) then
sumj = stdlib${ii}$_${ri}$dot( j-1, a( 1_${ik}$, j ), 1_${ik}$, x, 1_${ik}$ )
else if( j<n ) then
sumj = stdlib${ii}$_${ri}$dot( n-j, a( j+1, j ), 1_${ik}$, x( j+1 ), 1_${ik}$ )
end if
else
! otherwise, use in-line code for the dot product.
if( upper ) then
do i = 1, j - 1
sumj = sumj + ( a( i, j )*uscal )*x( i )
end do
else if( j<n ) then
do i = j + 1, n
sumj = sumj + ( a( i, j )*uscal )*x( i )
end do
end if
end if
if( uscal==tscal ) then
! compute x(j) := ( x(j) - sumj ) / a(j,j) if 1/a(j,j)
! was not used to scale the dotproduct.
x( j ) = x( j ) - sumj
xj = abs( x( j ) )
if( nounit ) then
tjjs = a( j, j )*tscal
else
tjjs = tscal
if( tscal==one )go to 150
end if
! compute x(j) = x(j) / a(j,j), scaling if necessary.
tjj = abs( tjjs )
if( tjj>smlnum ) then
! abs(a(j,j)) > smlnum:
if( tjj<one ) then
if( xj>tjj*bignum ) then
! scale x by 1/abs(x(j)).
rec = one / xj
call stdlib${ii}$_${ri}$scal( n, rec, x, 1_${ik}$ )
scale = scale*rec
xmax = xmax*rec
end if
end if
x( j ) = x( j ) / tjjs
else if( tjj>zero ) then
! 0 < abs(a(j,j)) <= smlnum:
if( xj>tjj*bignum ) then
! scale x by (1/abs(x(j)))*abs(a(j,j))*bignum.
rec = ( tjj*bignum ) / xj
call stdlib${ii}$_${ri}$scal( n, rec, x, 1_${ik}$ )
scale = scale*rec
xmax = xmax*rec
end if
x( j ) = x( j ) / tjjs
else
! a(j,j) = 0: set x(1:n) = 0, x(j) = 1, and
! scale = 0, and compute a solution to a**t*x = 0.
do i = 1, n
x( i ) = zero
end do
x( j ) = one
scale = zero
xmax = zero
end if
150 continue
else
! compute x(j) := x(j) / a(j,j) - sumj if the dot
! product has already been divided by 1/a(j,j).
x( j ) = x( j ) / tjjs - sumj
end if
xmax = max( xmax, abs( x( j ) ) )
end do loop_160
end if
scale = scale / tscal
end if
! scale the column norms by 1/tscal for return.
if( tscal/=one ) then
call stdlib${ii}$_${ri}$scal( n, one / tscal, cnorm, 1_${ik}$ )
end if
return
end subroutine stdlib${ii}$_${ri}$latrs
#:endif
#:endfor
pure module subroutine stdlib${ii}$_clatrs( uplo, trans, diag, normin, n, a, lda, x, scale,cnorm, info )
!! CLATRS solves one of the triangular systems
!! A * x = s*b, A**T * x = s*b, or A**H * x = s*b,
!! with scaling to prevent overflow. Here A is an upper or lower
!! triangular matrix, A**T denotes the transpose of A, A**H denotes the
!! conjugate transpose of A, x and b are n-element vectors, and s is a
!! scaling factor, usually less than or equal to 1, chosen so that the
!! components of x will be less than the overflow threshold. If the
!! unscaled problem will not cause overflow, the Level 2 BLAS routine
!! CTRSV is called. If the matrix A is singular (A(j,j) = 0 for some j),
!! then s is set to 0 and a non-trivial solution to A*x = 0 is returned.
! -- lapack auxiliary routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: diag, normin, trans, uplo
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: lda, n
real(sp), intent(out) :: scale
! Array Arguments
real(sp), intent(inout) :: cnorm(*)
complex(sp), intent(in) :: a(lda,*)
complex(sp), intent(inout) :: x(*)
! =====================================================================
! Local Scalars
logical(lk) :: notran, nounit, upper
integer(${ik}$) :: i, imax, j, jfirst, jinc, jlast
real(sp) :: bignum, grow, rec, smlnum, tjj, tmax, tscal, xbnd, xj, xmax
complex(sp) :: csumj, tjjs, uscal, zdum
! Intrinsic Functions
! Statement Functions
real(sp) :: cabs1, cabs2
! Statement Function Definitions
cabs1( zdum ) = abs( real( zdum,KIND=sp) ) + abs( aimag( zdum ) )
cabs2( zdum ) = abs( real( zdum,KIND=sp) / 2. ) +abs( aimag( zdum ) / 2. )
! Executable Statements
info = 0_${ik}$
upper = stdlib_lsame( uplo, 'U' )
notran = stdlib_lsame( trans, 'N' )
nounit = stdlib_lsame( diag, 'N' )
! test the input parameters.
if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then
info = -1_${ik}$
else if( .not.notran .and. .not.stdlib_lsame( trans, 'T' ) .and. .not.stdlib_lsame( &
trans, 'C' ) ) then
info = -2_${ik}$
else if( .not.nounit .and. .not.stdlib_lsame( diag, 'U' ) ) then
info = -3_${ik}$
else if( .not.stdlib_lsame( normin, 'Y' ) .and. .not.stdlib_lsame( normin, 'N' ) ) &
then
info = -4_${ik}$
else if( n<0_${ik}$ ) then
info = -5_${ik}$
else if( lda<max( 1_${ik}$, n ) ) then
info = -7_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'CLATRS', -info )
return
end if
! quick return if possible
if( n==0 )return
! determine machine dependent parameters to control overflow.
smlnum = stdlib${ii}$_slamch( 'SAFE MINIMUM' )
bignum = one / smlnum
call stdlib${ii}$_slabad( smlnum, bignum )
smlnum = smlnum / stdlib${ii}$_slamch( 'PRECISION' )
bignum = one / smlnum
scale = one
if( stdlib_lsame( normin, 'N' ) ) then
! compute the 1-norm of each column, not including the diagonal.
if( upper ) then
! a is upper triangular.
do j = 1, n
cnorm( j ) = stdlib${ii}$_scasum( j-1, a( 1_${ik}$, j ), 1_${ik}$ )
end do
else
! a is lower triangular.
do j = 1, n - 1
cnorm( j ) = stdlib${ii}$_scasum( n-j, a( j+1, j ), 1_${ik}$ )
end do
cnorm( n ) = zero
end if
end if
! scale the column norms by tscal if the maximum element in cnorm is
! greater than bignum/2.
imax = stdlib${ii}$_isamax( n, cnorm, 1_${ik}$ )
tmax = cnorm( imax )
if( tmax<=bignum*half ) then
tscal = one
else
tscal = half / ( smlnum*tmax )
call stdlib${ii}$_sscal( n, tscal, cnorm, 1_${ik}$ )
end if
! compute a bound on the computed solution vector to see if the
! level 2 blas routine stdlib${ii}$_ctrsv can be used.
xmax = zero
do j = 1, n
xmax = max( xmax, cabs2( x( j ) ) )
end do
xbnd = xmax
if( notran ) then
! compute the growth in a * x = b.
if( upper ) then
jfirst = n
jlast = 1_${ik}$
jinc = -1_${ik}$
else
jfirst = 1_${ik}$
jlast = n
jinc = 1_${ik}$
end if
if( tscal/=one ) then
grow = zero
go to 60
end if
if( nounit ) then
! a is non-unit triangular.
! compute grow = 1/g(j) and xbnd = 1/m(j).
! initially, g(0) = max{x(i), i=1,...,n}.
grow = half / max( xbnd, smlnum )
xbnd = grow
do j = jfirst, jlast, jinc
! exit the loop if the growth factor is too small.
if( grow<=smlnum )go to 60
tjjs = a( j, j )
tjj = cabs1( tjjs )
if( tjj>=smlnum ) then
! m(j) = g(j-1) / abs(a(j,j))
xbnd = min( xbnd, min( one, tjj )*grow )
else
! m(j) could overflow, set xbnd to 0.
xbnd = zero
end if
if( tjj+cnorm( j )>=smlnum ) then
! g(j) = g(j-1)*( 1 + cnorm(j) / abs(a(j,j)) )
grow = grow*( tjj / ( tjj+cnorm( j ) ) )
else
! g(j) could overflow, set grow to 0.
grow = zero
end if
end do
grow = xbnd
else
! a is unit triangular.
! compute grow = 1/g(j), where g(0) = max{x(i), i=1,...,n}.
grow = min( one, half / max( xbnd, smlnum ) )
do j = jfirst, jlast, jinc
! exit the loop if the growth factor is too small.
if( grow<=smlnum )go to 60
! g(j) = g(j-1)*( 1 + cnorm(j) )
grow = grow*( one / ( one+cnorm( j ) ) )
end do
end if
60 continue
else
! compute the growth in a**t * x = b or a**h * x = b.
if( upper ) then
jfirst = 1_${ik}$
jlast = n
jinc = 1_${ik}$
else
jfirst = n
jlast = 1_${ik}$
jinc = -1_${ik}$
end if
if( tscal/=one ) then
grow = zero
go to 90
end if
if( nounit ) then
! a is non-unit triangular.
! compute grow = 1/g(j) and xbnd = 1/m(j).
! initially, m(0) = max{x(i), i=1,...,n}.
grow = half / max( xbnd, smlnum )
xbnd = grow
do j = jfirst, jlast, jinc
! exit the loop if the growth factor is too small.
if( grow<=smlnum )go to 90
! g(j) = max( g(j-1), m(j-1)*( 1 + cnorm(j) ) )
xj = one + cnorm( j )
grow = min( grow, xbnd / xj )
tjjs = a( j, j )
tjj = cabs1( tjjs )
if( tjj>=smlnum ) then
! m(j) = m(j-1)*( 1 + cnorm(j) ) / abs(a(j,j))
if( xj>tjj )xbnd = xbnd*( tjj / xj )
else
! m(j) could overflow, set xbnd to 0.
xbnd = zero
end if
end do
grow = min( grow, xbnd )
else
! a is unit triangular.
! compute grow = 1/g(j), where g(0) = max{x(i), i=1,...,n}.
grow = min( one, half / max( xbnd, smlnum ) )
do j = jfirst, jlast, jinc
! exit the loop if the growth factor is too small.
if( grow<=smlnum )go to 90
! g(j) = ( 1 + cnorm(j) )*g(j-1)
xj = one + cnorm( j )
grow = grow / xj
end do
end if
90 continue
end if
if( ( grow*tscal )>smlnum ) then
! use the level 2 blas solve if the reciprocal of the bound on
! elements of x is not too small.
call stdlib${ii}$_ctrsv( uplo, trans, diag, n, a, lda, x, 1_${ik}$ )
else
! use a level 1 blas solve, scaling intermediate results.
if( xmax>bignum*half ) then
! scale x so that its components are less than or equal to
! bignum in absolute value.
scale = ( bignum*half ) / xmax
call stdlib${ii}$_csscal( n, scale, x, 1_${ik}$ )
xmax = bignum
else
xmax = xmax*two
end if
if( notran ) then
! solve a * x = b
loop_110: do j = jfirst, jlast, jinc
! compute x(j) = b(j) / a(j,j), scaling x if necessary.
xj = cabs1( x( j ) )
if( nounit ) then
tjjs = a( j, j )*tscal
else
tjjs = tscal
if( tscal==one )go to 105
end if
tjj = cabs1( tjjs )
if( tjj>smlnum ) then
! abs(a(j,j)) > smlnum:
if( tjj<one ) then
if( xj>tjj*bignum ) then
! scale x by 1/b(j).
rec = one / xj
call stdlib${ii}$_csscal( n, rec, x, 1_${ik}$ )
scale = scale*rec
xmax = xmax*rec
end if
end if
x( j ) = stdlib${ii}$_cladiv( x( j ), tjjs )
xj = cabs1( x( j ) )
else if( tjj>zero ) then
! 0 < abs(a(j,j)) <= smlnum:
if( xj>tjj*bignum ) then
! scale x by (1/abs(x(j)))*abs(a(j,j))*bignum
! to avoid overflow when dividing by a(j,j).
rec = ( tjj*bignum ) / xj
if( cnorm( j )>one ) then
! scale by 1/cnorm(j) to avoid overflow when
! multiplying x(j) times column j.
rec = rec / cnorm( j )
end if
call stdlib${ii}$_csscal( n, rec, x, 1_${ik}$ )
scale = scale*rec
xmax = xmax*rec
end if
x( j ) = stdlib${ii}$_cladiv( x( j ), tjjs )
xj = cabs1( x( j ) )
else
! a(j,j) = 0: set x(1:n) = 0, x(j) = 1, and
! scale = 0, and compute a solution to a*x = 0.
do i = 1, n
x( i ) = zero
end do
x( j ) = one
xj = one
scale = zero
xmax = zero
end if
105 continue
! scale x if necessary to avoid overflow when adding a
! multiple of column j of a.
if( xj>one ) then
rec = one / xj
if( cnorm( j )>( bignum-xmax )*rec ) then
! scale x by 1/(2*abs(x(j))).
rec = rec*half
call stdlib${ii}$_csscal( n, rec, x, 1_${ik}$ )
scale = scale*rec
end if
else if( xj*cnorm( j )>( bignum-xmax ) ) then
! scale x by 1/2.
call stdlib${ii}$_csscal( n, half, x, 1_${ik}$ )
scale = scale*half
end if
if( upper ) then
if( j>1_${ik}$ ) then
! compute the update
! x(1:j-1) := x(1:j-1) - x(j) * a(1:j-1,j)
call stdlib${ii}$_caxpy( j-1, -x( j )*tscal, a( 1_${ik}$, j ), 1_${ik}$, x,1_${ik}$ )
i = stdlib${ii}$_icamax( j-1, x, 1_${ik}$ )
xmax = cabs1( x( i ) )
end if
else
if( j<n ) then
! compute the update
! x(j+1:n) := x(j+1:n) - x(j) * a(j+1:n,j)
call stdlib${ii}$_caxpy( n-j, -x( j )*tscal, a( j+1, j ), 1_${ik}$,x( j+1 ), 1_${ik}$ )
i = j + stdlib${ii}$_icamax( n-j, x( j+1 ), 1_${ik}$ )
xmax = cabs1( x( i ) )
end if
end if
end do loop_110
else if( stdlib_lsame( trans, 'T' ) ) then
! solve a**t * x = b
loop_150: do j = jfirst, jlast, jinc
! compute x(j) = b(j) - sum a(k,j)*x(k).
! k<>j
xj = cabs1( x( j ) )
uscal = tscal
rec = one / max( xmax, one )
if( cnorm( j )>( bignum-xj )*rec ) then
! if x(j) could overflow, scale x by 1/(2*xmax).
rec = rec*half
if( nounit ) then
tjjs = a( j, j )*tscal
else
tjjs = tscal
end if
tjj = cabs1( tjjs )
if( tjj>one ) then
! divide by a(j,j) when scaling x if a(j,j) > 1.
rec = min( one, rec*tjj )
uscal = stdlib${ii}$_cladiv( uscal, tjjs )
end if
if( rec<one ) then
call stdlib${ii}$_csscal( n, rec, x, 1_${ik}$ )
scale = scale*rec
xmax = xmax*rec
end if
end if
csumj = zero
if( uscal==cmplx( one,KIND=sp) ) then
! if the scaling needed for a in the dot product is 1,
! call stdlib${ii}$_cdotu to perform the dot product.
if( upper ) then
csumj = stdlib${ii}$_cdotu( j-1, a( 1_${ik}$, j ), 1_${ik}$, x, 1_${ik}$ )
else if( j<n ) then
csumj = stdlib${ii}$_cdotu( n-j, a( j+1, j ), 1_${ik}$, x( j+1 ), 1_${ik}$ )
end if
else
! otherwise, use in-line code for the dot product.
if( upper ) then
do i = 1, j - 1
csumj = csumj + ( a( i, j )*uscal )*x( i )
end do
else if( j<n ) then
do i = j + 1, n
csumj = csumj + ( a( i, j )*uscal )*x( i )
end do
end if
end if
if( uscal==cmplx( tscal,KIND=sp) ) then
! compute x(j) := ( x(j) - csumj ) / a(j,j) if 1/a(j,j)
! was not used to scale the dotproduct.
x( j ) = x( j ) - csumj
xj = cabs1( x( j ) )
if( nounit ) then
tjjs = a( j, j )*tscal
else
tjjs = tscal
if( tscal==one )go to 145
end if
! compute x(j) = x(j) / a(j,j), scaling if necessary.
tjj = cabs1( tjjs )
if( tjj>smlnum ) then
! abs(a(j,j)) > smlnum:
if( tjj<one ) then
if( xj>tjj*bignum ) then
! scale x by 1/abs(x(j)).
rec = one / xj
call stdlib${ii}$_csscal( n, rec, x, 1_${ik}$ )
scale = scale*rec
xmax = xmax*rec
end if
end if
x( j ) = stdlib${ii}$_cladiv( x( j ), tjjs )
else if( tjj>zero ) then
! 0 < abs(a(j,j)) <= smlnum:
if( xj>tjj*bignum ) then
! scale x by (1/abs(x(j)))*abs(a(j,j))*bignum.
rec = ( tjj*bignum ) / xj
call stdlib${ii}$_csscal( n, rec, x, 1_${ik}$ )
scale = scale*rec
xmax = xmax*rec
end if
x( j ) = stdlib${ii}$_cladiv( x( j ), tjjs )
else
! a(j,j) = 0: set x(1:n) = 0, x(j) = 1, and
! scale = 0 and compute a solution to a**t *x = 0.
do i = 1, n
x( i ) = zero
end do
x( j ) = one
scale = zero
xmax = zero
end if
145 continue
else
! compute x(j) := x(j) / a(j,j) - csumj if the dot
! product has already been divided by 1/a(j,j).
x( j ) = stdlib${ii}$_cladiv( x( j ), tjjs ) - csumj
end if
xmax = max( xmax, cabs1( x( j ) ) )
end do loop_150
else
! solve a**h * x = b
loop_190: do j = jfirst, jlast, jinc
! compute x(j) = b(j) - sum a(k,j)*x(k).
! k<>j
xj = cabs1( x( j ) )
uscal = tscal
rec = one / max( xmax, one )
if( cnorm( j )>( bignum-xj )*rec ) then
! if x(j) could overflow, scale x by 1/(2*xmax).
rec = rec*half
if( nounit ) then
tjjs = conjg( a( j, j ) )*tscal
else
tjjs = tscal
end if
tjj = cabs1( tjjs )
if( tjj>one ) then
! divide by a(j,j) when scaling x if a(j,j) > 1.
rec = min( one, rec*tjj )
uscal = stdlib${ii}$_cladiv( uscal, tjjs )
end if
if( rec<one ) then
call stdlib${ii}$_csscal( n, rec, x, 1_${ik}$ )
scale = scale*rec
xmax = xmax*rec
end if
end if
csumj = zero
if( uscal==cmplx( one,KIND=sp) ) then
! if the scaling needed for a in the dot product is 1,
! call stdlib${ii}$_cdotc to perform the dot product.
if( upper ) then
csumj = stdlib${ii}$_cdotc( j-1, a( 1_${ik}$, j ), 1_${ik}$, x, 1_${ik}$ )
else if( j<n ) then
csumj = stdlib${ii}$_cdotc( n-j, a( j+1, j ), 1_${ik}$, x( j+1 ), 1_${ik}$ )
end if
else
! otherwise, use in-line code for the dot product.
if( upper ) then
do i = 1, j - 1
csumj = csumj + ( conjg( a( i, j ) )*uscal )*x( i )
end do
else if( j<n ) then
do i = j + 1, n
csumj = csumj + ( conjg( a( i, j ) )*uscal )*x( i )
end do
end if
end if
if( uscal==cmplx( tscal,KIND=sp) ) then
! compute x(j) := ( x(j) - csumj ) / a(j,j) if 1/a(j,j)
! was not used to scale the dotproduct.
x( j ) = x( j ) - csumj
xj = cabs1( x( j ) )
if( nounit ) then
tjjs = conjg( a( j, j ) )*tscal
else
tjjs = tscal
if( tscal==one )go to 185
end if
! compute x(j) = x(j) / a(j,j), scaling if necessary.
tjj = cabs1( tjjs )
if( tjj>smlnum ) then
! abs(a(j,j)) > smlnum:
if( tjj<one ) then
if( xj>tjj*bignum ) then
! scale x by 1/abs(x(j)).
rec = one / xj
call stdlib${ii}$_csscal( n, rec, x, 1_${ik}$ )
scale = scale*rec
xmax = xmax*rec
end if
end if
x( j ) = stdlib${ii}$_cladiv( x( j ), tjjs )
else if( tjj>zero ) then
! 0 < abs(a(j,j)) <= smlnum:
if( xj>tjj*bignum ) then
! scale x by (1/abs(x(j)))*abs(a(j,j))*bignum.
rec = ( tjj*bignum ) / xj
call stdlib${ii}$_csscal( n, rec, x, 1_${ik}$ )
scale = scale*rec
xmax = xmax*rec
end if
x( j ) = stdlib${ii}$_cladiv( x( j ), tjjs )
else
! a(j,j) = 0: set x(1:n) = 0, x(j) = 1, and
! scale = 0 and compute a solution to a**h *x = 0.
do i = 1, n
x( i ) = zero
end do
x( j ) = one
scale = zero
xmax = zero
end if
185 continue
else
! compute x(j) := x(j) / a(j,j) - csumj if the dot
! product has already been divided by 1/a(j,j).
x( j ) = stdlib${ii}$_cladiv( x( j ), tjjs ) - csumj
end if
xmax = max( xmax, cabs1( x( j ) ) )
end do loop_190
end if
scale = scale / tscal
end if
! scale the column norms by 1/tscal for return.
if( tscal/=one ) then
call stdlib${ii}$_sscal( n, one / tscal, cnorm, 1_${ik}$ )
end if
return
end subroutine stdlib${ii}$_clatrs
pure module subroutine stdlib${ii}$_zlatrs( uplo, trans, diag, normin, n, a, lda, x, scale,cnorm, info )
!! ZLATRS solves one of the triangular systems
!! A * x = s*b, A**T * x = s*b, or A**H * x = s*b,
!! with scaling to prevent overflow. Here A is an upper or lower
!! triangular matrix, A**T denotes the transpose of A, A**H denotes the
!! conjugate transpose of A, x and b are n-element vectors, and s is a
!! scaling factor, usually less than or equal to 1, chosen so that the
!! components of x will be less than the overflow threshold. If the
!! unscaled problem will not cause overflow, the Level 2 BLAS routine
!! ZTRSV is called. If the matrix A is singular (A(j,j) = 0 for some j),
!! then s is set to 0 and a non-trivial solution to A*x = 0 is returned.
! -- lapack auxiliary routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: diag, normin, trans, uplo
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: lda, n
real(dp), intent(out) :: scale
! Array Arguments
real(dp), intent(inout) :: cnorm(*)
complex(dp), intent(in) :: a(lda,*)
complex(dp), intent(inout) :: x(*)
! =====================================================================
! Local Scalars
logical(lk) :: notran, nounit, upper
integer(${ik}$) :: i, imax, j, jfirst, jinc, jlast
real(dp) :: bignum, grow, rec, smlnum, tjj, tmax, tscal, xbnd, xj, xmax
complex(dp) :: csumj, tjjs, uscal, zdum
! Intrinsic Functions
! Statement Functions
real(dp) :: cabs1, cabs2
! Statement Function Definitions
cabs1( zdum ) = abs( real( zdum,KIND=dp) ) + abs( aimag( zdum ) )
cabs2( zdum ) = abs( real( zdum,KIND=dp) / 2._dp ) +abs( aimag( zdum ) / 2._dp )
! Executable Statements
info = 0_${ik}$
upper = stdlib_lsame( uplo, 'U' )
notran = stdlib_lsame( trans, 'N' )
nounit = stdlib_lsame( diag, 'N' )
! test the input parameters.
if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then
info = -1_${ik}$
else if( .not.notran .and. .not.stdlib_lsame( trans, 'T' ) .and. .not.stdlib_lsame( &
trans, 'C' ) ) then
info = -2_${ik}$
else if( .not.nounit .and. .not.stdlib_lsame( diag, 'U' ) ) then
info = -3_${ik}$
else if( .not.stdlib_lsame( normin, 'Y' ) .and. .not.stdlib_lsame( normin, 'N' ) ) &
then
info = -4_${ik}$
else if( n<0_${ik}$ ) then
info = -5_${ik}$
else if( lda<max( 1_${ik}$, n ) ) then
info = -7_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'ZLATRS', -info )
return
end if
! quick return if possible
if( n==0 )return
! determine machine dependent parameters to control overflow.
smlnum = stdlib${ii}$_dlamch( 'SAFE MINIMUM' )
bignum = one / smlnum
call stdlib${ii}$_dlabad( smlnum, bignum )
smlnum = smlnum / stdlib${ii}$_dlamch( 'PRECISION' )
bignum = one / smlnum
scale = one
if( stdlib_lsame( normin, 'N' ) ) then
! compute the 1-norm of each column, not including the diagonal.
if( upper ) then
! a is upper triangular.
do j = 1, n
cnorm( j ) = stdlib${ii}$_dzasum( j-1, a( 1_${ik}$, j ), 1_${ik}$ )
end do
else
! a is lower triangular.
do j = 1, n - 1
cnorm( j ) = stdlib${ii}$_dzasum( n-j, a( j+1, j ), 1_${ik}$ )
end do
cnorm( n ) = zero
end if
end if
! scale the column norms by tscal if the maximum element in cnorm is
! greater than bignum/2.
imax = stdlib${ii}$_idamax( n, cnorm, 1_${ik}$ )
tmax = cnorm( imax )
if( tmax<=bignum*half ) then
tscal = one
else
tscal = half / ( smlnum*tmax )
call stdlib${ii}$_dscal( n, tscal, cnorm, 1_${ik}$ )
end if
! compute a bound on the computed solution vector to see if the
! level 2 blas routine stdlib${ii}$_ztrsv can be used.
xmax = zero
do j = 1, n
xmax = max( xmax, cabs2( x( j ) ) )
end do
xbnd = xmax
if( notran ) then
! compute the growth in a * x = b.
if( upper ) then
jfirst = n
jlast = 1_${ik}$
jinc = -1_${ik}$
else
jfirst = 1_${ik}$
jlast = n
jinc = 1_${ik}$
end if
if( tscal/=one ) then
grow = zero
go to 60
end if
if( nounit ) then
! a is non-unit triangular.
! compute grow = 1/g(j) and xbnd = 1/m(j).
! initially, g(0) = max{x(i), i=1,...,n}.
grow = half / max( xbnd, smlnum )
xbnd = grow
do j = jfirst, jlast, jinc
! exit the loop if the growth factor is too small.
if( grow<=smlnum )go to 60
tjjs = a( j, j )
tjj = cabs1( tjjs )
if( tjj>=smlnum ) then
! m(j) = g(j-1) / abs(a(j,j))
xbnd = min( xbnd, min( one, tjj )*grow )
else
! m(j) could overflow, set xbnd to 0.
xbnd = zero
end if
if( tjj+cnorm( j )>=smlnum ) then
! g(j) = g(j-1)*( 1 + cnorm(j) / abs(a(j,j)) )
grow = grow*( tjj / ( tjj+cnorm( j ) ) )
else
! g(j) could overflow, set grow to 0.
grow = zero
end if
end do
grow = xbnd
else
! a is unit triangular.
! compute grow = 1/g(j), where g(0) = max{x(i), i=1,...,n}.
grow = min( one, half / max( xbnd, smlnum ) )
do j = jfirst, jlast, jinc
! exit the loop if the growth factor is too small.
if( grow<=smlnum )go to 60
! g(j) = g(j-1)*( 1 + cnorm(j) )
grow = grow*( one / ( one+cnorm( j ) ) )
end do
end if
60 continue
else
! compute the growth in a**t * x = b or a**h * x = b.
if( upper ) then
jfirst = 1_${ik}$
jlast = n
jinc = 1_${ik}$
else
jfirst = n
jlast = 1_${ik}$
jinc = -1_${ik}$
end if
if( tscal/=one ) then
grow = zero
go to 90
end if
if( nounit ) then
! a is non-unit triangular.
! compute grow = 1/g(j) and xbnd = 1/m(j).
! initially, m(0) = max{x(i), i=1,...,n}.
grow = half / max( xbnd, smlnum )
xbnd = grow
do j = jfirst, jlast, jinc
! exit the loop if the growth factor is too small.
if( grow<=smlnum )go to 90
! g(j) = max( g(j-1), m(j-1)*( 1 + cnorm(j) ) )
xj = one + cnorm( j )
grow = min( grow, xbnd / xj )
tjjs = a( j, j )
tjj = cabs1( tjjs )
if( tjj>=smlnum ) then
! m(j) = m(j-1)*( 1 + cnorm(j) ) / abs(a(j,j))
if( xj>tjj )xbnd = xbnd*( tjj / xj )
else
! m(j) could overflow, set xbnd to 0.
xbnd = zero
end if
end do
grow = min( grow, xbnd )
else
! a is unit triangular.
! compute grow = 1/g(j), where g(0) = max{x(i), i=1,...,n}.
grow = min( one, half / max( xbnd, smlnum ) )
do j = jfirst, jlast, jinc
! exit the loop if the growth factor is too small.
if( grow<=smlnum )go to 90
! g(j) = ( 1 + cnorm(j) )*g(j-1)
xj = one + cnorm( j )
grow = grow / xj
end do
end if
90 continue
end if
if( ( grow*tscal )>smlnum ) then
! use the level 2 blas solve if the reciprocal of the bound on
! elements of x is not too small.
call stdlib${ii}$_ztrsv( uplo, trans, diag, n, a, lda, x, 1_${ik}$ )
else
! use a level 1 blas solve, scaling intermediate results.
if( xmax>bignum*half ) then
! scale x so that its components are less than or equal to
! bignum in absolute value.
scale = ( bignum*half ) / xmax
call stdlib${ii}$_zdscal( n, scale, x, 1_${ik}$ )
xmax = bignum
else
xmax = xmax*two
end if
if( notran ) then
! solve a * x = b
loop_120: do j = jfirst, jlast, jinc
! compute x(j) = b(j) / a(j,j), scaling x if necessary.
xj = cabs1( x( j ) )
if( nounit ) then
tjjs = a( j, j )*tscal
else
tjjs = tscal
if( tscal==one )go to 110
end if
tjj = cabs1( tjjs )
if( tjj>smlnum ) then
! abs(a(j,j)) > smlnum:
if( tjj<one ) then
if( xj>tjj*bignum ) then
! scale x by 1/b(j).
rec = one / xj
call stdlib${ii}$_zdscal( n, rec, x, 1_${ik}$ )
scale = scale*rec
xmax = xmax*rec
end if
end if
x( j ) = stdlib${ii}$_zladiv( x( j ), tjjs )
xj = cabs1( x( j ) )
else if( tjj>zero ) then
! 0 < abs(a(j,j)) <= smlnum:
if( xj>tjj*bignum ) then
! scale x by (1/abs(x(j)))*abs(a(j,j))*bignum
! to avoid overflow when dividing by a(j,j).
rec = ( tjj*bignum ) / xj
if( cnorm( j )>one ) then
! scale by 1/cnorm(j) to avoid overflow when
! multiplying x(j) times column j.
rec = rec / cnorm( j )
end if
call stdlib${ii}$_zdscal( n, rec, x, 1_${ik}$ )
scale = scale*rec
xmax = xmax*rec
end if
x( j ) = stdlib${ii}$_zladiv( x( j ), tjjs )
xj = cabs1( x( j ) )
else
! a(j,j) = 0: set x(1:n) = 0, x(j) = 1, and
! scale = 0, and compute a solution to a*x = 0.
do i = 1, n
x( i ) = zero
end do
x( j ) = one
xj = one
scale = zero
xmax = zero
end if
110 continue
! scale x if necessary to avoid overflow when adding a
! multiple of column j of a.
if( xj>one ) then
rec = one / xj
if( cnorm( j )>( bignum-xmax )*rec ) then
! scale x by 1/(2*abs(x(j))).
rec = rec*half
call stdlib${ii}$_zdscal( n, rec, x, 1_${ik}$ )
scale = scale*rec
end if
else if( xj*cnorm( j )>( bignum-xmax ) ) then
! scale x by 1/2.
call stdlib${ii}$_zdscal( n, half, x, 1_${ik}$ )
scale = scale*half
end if
if( upper ) then
if( j>1_${ik}$ ) then
! compute the update
! x(1:j-1) := x(1:j-1) - x(j) * a(1:j-1,j)
call stdlib${ii}$_zaxpy( j-1, -x( j )*tscal, a( 1_${ik}$, j ), 1_${ik}$, x,1_${ik}$ )
i = stdlib${ii}$_izamax( j-1, x, 1_${ik}$ )
xmax = cabs1( x( i ) )
end if
else
if( j<n ) then
! compute the update
! x(j+1:n) := x(j+1:n) - x(j) * a(j+1:n,j)
call stdlib${ii}$_zaxpy( n-j, -x( j )*tscal, a( j+1, j ), 1_${ik}$,x( j+1 ), 1_${ik}$ )
i = j + stdlib${ii}$_izamax( n-j, x( j+1 ), 1_${ik}$ )
xmax = cabs1( x( i ) )
end if
end if
end do loop_120
else if( stdlib_lsame( trans, 'T' ) ) then
! solve a**t * x = b
loop_170: do j = jfirst, jlast, jinc
! compute x(j) = b(j) - sum a(k,j)*x(k).
! k<>j
xj = cabs1( x( j ) )
uscal = tscal
rec = one / max( xmax, one )
if( cnorm( j )>( bignum-xj )*rec ) then
! if x(j) could overflow, scale x by 1/(2*xmax).
rec = rec*half
if( nounit ) then
tjjs = a( j, j )*tscal
else
tjjs = tscal
end if
tjj = cabs1( tjjs )
if( tjj>one ) then
! divide by a(j,j) when scaling x if a(j,j) > 1.
rec = min( one, rec*tjj )
uscal = stdlib${ii}$_zladiv( uscal, tjjs )
end if
if( rec<one ) then
call stdlib${ii}$_zdscal( n, rec, x, 1_${ik}$ )
scale = scale*rec
xmax = xmax*rec
end if
end if
csumj = zero
if( uscal==cmplx( one,KIND=dp) ) then
! if the scaling needed for a in the dot product is 1,
! call stdlib${ii}$_zdotu to perform the dot product.
if( upper ) then
csumj = stdlib${ii}$_zdotu( j-1, a( 1_${ik}$, j ), 1_${ik}$, x, 1_${ik}$ )
else if( j<n ) then
csumj = stdlib${ii}$_zdotu( n-j, a( j+1, j ), 1_${ik}$, x( j+1 ), 1_${ik}$ )
end if
else
! otherwise, use in-line code for the dot product.
if( upper ) then
do i = 1, j - 1
csumj = csumj + ( a( i, j )*uscal )*x( i )
end do
else if( j<n ) then
do i = j + 1, n
csumj = csumj + ( a( i, j )*uscal )*x( i )
end do
end if
end if
if( uscal==cmplx( tscal,KIND=dp) ) then
! compute x(j) := ( x(j) - csumj ) / a(j,j) if 1/a(j,j)
! was not used to scale the dotproduct.
x( j ) = x( j ) - csumj
xj = cabs1( x( j ) )
if( nounit ) then
tjjs = a( j, j )*tscal
else
tjjs = tscal
if( tscal==one )go to 160
end if
! compute x(j) = x(j) / a(j,j), scaling if necessary.
tjj = cabs1( tjjs )
if( tjj>smlnum ) then
! abs(a(j,j)) > smlnum:
if( tjj<one ) then
if( xj>tjj*bignum ) then
! scale x by 1/abs(x(j)).
rec = one / xj
call stdlib${ii}$_zdscal( n, rec, x, 1_${ik}$ )
scale = scale*rec
xmax = xmax*rec
end if
end if
x( j ) = stdlib${ii}$_zladiv( x( j ), tjjs )
else if( tjj>zero ) then
! 0 < abs(a(j,j)) <= smlnum:
if( xj>tjj*bignum ) then
! scale x by (1/abs(x(j)))*abs(a(j,j))*bignum.
rec = ( tjj*bignum ) / xj
call stdlib${ii}$_zdscal( n, rec, x, 1_${ik}$ )
scale = scale*rec
xmax = xmax*rec
end if
x( j ) = stdlib${ii}$_zladiv( x( j ), tjjs )
else
! a(j,j) = 0: set x(1:n) = 0, x(j) = 1, and
! scale = 0 and compute a solution to a**t *x = 0.
do i = 1, n
x( i ) = zero
end do
x( j ) = one
scale = zero
xmax = zero
end if
160 continue
else
! compute x(j) := x(j) / a(j,j) - csumj if the dot
! product has already been divided by 1/a(j,j).
x( j ) = stdlib${ii}$_zladiv( x( j ), tjjs ) - csumj
end if
xmax = max( xmax, cabs1( x( j ) ) )
end do loop_170
else
! solve a**h * x = b
loop_220: do j = jfirst, jlast, jinc
! compute x(j) = b(j) - sum a(k,j)*x(k).
! k<>j
xj = cabs1( x( j ) )
uscal = tscal
rec = one / max( xmax, one )
if( cnorm( j )>( bignum-xj )*rec ) then
! if x(j) could overflow, scale x by 1/(2*xmax).
rec = rec*half
if( nounit ) then
tjjs = conjg( a( j, j ) )*tscal
else
tjjs = tscal
end if
tjj = cabs1( tjjs )
if( tjj>one ) then
! divide by a(j,j) when scaling x if a(j,j) > 1.
rec = min( one, rec*tjj )
uscal = stdlib${ii}$_zladiv( uscal, tjjs )
end if
if( rec<one ) then
call stdlib${ii}$_zdscal( n, rec, x, 1_${ik}$ )
scale = scale*rec
xmax = xmax*rec
end if
end if
csumj = zero
if( uscal==cmplx( one,KIND=dp) ) then
! if the scaling needed for a in the dot product is 1,
! call stdlib${ii}$_zdotc to perform the dot product.
if( upper ) then
csumj = stdlib${ii}$_zdotc( j-1, a( 1_${ik}$, j ), 1_${ik}$, x, 1_${ik}$ )
else if( j<n ) then
csumj = stdlib${ii}$_zdotc( n-j, a( j+1, j ), 1_${ik}$, x( j+1 ), 1_${ik}$ )
end if
else
! otherwise, use in-line code for the dot product.
if( upper ) then
do i = 1, j - 1
csumj = csumj + ( conjg( a( i, j ) )*uscal )*x( i )
end do
else if( j<n ) then
do i = j + 1, n
csumj = csumj + ( conjg( a( i, j ) )*uscal )*x( i )
end do
end if
end if
if( uscal==cmplx( tscal,KIND=dp) ) then
! compute x(j) := ( x(j) - csumj ) / a(j,j) if 1/a(j,j)
! was not used to scale the dotproduct.
x( j ) = x( j ) - csumj
xj = cabs1( x( j ) )
if( nounit ) then
tjjs = conjg( a( j, j ) )*tscal
else
tjjs = tscal
if( tscal==one )go to 210
end if
! compute x(j) = x(j) / a(j,j), scaling if necessary.
tjj = cabs1( tjjs )
if( tjj>smlnum ) then
! abs(a(j,j)) > smlnum:
if( tjj<one ) then
if( xj>tjj*bignum ) then
! scale x by 1/abs(x(j)).
rec = one / xj
call stdlib${ii}$_zdscal( n, rec, x, 1_${ik}$ )
scale = scale*rec
