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