stdlib_lapack_solve_tri_comp.fypp Source File


Source Code

#: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