#:include "common.fypp" submodule(stdlib_lapack_solve) stdlib_lapack_solve_chol_comp implicit none contains #:for ik,it,ii in LINALG_INT_KINDS_TYPES pure module subroutine stdlib${ii}$_spocon( uplo, n, a, lda, anorm, rcond, work, iwork,info ) !! SPOCON estimates the reciprocal of the condition number (in the !! 1-norm) of a real symmetric positive definite matrix using the !! Cholesky factorization A = U**T*U or A = L*L**T computed by SPOTRF. !! An estimate is obtained for norm(inv(A)), and the reciprocal of the !! condition number is computed as RCOND = 1 / (ANORM * 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) :: uplo integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, n real(sp), intent(in) :: anorm real(sp), intent(out) :: rcond ! Array Arguments integer(${ik}$), intent(out) :: iwork(*) real(sp), intent(inout) :: a(lda,*) real(sp), intent(out) :: work(*) ! ===================================================================== ! Local Scalars logical(lk) :: upper character :: normin integer(${ik}$) :: ix, kase real(sp) :: ainvnm, scale, scalel, scaleu, smlnum ! Local Arrays integer(${ik}$) :: isave(3_${ik}$) ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ upper = stdlib_lsame( uplo, 'U' ) if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -4_${ik}$ else if( anorm<zero ) then info = -5_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'SPOCON', -info ) return end if ! quick return if possible rcond = zero if( n==0_${ik}$ ) then rcond = one return else if( anorm==zero ) then return end if smlnum = stdlib${ii}$_slamch( 'SAFE MINIMUM' ) ! estimate the 1-norm of inv(a). kase = 0_${ik}$ normin = 'N' 10 continue call stdlib${ii}$_slacn2( n, work( n+1 ), work, iwork, ainvnm, kase, isave ) if( kase/=0_${ik}$ ) then if( upper ) then ! multiply by inv(u**t). call stdlib${ii}$_slatrs( 'UPPER', 'TRANSPOSE', 'NON-UNIT', normin, n, a,lda, work, & scalel, work( 2_${ik}$*n+1 ), info ) normin = 'Y' ! multiply by inv(u). call stdlib${ii}$_slatrs( 'UPPER', 'NO TRANSPOSE', 'NON-UNIT', normin, n,a, lda, work, & scaleu, work( 2_${ik}$*n+1 ), info ) else ! multiply by inv(l). call stdlib${ii}$_slatrs( 'LOWER', 'NO TRANSPOSE', 'NON-UNIT', normin, n,a, lda, work, & scalel, work( 2_${ik}$*n+1 ), info ) normin = 'Y' ! multiply by inv(l**t). call stdlib${ii}$_slatrs( 'LOWER', 'TRANSPOSE', 'NON-UNIT', normin, n, a,lda, work, & scaleu, work( 2_${ik}$*n+1 ), info ) end if ! multiply by 1/scale if doing so will not cause overflow. scale = scalel*scaleu if( scale/=one ) then ix = stdlib${ii}$_isamax( n, work, 1_${ik}$ ) if( scale<abs( work( ix ) )*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 / ainvnm ) / anorm 20 continue return end subroutine stdlib${ii}$_spocon pure module subroutine stdlib${ii}$_dpocon( uplo, n, a, lda, anorm, rcond, work, iwork,info ) !! DPOCON estimates the reciprocal of the condition number (in the !! 1-norm) of a real symmetric positive definite matrix using the !! Cholesky factorization A = U**T*U or A = L*L**T computed by DPOTRF. !! An estimate is obtained for norm(inv(A)), and the reciprocal of the !! condition number is computed as RCOND = 1 / (ANORM * 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) :: uplo integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, n real(dp), intent(in) :: anorm real(dp), intent(out) :: rcond ! Array Arguments integer(${ik}$), intent(out) :: iwork(*) real(dp), intent(inout) :: a(lda,*) real(dp), intent(out) :: work(*) ! ===================================================================== ! Local Scalars logical(lk) :: upper character :: normin integer(${ik}$) :: ix, kase real(dp) :: ainvnm, scale, scalel, scaleu, smlnum ! Local Arrays integer(${ik}$) :: isave(3_${ik}$) ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ upper = stdlib_lsame( uplo, 'U' ) if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -4_${ik}$ else if( anorm<zero ) then info = -5_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DPOCON', -info ) return end if ! quick return if possible rcond = zero if( n==0_${ik}$ ) then rcond = one return else if( anorm==zero ) then return end if smlnum = stdlib${ii}$_dlamch( 'SAFE MINIMUM' ) ! estimate the 1-norm of inv(a). kase = 0_${ik}$ normin = 'N' 10 continue call stdlib${ii}$_dlacn2( n, work( n+1 ), work, iwork, ainvnm, kase, isave ) if( kase/=0_${ik}$ ) then if( upper ) then ! multiply by inv(u**t). call stdlib${ii}$_dlatrs( 'UPPER', 'TRANSPOSE', 'NON-UNIT', normin, n, a,lda, work, & scalel, work( 2_${ik}$*n+1 ), info ) normin = 'Y' ! multiply by inv(u). call stdlib${ii}$_dlatrs( 'UPPER', 'NO TRANSPOSE', 'NON-UNIT', normin, n,a, lda, work, & scaleu, work( 2_${ik}$*n+1 ), info ) else ! multiply by inv(l). call stdlib${ii}$_dlatrs( 'LOWER', 'NO TRANSPOSE', 'NON-UNIT', normin, n,a, lda, work, & scalel, work( 2_${ik}$*n+1 ), info ) normin = 'Y' ! multiply by inv(l**t). call stdlib${ii}$_dlatrs( 'LOWER', 'TRANSPOSE', 'NON-UNIT', normin, n, a,lda, work, & scaleu, work( 2_${ik}$*n+1 ), info ) end if ! multiply by 1/scale if doing so will not cause overflow. scale = scalel*scaleu if( scale/=one ) then ix = stdlib${ii}$_idamax( n, work, 1_${ik}$ ) if( scale<abs( work( ix ) )*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 / ainvnm ) / anorm 20 continue return end subroutine stdlib${ii}$_dpocon #:for rk,rt,ri in REAL_KINDS_TYPES #:if not rk in ["sp","dp"] pure module subroutine stdlib${ii}$_${ri}$pocon( uplo, n, a, lda, anorm, rcond, work, iwork,info ) !! DPOCON: estimates the reciprocal of the condition number (in the !! 1-norm) of a real symmetric positive definite matrix using the !! Cholesky factorization A = U**T*U or A = L*L**T computed by DPOTRF. !! An estimate is obtained for norm(inv(A)), and the reciprocal of the !! condition number is computed as RCOND = 1 / (ANORM * 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) :: uplo integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, n real(${rk}$), intent(in) :: anorm real(${rk}$), intent(out) :: rcond ! Array Arguments integer(${ik}$), intent(out) :: iwork(*) real(${rk}$), intent(inout) :: a(lda,*) real(${rk}$), intent(out) :: work(*) ! ===================================================================== ! Local Scalars logical(lk) :: upper character :: normin integer(${ik}$) :: ix, kase real(${rk}$) :: ainvnm, scale, scalel, scaleu, smlnum ! Local Arrays integer(${ik}$) :: isave(3_${ik}$) ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ upper = stdlib_lsame( uplo, 'U' ) if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -4_${ik}$ else if( anorm<zero ) then info = -5_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DPOCON', -info ) return end if ! quick return if possible rcond = zero if( n==0_${ik}$ ) then rcond = one return else if( anorm==zero ) then return end if smlnum = stdlib${ii}$_${ri}$lamch( 'SAFE MINIMUM' ) ! estimate the 1-norm of inv(a). kase = 0_${ik}$ normin = 'N' 10 continue call stdlib${ii}$_${ri}$lacn2( n, work( n+1 ), work, iwork, ainvnm, kase, isave ) if( kase/=0_${ik}$ ) then if( upper ) then ! multiply by inv(u**t). call stdlib${ii}$_${ri}$latrs( 'UPPER', 'TRANSPOSE', 'NON-UNIT', normin, n, a,lda, work, & scalel, work( 2_${ik}$*n+1 ), info ) normin = 'Y' ! multiply by inv(u). call stdlib${ii}$_${ri}$latrs( 'UPPER', 'NO TRANSPOSE', 'NON-UNIT', normin, n,a, lda, work, & scaleu, work( 2_${ik}$*n+1 ), info ) else ! multiply by inv(l). call stdlib${ii}$_${ri}$latrs( 'LOWER', 'NO TRANSPOSE', 'NON-UNIT', normin, n,a, lda, work, & scalel, work( 2_${ik}$*n+1 ), info ) normin = 'Y' ! multiply by inv(l**t). call stdlib${ii}$_${ri}$latrs( 'LOWER', 'TRANSPOSE', 'NON-UNIT', normin, n, a,lda, work, & scaleu, work( 2_${ik}$*n+1 ), info ) end if ! multiply by 1/scale if doing so will not cause overflow. scale = scalel*scaleu if( scale/=one ) then ix = stdlib${ii}$_i${ri}$amax( n, work, 1_${ik}$ ) if( scale<abs( work( ix ) )*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 / ainvnm ) / anorm 20 continue return end subroutine stdlib${ii}$_${ri}$pocon #:endif #:endfor pure module subroutine stdlib${ii}$_cpocon( uplo, n, a, lda, anorm, rcond, work, rwork,info ) !! CPOCON estimates the reciprocal of the condition number (in the !! 1-norm) of a complex Hermitian positive definite matrix using the !! Cholesky factorization A = U**H*U or A = L*L**H computed by CPOTRF. !! An estimate is obtained for norm(inv(A)), and the reciprocal of the !! condition number is computed as RCOND = 1 / (ANORM * 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) :: uplo integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, n real(sp), intent(in) :: anorm real(sp), intent(out) :: rcond ! Array Arguments real(sp), intent(out) :: rwork(*) complex(sp), intent(inout) :: a(lda,*) complex(sp), intent(out) :: work(*) ! ===================================================================== ! Local Scalars logical(lk) :: upper character :: normin integer(${ik}$) :: ix, kase real(sp) :: ainvnm, scale, scalel, scaleu, smlnum 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' ) if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -4_${ik}$ else if( anorm<zero ) then info = -5_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'CPOCON', -info ) return end if ! quick return if possible rcond = zero if( n==0_${ik}$ ) then rcond = one return else if( anorm==zero ) then return end if smlnum = stdlib${ii}$_slamch( 'SAFE MINIMUM' ) ! estimate the 1-norm of inv(a). kase = 0_${ik}$ normin = 'N' 10 continue call stdlib${ii}$_clacn2( n, work( n+1 ), work, ainvnm, kase, isave ) if( kase/=0_${ik}$ ) then if( upper ) then ! multiply by inv(u**h). call stdlib${ii}$_clatrs( 'UPPER', 'CONJUGATE TRANSPOSE', 'NON-UNIT',normin, n, a, lda,& work, scalel, rwork, info ) normin = 'Y' ! multiply by inv(u). call stdlib${ii}$_clatrs( 'UPPER', 'NO TRANSPOSE', 'NON-UNIT', normin, n,a, lda, work, & scaleu, rwork, info ) else ! multiply by inv(l). call stdlib${ii}$_clatrs( 'LOWER', 'NO TRANSPOSE', 'NON-UNIT', normin, n,a, lda, work, & scalel, rwork, info ) normin = 'Y' ! multiply by inv(l**h). call stdlib${ii}$_clatrs( 'LOWER', 'CONJUGATE TRANSPOSE', 'NON-UNIT',normin, n, a, lda,& work, scaleu, rwork, info ) end if ! multiply by 1/scale if doing so will not cause overflow. scale = scalel*scaleu if( scale/=one ) then ix = stdlib${ii}$_icamax( n, work, 1_${ik}$ ) if( scale<cabs1( work( ix ) )*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 / ainvnm ) / anorm 20 continue return end subroutine stdlib${ii}$_cpocon pure module subroutine stdlib${ii}$_zpocon( uplo, n, a, lda, anorm, rcond, work, rwork,info ) !! ZPOCON estimates the reciprocal of the condition number (in the !! 1-norm) of a complex Hermitian positive definite matrix using the !! Cholesky factorization A = U**H*U or A = L*L**H computed by ZPOTRF. !! An estimate is obtained for norm(inv(A)), and the reciprocal of the !! condition number is computed as RCOND = 1 / (ANORM * 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) :: uplo integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, n real(dp), intent(in) :: anorm real(dp), intent(out) :: rcond ! Array Arguments real(dp), intent(out) :: rwork(*) complex(dp), intent(inout) :: a(lda,*) complex(dp), intent(out) :: work(*) ! ===================================================================== ! Local Scalars logical(lk) :: upper character :: normin integer(${ik}$) :: ix, kase real(dp) :: ainvnm, scale, scalel, scaleu, smlnum 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' ) if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -4_${ik}$ else if( anorm<zero ) then info = -5_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'ZPOCON', -info ) return end if ! quick return if possible rcond = zero if( n==0_${ik}$ ) then rcond = one return else if( anorm==zero ) then return end if smlnum = stdlib${ii}$_dlamch( 'SAFE MINIMUM' ) ! estimate the 1-norm of inv(a). kase = 0_${ik}$ normin = 'N' 10 continue call stdlib${ii}$_zlacn2( n, work( n+1 ), work, ainvnm, kase, isave ) if( kase/=0_${ik}$ ) then if( upper ) then ! multiply by inv(u**h). call stdlib${ii}$_zlatrs( 'UPPER', 'CONJUGATE TRANSPOSE', 'NON-UNIT',normin, n, a, lda,& work, scalel, rwork, info ) normin = 'Y' ! multiply by inv(u). call stdlib${ii}$_zlatrs( 'UPPER', 'NO TRANSPOSE', 'NON-UNIT', normin, n,a, lda, work, & scaleu, rwork, info ) else ! multiply by inv(l). call stdlib${ii}$_zlatrs( 'LOWER', 'NO TRANSPOSE', 'NON-UNIT', normin, n,a, lda, work, & scalel, rwork, info ) normin = 'Y' ! multiply by inv(l**h). call stdlib${ii}$_zlatrs( 'LOWER', 'CONJUGATE TRANSPOSE', 'NON-UNIT',normin, n, a, lda,& work, scaleu, rwork, info ) end if ! multiply by 1/scale if doing so will not cause overflow. scale = scalel*scaleu if( scale/=one ) then ix = stdlib${ii}$_izamax( n, work, 1_${ik}$ ) if( scale<cabs1( work( ix ) )*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 / ainvnm ) / anorm 20 continue return end subroutine stdlib${ii}$_zpocon #:for ck,ct,ci in CMPLX_KINDS_TYPES #:if not ck in ["sp","dp"] pure module subroutine stdlib${ii}$_${ci}$pocon( uplo, n, a, lda, anorm, rcond, work, rwork,info ) !! ZPOCON: estimates the reciprocal of the condition number (in the !! 1-norm) of a complex Hermitian positive definite matrix using the !! Cholesky factorization A = U**H*U or A = L*L**H computed by ZPOTRF. !! An estimate is obtained for norm(inv(A)), and the reciprocal of the !! condition number is computed as RCOND = 1 / (ANORM * 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) :: uplo integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, n real(${ck}$), intent(in) :: anorm real(${ck}$), intent(out) :: rcond ! Array Arguments real(${ck}$), intent(out) :: rwork(*) complex(${ck}$), intent(inout) :: a(lda,*) complex(${ck}$), intent(out) :: work(*) ! ===================================================================== ! Local Scalars logical(lk) :: upper character :: normin integer(${ik}$) :: ix, kase real(${ck}$) :: ainvnm, scale, scalel, scaleu, smlnum 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' ) if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -4_${ik}$ else if( anorm<zero ) then info = -5_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'ZPOCON', -info ) return end if ! quick return if possible rcond = zero if( n==0_${ik}$ ) then rcond = one return else if( anorm==zero ) then return end if smlnum = stdlib${ii}$_${c2ri(ci)}$lamch( 'SAFE MINIMUM' ) ! estimate the 1-norm of inv(a). kase = 0_${ik}$ normin = 'N' 10 continue call stdlib${ii}$_${ci}$lacn2( n, work( n+1 ), work, ainvnm, kase, isave ) if( kase/=0_${ik}$ ) then if( upper ) then ! multiply by inv(u**h). call stdlib${ii}$_${ci}$latrs( 'UPPER', 'CONJUGATE TRANSPOSE', 'NON-UNIT',normin, n, a, lda,& work, scalel, rwork, info ) normin = 'Y' ! multiply by inv(u). call stdlib${ii}$_${ci}$latrs( 'UPPER', 'NO TRANSPOSE', 'NON-UNIT', normin, n,a, lda, work, & scaleu, rwork, info ) else ! multiply by inv(l). call stdlib${ii}$_${ci}$latrs( 'LOWER', 'NO TRANSPOSE', 'NON-UNIT', normin, n,a, lda, work, & scalel, rwork, info ) normin = 'Y' ! multiply by inv(l**h). call stdlib${ii}$_${ci}$latrs( 'LOWER', 'CONJUGATE TRANSPOSE', 'NON-UNIT',normin, n, a, lda,& work, scaleu, rwork, info ) end if ! multiply by 1/scale if doing so will not cause overflow. scale = scalel*scaleu if( scale/=one ) then ix = stdlib${ii}$_i${ci}$amax( n, work, 1_${ik}$ ) if( scale<cabs1( work( ix ) )*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 / ainvnm ) / anorm 20 continue return end subroutine stdlib${ii}$_${ci}$pocon #:endif #:endfor pure module subroutine stdlib${ii}$_spotrf( uplo, n, a, lda, info ) !! SPOTRF computes the Cholesky factorization of a real symmetric !! positive definite matrix A. !! The factorization has the form !! A = U**T * U, if UPLO = 'U', or !! A = L * L**T, if UPLO = 'L', !! where U is an upper triangular matrix and L is lower triangular. !! This is the block version of the algorithm, calling Level 3 BLAS. ! -- 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) :: uplo integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, n ! Array Arguments real(sp), intent(inout) :: a(lda,*) ! ===================================================================== ! Local Scalars logical(lk) :: upper integer(${ik}$) :: j, jb, nb ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ upper = stdlib_lsame( uplo, 'U' ) if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -4_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'SPOTRF', -info ) return end if ! quick return if possible if( n==0 )return ! determine the block size for this environment. nb = stdlib${ii}$_ilaenv( 1_${ik}$, 'SPOTRF', uplo, n, -1_${ik}$, -1_${ik}$, -1_${ik}$ ) if( nb<=1_${ik}$ .or. nb>=n ) then ! use unblocked code. call stdlib${ii}$_spotrf2( uplo, n, a, lda, info ) else ! use blocked code. if( upper ) then ! compute the cholesky factorization a = u**t*u. do j = 1, n, nb ! update and factorize the current diagonal block and test ! for non-positive-definiteness. jb = min( nb, n-j+1 ) call stdlib${ii}$_ssyrk( 'UPPER', 'TRANSPOSE', jb, j-1, -one,a( 1_${ik}$, j ), lda, one, a(& j, j ), lda ) call stdlib${ii}$_spotrf2( 'UPPER', jb, a( j, j ), lda, info ) if( info/=0 )go to 30 if( j+jb<=n ) then ! compute the current block row. call stdlib${ii}$_sgemm( 'TRANSPOSE', 'NO TRANSPOSE', jb, n-j-jb+1,j-1, -one, a( & 1_${ik}$, j ), lda, a( 1_${ik}$, j+jb ),lda, one, a( j, j+jb ), lda ) call stdlib${ii}$_strsm( 'LEFT', 'UPPER', 'TRANSPOSE', 'NON-UNIT',jb, n-j-jb+1, & one, a( j, j ), lda,a( j, j+jb ), lda ) end if end do else ! compute the cholesky factorization a = l*l**t. do j = 1, n, nb ! update and factorize the current diagonal block and test ! for non-positive-definiteness. jb = min( nb, n-j+1 ) call stdlib${ii}$_ssyrk( 'LOWER', 'NO TRANSPOSE', jb, j-1, -one,a( j, 1_${ik}$ ), lda, one,& a( j, j ), lda ) call stdlib${ii}$_spotrf2( 'LOWER', jb, a( j, j ), lda, info ) if( info/=0 )go to 30 if( j+jb<=n ) then ! compute the current block column. call stdlib${ii}$_sgemm( 'NO TRANSPOSE', 'TRANSPOSE', n-j-jb+1, jb,j-1, -one, a( & j+jb, 1_${ik}$ ), lda, a( j, 1_${ik}$ ),lda, one, a( j+jb, j ), lda ) call stdlib${ii}$_strsm( 'RIGHT', 'LOWER', 'TRANSPOSE', 'NON-UNIT',n-j-jb+1, jb, & one, a( j, j ), lda,a( j+jb, j ), lda ) end if end do end if end if go to 40 30 continue info = info + j - 1_${ik}$ 40 continue return end subroutine stdlib${ii}$_spotrf pure module subroutine stdlib${ii}$_dpotrf( uplo, n, a, lda, info ) !! DPOTRF computes the Cholesky factorization of a real symmetric !! positive definite matrix A. !! The factorization has the form !! A = U**T * U, if UPLO = 'U', or !! A = L * L**T, if UPLO = 'L', !! where U is an upper triangular matrix and L is lower triangular. !! This is the block version of the algorithm, calling Level 3 BLAS. ! -- 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) :: uplo integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, n ! Array Arguments real(dp), intent(inout) :: a(lda,*) ! ===================================================================== ! Local Scalars logical(lk) :: upper integer(${ik}$) :: j, jb, nb ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ upper = stdlib_lsame( uplo, 'U' ) if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -4_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DPOTRF', -info ) return end if ! quick return if possible if( n==0 )return ! determine the block size for this environment. nb = stdlib${ii}$_ilaenv( 1_${ik}$, 'DPOTRF', uplo, n, -1_${ik}$, -1_${ik}$, -1_${ik}$ ) if( nb<=1_${ik}$ .or. nb>=n ) then ! use unblocked code. call stdlib${ii}$_dpotrf2( uplo, n, a, lda, info ) else ! use blocked code. if( upper ) then ! compute the cholesky factorization a = u**t*u. do j = 1, n, nb ! update and factorize the current diagonal block and test ! for non-positive-definiteness. jb = min( nb, n-j+1 ) call stdlib${ii}$_dsyrk( 'UPPER', 'TRANSPOSE', jb, j-1, -one,a( 1_${ik}$, j ), lda, one, a(& j, j ), lda ) call stdlib${ii}$_dpotrf2( 'UPPER', jb, a( j, j ), lda, info ) if( info/=0 )go to 30 if( j+jb<=n ) then ! compute the current block row. call stdlib${ii}$_dgemm( 'TRANSPOSE', 'NO TRANSPOSE', jb, n-j-jb+1,j-1, -one, a( & 1_${ik}$, j ), lda, a( 1_${ik}$, j+jb ),lda, one, a( j, j+jb ), lda ) call stdlib${ii}$_dtrsm( 'LEFT', 'UPPER', 'TRANSPOSE', 'NON-UNIT',jb, n-j-jb+1, & one, a( j, j ), lda,a( j, j+jb ), lda ) end if end do else ! compute the cholesky factorization a = l*l**t. do j = 1, n, nb ! update and factorize the current diagonal block and test ! for non-positive-definiteness. jb = min( nb, n-j+1 ) call stdlib${ii}$_dsyrk( 'LOWER', 'NO TRANSPOSE', jb, j-1, -one,a( j, 1_${ik}$ ), lda, one,& a( j, j ), lda ) call stdlib${ii}$_dpotrf2( 'LOWER', jb, a( j, j ), lda, info ) if( info/=0 )go to 30 if( j+jb<=n ) then ! compute the current block column. call stdlib${ii}$_dgemm( 'NO TRANSPOSE', 'TRANSPOSE', n-j-jb+1, jb,j-1, -one, a( & j+jb, 1_${ik}$ ), lda, a( j, 1_${ik}$ ),lda, one, a( j+jb, j ), lda ) call stdlib${ii}$_dtrsm( 'RIGHT', 'LOWER', 'TRANSPOSE', 'NON-UNIT',n-j-jb+1, jb, & one, a( j, j ), lda,a( j+jb, j ), lda ) end if end do end if end if go to 40 30 continue info = info + j - 1_${ik}$ 40 continue return end subroutine stdlib${ii}$_dpotrf #:for rk,rt,ri in REAL_KINDS_TYPES #:if not rk in ["sp","dp"] pure module subroutine stdlib${ii}$_${ri}$potrf( uplo, n, a, lda, info ) !! DPOTRF: computes the Cholesky factorization of a real symmetric !! positive definite matrix A. !! The factorization has the form !! A = U**T * U, if UPLO = 'U', or !! A = L * L**T, if UPLO = 'L', !! where U is an upper triangular matrix and L is lower triangular. !! This is the block version of the algorithm, calling Level 3 BLAS. ! -- 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) :: uplo integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, n ! Array Arguments real(${rk}$), intent(inout) :: a(lda,*) ! ===================================================================== ! Local Scalars logical(lk) :: upper integer(${ik}$) :: j, jb, nb ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ upper = stdlib_lsame( uplo, 'U' ) if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -4_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DPOTRF', -info ) return end if ! quick return if possible if( n==0 )return ! determine the block size for this environment. nb = stdlib${ii}$_ilaenv( 1_${ik}$, 'DPOTRF', uplo, n, -1_${ik}$, -1_${ik}$, -1_${ik}$ ) if( nb<=1_${ik}$ .or. nb>=n ) then ! use unblocked code. call stdlib${ii}$_${ri}$potrf2( uplo, n, a, lda, info ) else ! use blocked code. if( upper ) then ! compute the cholesky factorization a = u**t*u. do j = 1, n, nb ! update and factorize the current diagonal block and test ! for non-positive-definiteness. jb = min( nb, n-j+1 ) call stdlib${ii}$_${ri}$syrk( 'UPPER', 'TRANSPOSE', jb, j-1, -one,a( 1_${ik}$, j ), lda, one, a(& j, j ), lda ) call stdlib${ii}$_${ri}$potrf2( 'UPPER', jb, a( j, j ), lda, info ) if( info/=0 )go to 30 if( j+jb<=n ) then ! compute the current block row. call stdlib${ii}$_${ri}$gemm( 'TRANSPOSE', 'NO TRANSPOSE', jb, n-j-jb+1,j-1, -one, a( & 1_${ik}$, j ), lda, a( 1_${ik}$, j+jb ),lda, one, a( j, j+jb ), lda ) call stdlib${ii}$_${ri}$trsm( 'LEFT', 'UPPER', 'TRANSPOSE', 'NON-UNIT',jb, n-j-jb+1, & one, a( j, j ), lda,a( j, j+jb ), lda ) end if end do else ! compute the cholesky factorization a = l*l**t. do j = 1, n, nb ! update and factorize the current diagonal block and test ! for non-positive-definiteness. jb = min( nb, n-j+1 ) call stdlib${ii}$_${ri}$syrk( 'LOWER', 'NO TRANSPOSE', jb, j-1, -one,a( j, 1_${ik}$ ), lda, one,& a( j, j ), lda ) call stdlib${ii}$_${ri}$potrf2( 'LOWER', jb, a( j, j ), lda, info ) if( info/=0 )go to 30 if( j+jb<=n ) then ! compute the current block column. call stdlib${ii}$_${ri}$gemm( 'NO TRANSPOSE', 'TRANSPOSE', n-j-jb+1, jb,j-1, -one, a( & j+jb, 1_${ik}$ ), lda, a( j, 1_${ik}$ ),lda, one, a( j+jb, j ), lda ) call stdlib${ii}$_${ri}$trsm( 'RIGHT', 'LOWER', 'TRANSPOSE', 'NON-UNIT',n-j-jb+1, jb, & one, a( j, j ), lda,a( j+jb, j ), lda ) end if end do end if end if go to 40 30 continue info = info + j - 1_${ik}$ 40 continue return end subroutine stdlib${ii}$_${ri}$potrf #:endif #:endfor pure module subroutine stdlib${ii}$_cpotrf( uplo, n, a, lda, info ) !! CPOTRF computes the Cholesky factorization of a complex Hermitian !! positive definite matrix A. !! The factorization has the form !! A = U**H * U, if UPLO = 'U', or !! A = L * L**H, if UPLO = 'L', !! where U is an upper triangular matrix and L is lower triangular. !! This is the block version of the algorithm, calling Level 3 BLAS. ! -- 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) :: uplo integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, n ! Array Arguments complex(sp), intent(inout) :: a(lda,*) ! ===================================================================== ! Local Scalars logical(lk) :: upper integer(${ik}$) :: j, jb, nb ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ upper = stdlib_lsame( uplo, 'U' ) if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -4_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'CPOTRF', -info ) return end if ! quick return if possible if( n==0 )return ! determine the block size for this environment. nb = stdlib${ii}$_ilaenv( 1_${ik}$, 'CPOTRF', uplo, n, -1_${ik}$, -1_${ik}$, -1_${ik}$ ) if( nb<=1_${ik}$ .or. nb>=n ) then ! use unblocked code. call stdlib${ii}$_cpotrf2( uplo, n, a, lda, info ) else ! use blocked code. if( upper ) then ! compute the cholesky factorization a = u**h *u. do j = 1, n, nb ! update and factorize the current diagonal block and test ! for non-positive-definiteness. jb = min( nb, n-j+1 ) call stdlib${ii}$_cherk( 'UPPER', 'CONJUGATE TRANSPOSE', jb, j-1,-one, a( 1_${ik}$, j ), & lda, one, a( j, j ), lda ) call stdlib${ii}$_cpotrf2( 'UPPER', jb, a( j, j ), lda, info ) if( info/=0 )go to 30 if( j+jb<=n ) then ! compute the current block row. call stdlib${ii}$_cgemm( 'CONJUGATE TRANSPOSE', 'NO TRANSPOSE', jb,n-j-jb+1, j-1,& -cone, a( 1_${ik}$, j ), lda,a( 1_${ik}$, j+jb ), lda, cone, a( j, j+jb ),lda ) call stdlib${ii}$_ctrsm( 'LEFT', 'UPPER', 'CONJUGATE TRANSPOSE','NON-UNIT', jb, & n-j-jb+1, cone, a( j, j ),lda, a( j, j+jb ), lda ) end if end do else ! compute the cholesky factorization a = l*l**h. do j = 1, n, nb ! update and factorize the current diagonal block and test ! for non-positive-definiteness. jb = min( nb, n-j+1 ) call stdlib${ii}$_cherk( 'LOWER', 'NO TRANSPOSE', jb, j-1, -one,a( j, 1_${ik}$ ), lda, one,& a( j, j ), lda ) call stdlib${ii}$_cpotrf2( 'LOWER', jb, a( j, j ), lda, info ) if( info/=0 )go to 30 if( j+jb<=n ) then ! compute the current block column. call stdlib${ii}$_cgemm( 'NO TRANSPOSE', 'CONJUGATE TRANSPOSE',n-j-jb+1, jb, j-1,& -cone, a( j+jb, 1_${ik}$ ),lda, a( j, 1_${ik}$ ), lda, cone, a( j+jb, j ),lda ) call stdlib${ii}$_ctrsm( 'RIGHT', 'LOWER', 'CONJUGATE TRANSPOSE','NON-UNIT', n-j-& jb+1, jb, cone, a( j, j ),lda, a( j+jb, j ), lda ) end if end do end if end if go to 40 30 continue info = info + j - 1_${ik}$ 40 continue return end subroutine stdlib${ii}$_cpotrf pure module subroutine stdlib${ii}$_zpotrf( uplo, n, a, lda, info ) !! ZPOTRF computes the Cholesky factorization of a complex Hermitian !! positive definite matrix A. !! The factorization has the form !! A = U**H * U, if UPLO = 'U', or !! A = L * L**H, if UPLO = 'L', !! where U is an upper triangular matrix and L is lower triangular. !! This is the block version of the algorithm, calling Level 3 BLAS. ! -- 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) :: uplo integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, n ! Array Arguments complex(dp), intent(inout) :: a(lda,*) ! ===================================================================== ! Local Scalars logical(lk) :: upper integer(${ik}$) :: j, jb, nb ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ upper = stdlib_lsame( uplo, 'U' ) if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -4_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'ZPOTRF', -info ) return end if ! quick return if possible if( n==0 )return ! determine the block size for this environment. nb = stdlib${ii}$_ilaenv( 1_${ik}$, 'ZPOTRF', uplo, n, -1_${ik}$, -1_${ik}$, -1_${ik}$ ) if( nb<=1_${ik}$ .or. nb>=n ) then ! use unblocked code. call stdlib${ii}$_zpotrf2( uplo, n, a, lda, info ) else ! use blocked code. if( upper ) then ! compute the cholesky factorization a = u**h *u. do j = 1, n, nb ! update and factorize the current diagonal block and test ! for non-positive-definiteness. jb = min( nb, n-j+1 ) call stdlib${ii}$_zherk( 'UPPER', 'CONJUGATE TRANSPOSE', jb, j-1,-one, a( 1_${ik}$, j ), & lda, one, a( j, j ), lda ) call stdlib${ii}$_zpotrf2( 'UPPER', jb, a( j, j ), lda, info ) if( info/=0 )go to 30 if( j+jb<=n ) then ! compute the current block row. call stdlib${ii}$_zgemm( 'CONJUGATE TRANSPOSE', 'NO TRANSPOSE', jb,n-j-jb+1, j-1,& -cone, a( 1_${ik}$, j ), lda,a( 1_${ik}$, j+jb ), lda, cone, a( j, j+jb ),lda ) call stdlib${ii}$_ztrsm( 'LEFT', 'UPPER', 'CONJUGATE TRANSPOSE','NON-UNIT', jb, & n-j-jb+1, cone, a( j, j ),lda, a( j, j+jb ), lda ) end if end do else ! compute the cholesky factorization a = l*l**h. do j = 1, n, nb ! update and factorize the current diagonal block and test ! for non-positive-definiteness. jb = min( nb, n-j+1 ) call stdlib${ii}$_zherk( 'LOWER', 'NO TRANSPOSE', jb, j-1, -one,a( j, 1_${ik}$ ), lda, one,& a( j, j ), lda ) call stdlib${ii}$_zpotrf2( 'LOWER', jb, a( j, j ), lda, info ) if( info/=0 )go to 30 if( j+jb<=n ) then ! compute the current block column. call stdlib${ii}$_zgemm( 'NO TRANSPOSE', 'CONJUGATE TRANSPOSE',n-j-jb+1, jb, j-1,& -cone, a( j+jb, 1_${ik}$ ),lda, a( j, 1_${ik}$ ), lda, cone, a( j+jb, j ),lda ) call stdlib${ii}$_ztrsm( 'RIGHT', 'LOWER', 'CONJUGATE TRANSPOSE','NON-UNIT', n-j-& jb+1, jb, cone, a( j, j ),lda, a( j+jb, j ), lda ) end if end do end if end if go to 40 30 continue info = info + j - 1_${ik}$ 40 continue return end subroutine stdlib${ii}$_zpotrf #:for ck,ct,ci in CMPLX_KINDS_TYPES #:if not ck in ["sp","dp"] pure module subroutine stdlib${ii}$_${ci}$potrf( uplo, n, a, lda, info ) !! ZPOTRF: computes the Cholesky factorization of a complex Hermitian !! positive definite matrix A. !! The factorization has the form !! A = U**H * U, if UPLO = 'U', or !! A = L * L**H, if UPLO = 'L', !! where U is an upper triangular matrix and L is lower triangular. !! This is the block version of the algorithm, calling Level 3 BLAS. ! -- 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) :: uplo integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, n ! Array Arguments complex(${ck}$), intent(inout) :: a(lda,*) ! ===================================================================== ! Local Scalars logical(lk) :: upper integer(${ik}$) :: j, jb, nb ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ upper = stdlib_lsame( uplo, 'U' ) if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -4_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'ZPOTRF', -info ) return end if ! quick return if possible if( n==0 )return ! determine the block size for this environment. nb = stdlib${ii}$_ilaenv( 1_${ik}$, 'ZPOTRF', uplo, n, -1_${ik}$, -1_${ik}$, -1_${ik}$ ) if( nb<=1_${ik}$ .or. nb>=n ) then ! use unblocked code. call stdlib${ii}$_${ci}$potrf2( uplo, n, a, lda, info ) else ! use blocked code. if( upper ) then ! compute the cholesky factorization a = u**h *u. do j = 1, n, nb ! update and factorize the current diagonal block and test ! for non-positive-definiteness. jb = min( nb, n-j+1 ) call stdlib${ii}$_${ci}$herk( 'UPPER', 'CONJUGATE TRANSPOSE', jb, j-1,-one, a( 1_${ik}$, j ), & lda, one, a( j, j ), lda ) call stdlib${ii}$_${ci}$potrf2( 'UPPER', jb, a( j, j ), lda, info ) if( info/=0 )go to 30 if( j+jb<=n ) then ! compute the current block row. call stdlib${ii}$_${ci}$gemm( 'CONJUGATE TRANSPOSE', 'NO TRANSPOSE', jb,n-j-jb+1, j-1,& -cone, a( 1_${ik}$, j ), lda,a( 1_${ik}$, j+jb ), lda, cone, a( j, j+jb ),lda ) call stdlib${ii}$_${ci}$trsm( 'LEFT', 'UPPER', 'CONJUGATE TRANSPOSE','NON-UNIT', jb, & n-j-jb+1, cone, a( j, j ),lda, a( j, j+jb ), lda ) end if end do else ! compute the cholesky factorization a = l*l**h. do j = 1, n, nb ! update and factorize the current diagonal block and test ! for non-positive-definiteness. jb = min( nb, n-j+1 ) call stdlib${ii}$_${ci}$herk( 'LOWER', 'NO TRANSPOSE', jb, j-1, -one,a( j, 1_${ik}$ ), lda, one,& a( j, j ), lda ) call stdlib${ii}$_${ci}$potrf2( 'LOWER', jb, a( j, j ), lda, info ) if( info/=0 )go to 30 if( j+jb<=n ) then ! compute the current block column. call stdlib${ii}$_${ci}$gemm( 'NO TRANSPOSE', 'CONJUGATE TRANSPOSE',n-j-jb+1, jb, j-1,& -cone, a( j+jb, 1_${ik}$ ),lda, a( j, 1_${ik}$ ), lda, cone, a( j+jb, j ),lda ) call stdlib${ii}$_${ci}$trsm( 'RIGHT', 'LOWER', 'CONJUGATE TRANSPOSE','NON-UNIT', n-j-& jb+1, jb, cone, a( j, j ),lda, a( j+jb, j ), lda ) end if end do end if end if go to 40 30 continue info = info + j - 1_${ik}$ 40 continue return end subroutine stdlib${ii}$_${ci}$potrf #:endif #:endfor pure recursive module subroutine stdlib${ii}$_spotrf2( uplo, n, a, lda, info ) !! SPOTRF2 computes the Cholesky factorization of a real symmetric !! positive definite matrix A using the recursive algorithm. !! The factorization has the form !! A = U**T * U, if UPLO = 'U', or !! A = L * L**T, if UPLO = 'L', !! where U is an upper triangular matrix and L is lower triangular. !! This is the recursive version of the algorithm. It divides !! the matrix into four submatrices: !! [ A11 | A12 ] where A11 is n1 by n1 and A22 is n2 by n2 !! A = [ -----|----- ] with n1 = n/2 !! [ A21 | A22 ] n2 = n-n1 !! The subroutine calls itself to factor A11. Update and scale A21 !! or A12, update A22 then call itself to factor A22. ! -- 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) :: uplo integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, n ! Array Arguments real(sp), intent(inout) :: a(lda,*) ! ===================================================================== ! Local Scalars logical(lk) :: upper integer(${ik}$) :: n1, n2, iinfo ! Intrinsic Functions ! Executable Statements ! test the input parameters info = 0_${ik}$ upper = stdlib_lsame( uplo, 'U' ) if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -4_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'SPOTRF2', -info ) return end if ! quick return if possible if( n==0 )return ! n=1 case if( n==1_${ik}$ ) then ! test for non-positive-definiteness if( a( 1_${ik}$, 1_${ik}$ )<=zero.or.stdlib${ii}$_sisnan( a( 1_${ik}$, 1_${ik}$ ) ) ) then info = 1_${ik}$ return end if ! factor a( 1_${ik}$, 1_${ik}$ ) = sqrt( a( 1_${ik}$, 1_${ik}$ ) ) ! use recursive code else n1 = n/2_${ik}$ n2 = n-n1 ! factor a11 call stdlib${ii}$_spotrf2( uplo, n1, a( 1_${ik}$, 1_${ik}$ ), lda, iinfo ) if ( iinfo/=0_${ik}$ ) then info = iinfo return end if ! compute the cholesky factorization a = u**t*u if( upper ) then ! update and scale a12 call stdlib${ii}$_strsm( 'L', 'U', 'T', 'N', n1, n2, one,a( 1_${ik}$, 1_${ik}$ ), lda, a( 1_${ik}$, n1+1 ), & lda ) ! update and factor a22 call stdlib${ii}$_ssyrk( uplo, 'T', n2, n1, -one, a( 1_${ik}$, n1+1 ), lda,one, a( n1+1, n1+1 & ), lda ) call stdlib${ii}$_spotrf2( uplo, n2, a( n1+1, n1+1 ), lda, iinfo ) if ( iinfo/=0_${ik}$ ) then info = iinfo + n1 return end if ! compute the cholesky factorization a = l*l**t else ! update and scale a21 call stdlib${ii}$_strsm( 'R', 'L', 'T', 'N', n2, n1, one,a( 1_${ik}$, 1_${ik}$ ), lda, a( n1+1, 1_${ik}$ ), & lda ) ! update and factor a22 call stdlib${ii}$_ssyrk( uplo, 'N', n2, n1, -one, a( n1+1, 1_${ik}$ ), lda,one, a( n1+1, n1+1 & ), lda ) call stdlib${ii}$_spotrf2( uplo, n2, a( n1+1, n1+1 ), lda, iinfo ) if ( iinfo/=0_${ik}$ ) then info = iinfo + n1 return end if end if end if return end subroutine stdlib${ii}$_spotrf2 pure recursive module subroutine stdlib${ii}$_dpotrf2( uplo, n, a, lda, info ) !! DPOTRF2 computes the Cholesky factorization of a real symmetric !! positive definite matrix A using the recursive algorithm. !! The factorization has the form !! A = U**T * U, if UPLO = 'U', or !! A = L * L**T, if UPLO = 'L', !! where U is an upper triangular matrix and L is lower triangular. !! This is the recursive version of the algorithm. It divides !! the matrix into four submatrices: !! [ A11 | A12 ] where A11 is n1 by n1 and A22 is n2 by n2 !! A = [ -----|----- ] with n1 = n/2 !! [ A21 | A22 ] n2 = n-n1 !! The subroutine calls itself to factor A11. Update and scale A21 !! or A12, update A22 then calls itself to factor A22. ! -- 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) :: uplo integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, n ! Array Arguments real(dp), intent(inout) :: a(lda,*) ! ===================================================================== ! Local Scalars logical(lk) :: upper integer(${ik}$) :: n1, n2, iinfo ! Intrinsic Functions ! Executable Statements ! test the input parameters info = 0_${ik}$ upper = stdlib_lsame( uplo, 'U' ) if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -4_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DPOTRF2', -info ) return end if ! quick return if possible if( n==0 )return ! n=1 case if( n==1_${ik}$ ) then ! test for non-positive-definiteness if( a( 1_${ik}$, 1_${ik}$ )<=zero.or.stdlib${ii}$_disnan( a( 1_${ik}$, 1_${ik}$ ) ) ) then info = 1_${ik}$ return end if ! factor a( 1_${ik}$, 1_${ik}$ ) = sqrt( a( 1_${ik}$, 1_${ik}$ ) ) ! use recursive code else n1 = n/2_${ik}$ n2 = n-n1 ! factor a11 call stdlib${ii}$_dpotrf2( uplo, n1, a( 1_${ik}$, 1_${ik}$ ), lda, iinfo ) if ( iinfo/=0_${ik}$ ) then info = iinfo return end if ! compute the cholesky factorization a = u**t*u if( upper ) then ! update and scale a12 call stdlib${ii}$_dtrsm( 'L', 'U', 'T', 'N', n1, n2, one,a( 1_${ik}$, 1_${ik}$ ), lda, a( 1_${ik}$, n1+1 ), & lda ) ! update and factor a22 call stdlib${ii}$_dsyrk( uplo, 'T', n2, n1, -one, a( 1_${ik}$, n1+1 ), lda,one, a( n1+1, n1+1 & ), lda ) call stdlib${ii}$_dpotrf2( uplo, n2, a( n1+1, n1+1 ), lda, iinfo ) if ( iinfo/=0_${ik}$ ) then info = iinfo + n1 return end if ! compute the cholesky factorization a = l*l**t else ! update and scale a21 call stdlib${ii}$_dtrsm( 'R', 'L', 'T', 'N', n2, n1, one,a( 1_${ik}$, 1_${ik}$ ), lda, a( n1+1, 1_${ik}$ ), & lda ) ! update and factor a22 call stdlib${ii}$_dsyrk( uplo, 'N', n2, n1, -one, a( n1+1, 1_${ik}$ ), lda,one, a( n1+1, n1+1 & ), lda ) call stdlib${ii}$_dpotrf2( uplo, n2, a( n1+1, n1+1 ), lda, iinfo ) if ( iinfo/=0_${ik}$ ) then info = iinfo + n1 return end if end if end if return end subroutine stdlib${ii}$_dpotrf2 #:for rk,rt,ri in REAL_KINDS_TYPES #:if not rk in ["sp","dp"] pure recursive module subroutine stdlib${ii}$_${ri}$potrf2( uplo, n, a, lda, info ) !! DPOTRF2: computes the Cholesky factorization of a real symmetric !! positive definite matrix A using the recursive algorithm. !! The factorization has the form !! A = U**T * U, if UPLO = 'U', or !! A = L * L**T, if UPLO = 'L', !! where U is an upper triangular matrix and L is lower triangular. !! This is the recursive version of the algorithm. It divides !! the matrix into four submatrices: !! [ A11 | A12 ] where A11 is n1 by n1 and A22 is n2 by n2 !! A = [ -----|----- ] with n1 = n/2 !! [ A21 | A22 ] n2 = n-n1 !! The subroutine calls itself to factor A11. Update and scale A21 !! or A12, update A22 then calls itself to factor A22. ! -- 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) :: uplo integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, n ! Array Arguments real(${rk}$), intent(inout) :: a(lda,*) ! ===================================================================== ! Local Scalars logical(lk) :: upper integer(${ik}$) :: n1, n2, iinfo ! Intrinsic Functions ! Executable Statements ! test the input parameters info = 0_${ik}$ upper = stdlib_lsame( uplo, 'U' ) if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -4_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DPOTRF2', -info ) return end if ! quick return if possible if( n==0 )return ! n=1 case if( n==1_${ik}$ ) then ! test for non-positive-definiteness if( a( 1_${ik}$, 1_${ik}$ )<=zero.or.stdlib${ii}$_${ri}$isnan( a( 1_${ik}$, 1_${ik}$ ) ) ) then info = 1_${ik}$ return end if ! factor a( 1_${ik}$, 1_${ik}$ ) = sqrt( a( 1_${ik}$, 1_${ik}$ ) ) ! use recursive code else n1 = n/2_${ik}$ n2 = n-n1 ! factor a11 call stdlib${ii}$_${ri}$potrf2( uplo, n1, a( 1_${ik}$, 1_${ik}$ ), lda, iinfo ) if ( iinfo/=0_${ik}$ ) then info = iinfo return end if ! compute the cholesky factorization a = u**t*u if( upper ) then ! update and scale a12 call stdlib${ii}$_${ri}$trsm( 'L', 'U', 'T', 'N', n1, n2, one,a( 1_${ik}$, 1_${ik}$ ), lda, a( 1_${ik}$, n1+1 ), & lda ) ! update and factor a22 call stdlib${ii}$_${ri}$syrk( uplo, 'T', n2, n1, -one, a( 1_${ik}$, n1+1 ), lda,one, a( n1+1, n1+1 & ), lda ) call stdlib${ii}$_${ri}$potrf2( uplo, n2, a( n1+1, n1+1 ), lda, iinfo ) if ( iinfo/=0_${ik}$ ) then info = iinfo + n1 return end if ! compute the cholesky factorization a = l*l**t else ! update and scale a21 call stdlib${ii}$_${ri}$trsm( 'R', 'L', 'T', 'N', n2, n1, one,a( 1_${ik}$, 1_${ik}$ ), lda, a( n1+1, 1_${ik}$ ), & lda ) ! update and factor a22 call stdlib${ii}$_${ri}$syrk( uplo, 'N', n2, n1, -one, a( n1+1, 1_${ik}$ ), lda,one, a( n1+1, n1+1 & ), lda ) call stdlib${ii}$_${ri}$potrf2( uplo, n2, a( n1+1, n1+1 ), lda, iinfo ) if ( iinfo/=0_${ik}$ ) then info = iinfo + n1 return end if end if end if return end subroutine stdlib${ii}$_${ri}$potrf2 #:endif #:endfor pure recursive module subroutine stdlib${ii}$_cpotrf2( uplo, n, a, lda, info ) !! CPOTRF2 computes the Cholesky factorization of a Hermitian !! positive definite matrix A using the recursive algorithm. !! The factorization has the form !! A = U**H * U, if UPLO = 'U', or !! A = L * L**H, if UPLO = 'L', !! where U is an upper triangular matrix and L is lower triangular. !! This is the recursive version of the algorithm. It divides !! the matrix into four submatrices: !! [ A11 | A12 ] where A11 is n1 by n1 and A22 is n2 by n2 !! A = [ -----|----- ] with n1 = n/2 !! [ A21 | A22 ] n2 = n-n1 !! The subroutine calls itself to factor A11. Update and scale A21 !! or A12, update A22 then calls itself to factor A22. ! -- 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) :: uplo integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, n ! Array Arguments complex(sp), intent(inout) :: a(lda,*) ! ===================================================================== ! Local Scalars logical(lk) :: upper integer(${ik}$) :: n1, n2, iinfo real(sp) :: ajj ! Intrinsic Functions ! Executable Statements ! test the input parameters info = 0_${ik}$ upper = stdlib_lsame( uplo, 'U' ) if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -4_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'CPOTRF2', -info ) return end if ! quick return if possible if( n==0 )return ! n=1 case if( n==1_${ik}$ ) then ! test for non-positive-definiteness ajj = real( a( 1_${ik}$, 1_${ik}$ ),KIND=sp) if( ajj<=zero.or.stdlib${ii}$_sisnan( ajj ) ) then info = 1_${ik}$ return end if ! factor a( 1_${ik}$, 1_${ik}$ ) = sqrt( ajj ) ! use recursive code else n1 = n/2_${ik}$ n2 = n-n1 ! factor a11 call stdlib${ii}$_cpotrf2( uplo, n1, a( 1_${ik}$, 1_${ik}$ ), lda, iinfo ) if ( iinfo/=0_${ik}$ ) then info = iinfo return end if ! compute the cholesky factorization a = u**h*u if( upper ) then ! update and scale a12 call stdlib${ii}$_ctrsm( 'L', 'U', 'C', 'N', n1, n2, cone,a( 1_${ik}$, 1_${ik}$ ), lda, a( 1_${ik}$, n1+1 ),& lda ) ! update and factor a22 call stdlib${ii}$_cherk( uplo, 'C', n2, n1, -one, a( 1_${ik}$, n1+1 ), lda,one, a( n1+1, n1+1 & ), lda ) call stdlib${ii}$_cpotrf2( uplo, n2, a( n1+1, n1+1 ), lda, iinfo ) if ( iinfo/=0_${ik}$ ) then info = iinfo + n1 return end if ! compute the cholesky factorization a = l*l**h else ! update and scale a21 call stdlib${ii}$_ctrsm( 'R', 'L', 'C', 'N', n2, n1, cone,a( 1_${ik}$, 1_${ik}$ ), lda, a( n1+1, 1_${ik}$ ),& lda ) ! update and factor a22 call stdlib${ii}$_cherk( uplo, 'N', n2, n1, -one, a( n1+1, 1_${ik}$ ), lda,one, a( n1+1, n1+1 & ), lda ) call stdlib${ii}$_cpotrf2( uplo, n2, a( n1+1, n1+1 ), lda, iinfo ) if ( iinfo/=0_${ik}$ ) then info = iinfo + n1 return end if end if end if return end subroutine stdlib${ii}$_cpotrf2 pure recursive module subroutine stdlib${ii}$_zpotrf2( uplo, n, a, lda, info ) !! ZPOTRF2 computes the Cholesky factorization of a Hermitian !! positive definite matrix A using the recursive algorithm. !! The factorization has the form !! A = U**H * U, if UPLO = 'U', or !! A = L * L**H, if UPLO = 'L', !! where U is an upper triangular matrix and L is lower triangular. !! This is the recursive version of the algorithm. It divides !! the matrix into four submatrices: !! [ A11 | A12 ] where A11 is n1 by n1 and A22 is n2 by n2 !! A = [ -----|----- ] with n1 = n/2 !! [ A21 | A22 ] n2 = n-n1 !! The subroutine calls itself to factor A11. Update and scale A21 !! or A12, update A22 then call itself to factor A22. ! -- 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) :: uplo integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, n ! Array Arguments complex(dp), intent(inout) :: a(lda,*) ! ===================================================================== ! Local Scalars logical(lk) :: upper integer(${ik}$) :: n1, n2, iinfo real(dp) :: ajj ! Intrinsic Functions ! Executable Statements ! test the input parameters info = 0_${ik}$ upper = stdlib_lsame( uplo, 'U' ) if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -4_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'ZPOTRF2', -info ) return end if ! quick return if possible if( n==0 )return ! n=1 case if( n==1_${ik}$ ) then ! test for non-positive-definiteness ajj = real( a( 1_${ik}$, 1_${ik}$ ),KIND=dp) if( ajj<=zero.or.stdlib${ii}$_disnan( ajj ) ) then info = 1_${ik}$ return end if ! factor a( 1_${ik}$, 1_${ik}$ ) = sqrt( ajj ) ! use recursive code else n1 = n/2_${ik}$ n2 = n-n1 ! factor a11 call stdlib${ii}$_zpotrf2( uplo, n1, a( 1_${ik}$, 1_${ik}$ ), lda, iinfo ) if ( iinfo/=0_${ik}$ ) then info = iinfo return end if ! compute the cholesky factorization a = u**h*u if( upper ) then ! update and scale a12 call stdlib${ii}$_ztrsm( 'L', 'U', 'C', 'N', n1, n2, cone,a( 1_${ik}$, 1_${ik}$ ), lda, a( 1_${ik}$, n1+1 ),& lda ) ! update and factor a22 call stdlib${ii}$_zherk( uplo, 'C', n2, n1, -one, a( 1_${ik}$, n1+1 ), lda,one, a( n1+1, n1+1 & ), lda ) call stdlib${ii}$_zpotrf2( uplo, n2, a( n1+1, n1+1 ), lda, iinfo ) if ( iinfo/=0_${ik}$ ) then info = iinfo + n1 return end if ! compute the cholesky factorization a = l*l**h else ! update and scale a21 call stdlib${ii}$_ztrsm( 'R', 'L', 'C', 'N', n2, n1, cone,a( 1_${ik}$, 1_${ik}$ ), lda, a( n1+1, 1_${ik}$ ),& lda ) ! update and factor a22 call stdlib${ii}$_zherk( uplo, 'N', n2, n1, -one, a( n1+1, 1_${ik}$ ), lda,one, a( n1+1, n1+1 & ), lda ) call stdlib${ii}$_zpotrf2( uplo, n2, a( n1+1, n1+1 ), lda, iinfo ) if ( iinfo/=0_${ik}$ ) then info = iinfo + n1 return end if end if end if return end subroutine stdlib${ii}$_zpotrf2 #:for ck,ct,ci in CMPLX_KINDS_TYPES #:if not ck in ["sp","dp"] pure recursive module subroutine stdlib${ii}$_${ci}$potrf2( uplo, n, a, lda, info ) !! ZPOTRF2: computes the Cholesky factorization of a Hermitian !! positive definite matrix A using the recursive algorithm. !! The factorization has the form !! A = U**H * U, if UPLO = 'U', or !! A = L * L**H, if UPLO = 'L', !! where U is an upper triangular matrix and L is lower triangular. !! This is the recursive version of the algorithm. It divides !! the matrix into four submatrices: !! [ A11 | A12 ] where A11 is n1 by n1 and A22 is n2 by n2 !! A = [ -----|----- ] with n1 = n/2 !! [ A21 | A22 ] n2 = n-n1 !! The subroutine calls itself to factor A11. Update and scale A21 !! or A12, update A22 then call itself to factor A22. ! -- 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) :: uplo integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, n ! Array Arguments complex(${ck}$), intent(inout) :: a(lda,*) ! ===================================================================== ! Local Scalars logical(lk) :: upper integer(${ik}$) :: n1, n2, iinfo real(${ck}$) :: ajj ! Intrinsic Functions ! Executable Statements ! test the input parameters info = 0_${ik}$ upper = stdlib_lsame( uplo, 'U' ) if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -4_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'ZPOTRF2', -info ) return end if ! quick return if possible if( n==0 )return ! n=1 case if( n==1_${ik}$ ) then ! test for non-positive-definiteness ajj = real( a( 1_${ik}$, 1_${ik}$ ),KIND=${ck}$) if( ajj<=zero.or.stdlib${ii}$_${c2ri(ci)}$isnan( ajj ) ) then info = 1_${ik}$ return end if ! factor a( 1_${ik}$, 1_${ik}$ ) = sqrt( ajj ) ! use recursive code else n1 = n/2_${ik}$ n2 = n-n1 ! factor a11 call stdlib${ii}$_${ci}$potrf2( uplo, n1, a( 1_${ik}$, 1_${ik}$ ), lda, iinfo ) if ( iinfo/=0_${ik}$ ) then info = iinfo return end if ! compute the cholesky factorization a = u**h*u if( upper ) then ! update and scale a12 call stdlib${ii}$_${ci}$trsm( 'L', 'U', 'C', 'N', n1, n2, cone,a( 1_${ik}$, 1_${ik}$ ), lda, a( 1_${ik}$, n1+1 ),& lda ) ! update and factor a22 call stdlib${ii}$_${ci}$herk( uplo, 'C', n2, n1, -one, a( 1_${ik}$, n1+1 ), lda,one, a( n1+1, n1+1 & ), lda ) call stdlib${ii}$_${ci}$potrf2( uplo, n2, a( n1+1, n1+1 ), lda, iinfo ) if ( iinfo/=0_${ik}$ ) then info = iinfo + n1 return end if ! compute the cholesky factorization a = l*l**h else ! update and scale a21 call stdlib${ii}$_${ci}$trsm( 'R', 'L', 'C', 'N', n2, n1, cone,a( 1_${ik}$, 1_${ik}$ ), lda, a( n1+1, 1_${ik}$ ),& lda ) ! update and factor a22 call stdlib${ii}$_${ci}$herk( uplo, 'N', n2, n1, -one, a( n1+1, 1_${ik}$ ), lda,one, a( n1+1, n1+1 & ), lda ) call stdlib${ii}$_${ci}$potrf2( uplo, n2, a( n1+1, n1+1 ), lda, iinfo ) if ( iinfo/=0_${ik}$ ) then info = iinfo + n1 return end if end if end if return end subroutine stdlib${ii}$_${ci}$potrf2 #:endif #:endfor pure module subroutine stdlib${ii}$_spotf2( uplo, n, a, lda, info ) !! SPOTF2 computes the Cholesky factorization of a real symmetric !! positive definite matrix A. !! The factorization has the form !! A = U**T * U , if UPLO = 'U', or !! A = L * L**T, if UPLO = 'L', !! where U is an upper triangular matrix and L is lower triangular. !! This is the unblocked version of the algorithm, calling Level 2 BLAS. ! -- 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) :: uplo integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, n ! Array Arguments real(sp), intent(inout) :: a(lda,*) ! ===================================================================== ! Local Scalars logical(lk) :: upper integer(${ik}$) :: j real(sp) :: ajj ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ upper = stdlib_lsame( uplo, 'U' ) if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -4_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'SPOTF2', -info ) return end if ! quick return if possible if( n==0 )return if( upper ) then ! compute the cholesky factorization a = u**t *u. do j = 1, n ! compute u(j,j) and test for non-positive-definiteness. ajj = a( j, j ) - stdlib${ii}$_sdot( j-1, a( 1_${ik}$, j ), 1_${ik}$, a( 1_${ik}$, j ), 1_${ik}$ ) if( ajj<=zero.or.stdlib${ii}$_sisnan( ajj ) ) then a( j, j ) = ajj go to 30 end if ajj = sqrt( ajj ) a( j, j ) = ajj ! compute elements j+1:n of row j. if( j<n ) then call stdlib${ii}$_sgemv( 'TRANSPOSE', j-1, n-j, -one, a( 1_${ik}$, j+1 ),lda, a( 1_${ik}$, j ), 1_${ik}$,& one, a( j, j+1 ), lda ) call stdlib${ii}$_sscal( n-j, one / ajj, a( j, j+1 ), lda ) end if end do else ! compute the cholesky factorization a = l*l**t. do j = 1, n ! compute l(j,j) and test for non-positive-definiteness. ajj = a( j, j ) - stdlib${ii}$_sdot( j-1, a( j, 1_${ik}$ ), lda, a( j, 1_${ik}$ ),lda ) if( ajj<=zero.or.stdlib${ii}$_sisnan( ajj ) ) then a( j, j ) = ajj go to 30 end if ajj = sqrt( ajj ) a( j, j ) = ajj ! compute elements j+1:n of column j. if( j<n ) then call stdlib${ii}$_sgemv( 'NO TRANSPOSE', n-j, j-1, -one, a( j+1, 1_${ik}$ ),lda, a( j, 1_${ik}$ ),& lda, one, a( j+1, j ), 1_${ik}$ ) call stdlib${ii}$_sscal( n-j, one / ajj, a( j+1, j ), 1_${ik}$ ) end if end do end if go to 40 30 continue info = j 40 continue return end subroutine stdlib${ii}$_spotf2 pure module subroutine stdlib${ii}$_dpotf2( uplo, n, a, lda, info ) !! DPOTF2 computes the Cholesky factorization of a real symmetric !! positive definite matrix A. !! The factorization has the form !! A = U**T * U , if UPLO = 'U', or !! A = L * L**T, if UPLO = 'L', !! where U is an upper triangular matrix and L is lower triangular. !! This is the unblocked version of the algorithm, calling Level 2 BLAS. ! -- 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) :: uplo integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, n ! Array Arguments real(dp), intent(inout) :: a(lda,*) ! ===================================================================== ! Local Scalars logical(lk) :: upper integer(${ik}$) :: j real(dp) :: ajj ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ upper = stdlib_lsame( uplo, 'U' ) if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -4_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DPOTF2', -info ) return end if ! quick return if possible if( n==0 )return if( upper ) then ! compute the cholesky factorization a = u**t *u. do j = 1, n ! compute u(j,j) and test for non-positive-definiteness. ajj = a( j, j ) - stdlib${ii}$_ddot( j-1, a( 1_${ik}$, j ), 1_${ik}$, a( 1_${ik}$, j ), 1_${ik}$ ) if( ajj<=zero.or.stdlib${ii}$_disnan( ajj ) ) then a( j, j ) = ajj go to 30 end if ajj = sqrt( ajj ) a( j, j ) = ajj ! compute elements j+1:n of row j. if( j<n ) then call stdlib${ii}$_dgemv( 'TRANSPOSE', j-1, n-j, -one, a( 1_${ik}$, j+1 ),lda, a( 1_${ik}$, j ), 1_${ik}$,& one, a( j, j+1 ), lda ) call stdlib${ii}$_dscal( n-j, one / ajj, a( j, j+1 ), lda ) end if end do else ! compute the cholesky factorization a = l*l**t. do j = 1, n ! compute l(j,j) and test for non-positive-definiteness. ajj = a( j, j ) - stdlib${ii}$_ddot( j-1, a( j, 1_${ik}$ ), lda, a( j, 1_${ik}$ ),lda ) if( ajj<=zero.or.stdlib${ii}$_disnan( ajj ) ) then a( j, j ) = ajj go to 30 end if ajj = sqrt( ajj ) a( j, j ) = ajj ! compute elements j+1:n of column j. if( j<n ) then call stdlib${ii}$_dgemv( 'NO TRANSPOSE', n-j, j-1, -one, a( j+1, 1_${ik}$ ),lda, a( j, 1_${ik}$ ),& lda, one, a( j+1, j ), 1_${ik}$ ) call stdlib${ii}$_dscal( n-j, one / ajj, a( j+1, j ), 1_${ik}$ ) end if end do end if go to 40 30 continue info = j 40 continue return end subroutine stdlib${ii}$_dpotf2 #:for rk,rt,ri in REAL_KINDS_TYPES #:if not rk in ["sp","dp"] pure module subroutine stdlib${ii}$_${ri}$potf2( uplo, n, a, lda, info ) !! DPOTF2: computes the Cholesky factorization of a real symmetric !! positive definite matrix A. !! The factorization has the form !! A = U**T * U , if UPLO = 'U', or !! A = L * L**T, if UPLO = 'L', !! where U is an upper triangular matrix and L is lower triangular. !! This is the unblocked version of the algorithm, calling Level 2 BLAS. ! -- 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) :: uplo integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, n ! Array Arguments real(${rk}$), intent(inout) :: a(lda,*) ! ===================================================================== ! Local Scalars logical(lk) :: upper integer(${ik}$) :: j real(${rk}$) :: ajj ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ upper = stdlib_lsame( uplo, 'U' ) if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -4_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DPOTF2', -info ) return end if ! quick return if possible if( n==0 )return if( upper ) then ! compute the cholesky factorization a = u**t *u. do j = 1, n ! compute u(j,j) and test for non-positive-definiteness. ajj = a( j, j ) - stdlib${ii}$_${ri}$dot( j-1, a( 1_${ik}$, j ), 1_${ik}$, a( 1_${ik}$, j ), 1_${ik}$ ) if( ajj<=zero.or.stdlib${ii}$_${ri}$isnan( ajj ) ) then a( j, j ) = ajj go to 30 end if ajj = sqrt( ajj ) a( j, j ) = ajj ! compute elements j+1:n of row j. if( j<n ) then call stdlib${ii}$_${ri}$gemv( 'TRANSPOSE', j-1, n-j, -one, a( 1_${ik}$, j+1 ),lda, a( 1_${ik}$, j ), 1_${ik}$,& one, a( j, j+1 ), lda ) call stdlib${ii}$_${ri}$scal( n-j, one / ajj, a( j, j+1 ), lda ) end if end do else ! compute the cholesky factorization a = l*l**t. do j = 1, n ! compute l(j,j) and test for non-positive-definiteness. ajj = a( j, j ) - stdlib${ii}$_${ri}$dot( j-1, a( j, 1_${ik}$ ), lda, a( j, 1_${ik}$ ),lda ) if( ajj<=zero.or.stdlib${ii}$_${ri}$isnan( ajj ) ) then a( j, j ) = ajj go to 30 end if ajj = sqrt( ajj ) a( j, j ) = ajj ! compute elements j+1:n of column j. if( j<n ) then call stdlib${ii}$_${ri}$gemv( 'NO TRANSPOSE', n-j, j-1, -one, a( j+1, 1_${ik}$ ),lda, a( j, 1_${ik}$ ),& lda, one, a( j+1, j ), 1_${ik}$ ) call stdlib${ii}$_${ri}$scal( n-j, one / ajj, a( j+1, j ), 1_${ik}$ ) end if end do end if go to 40 30 continue info = j 40 continue return end subroutine stdlib${ii}$_${ri}$potf2 #:endif #:endfor pure module subroutine stdlib${ii}$_cpotf2( uplo, n, a, lda, info ) !! CPOTF2 computes the Cholesky factorization of a complex Hermitian !! positive definite matrix A. !! The factorization has the form !! A = U**H * U , if UPLO = 'U', or !! A = L * L**H, if UPLO = 'L', !! where U is an upper triangular matrix and L is lower triangular. !! This is the unblocked version of the algorithm, calling Level 2 BLAS. ! -- 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) :: uplo integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, n ! Array Arguments complex(sp), intent(inout) :: a(lda,*) ! ===================================================================== ! Local Scalars logical(lk) :: upper integer(${ik}$) :: j real(sp) :: ajj ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ upper = stdlib_lsame( uplo, 'U' ) if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -4_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'CPOTF2', -info ) return end if ! quick return if possible if( n==0 )return if( upper ) then ! compute the cholesky factorization a = u**h *u. do j = 1, n ! compute u(j,j) and test for non-positive-definiteness. ajj = real( real( a( j, j ),KIND=sp) - stdlib${ii}$_cdotc( j-1, a( 1_${ik}$, j ), 1_${ik}$,a( 1_${ik}$, j ),& 1_${ik}$ ),KIND=sp) if( ajj<=zero.or.stdlib${ii}$_sisnan( ajj ) ) then a( j, j ) = ajj go to 30 end if ajj = sqrt( ajj ) a( j, j ) = ajj ! compute elements j+1:n of row j. if( j<n ) then call stdlib${ii}$_clacgv( j-1, a( 1_${ik}$, j ), 1_${ik}$ ) call stdlib${ii}$_cgemv( 'TRANSPOSE', j-1, n-j, -cone, a( 1_${ik}$, j+1 ),lda, a( 1_${ik}$, j ), & 1_${ik}$, cone, a( j, j+1 ), lda ) call stdlib${ii}$_clacgv( j-1, a( 1_${ik}$, j ), 1_${ik}$ ) call stdlib${ii}$_csscal( n-j, one / ajj, a( j, j+1 ), lda ) end if end do else ! compute the cholesky factorization a = l*l**h. do j = 1, n ! compute l(j,j) and test for non-positive-definiteness. ajj = real( real( a( j, j ),KIND=sp) - stdlib${ii}$_cdotc( j-1, a( j, 1_${ik}$ ), lda,a( j, 1_${ik}$ & ), lda ),KIND=sp) if( ajj<=zero.or.stdlib${ii}$_sisnan( ajj ) ) then a( j, j ) = ajj go to 30 end if ajj = sqrt( ajj ) a( j, j ) = ajj ! compute elements j+1:n of column j. if( j<n ) then call stdlib${ii}$_clacgv( j-1, a( j, 1_${ik}$ ), lda ) call stdlib${ii}$_cgemv( 'NO TRANSPOSE', n-j, j-1, -cone, a( j+1, 1_${ik}$ ),lda, a( j, 1_${ik}$ )& , lda, cone, a( j+1, j ), 1_${ik}$ ) call stdlib${ii}$_clacgv( j-1, a( j, 1_${ik}$ ), lda ) call stdlib${ii}$_csscal( n-j, one / ajj, a( j+1, j ), 1_${ik}$ ) end if end do end if go to 40 30 continue info = j 40 continue return end subroutine stdlib${ii}$_cpotf2 pure module subroutine stdlib${ii}$_zpotf2( uplo, n, a, lda, info ) !! ZPOTF2 computes the Cholesky factorization of a complex Hermitian !! positive definite matrix A. !! The factorization has the form !! A = U**H * U , if UPLO = 'U', or !! A = L * L**H, if UPLO = 'L', !! where U is an upper triangular matrix and L is lower triangular. !! This is the unblocked version of the algorithm, calling Level 2 BLAS. ! -- 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) :: uplo integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, n ! Array Arguments complex(dp), intent(inout) :: a(lda,*) ! ===================================================================== ! Local Scalars logical(lk) :: upper integer(${ik}$) :: j real(dp) :: ajj ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ upper = stdlib_lsame( uplo, 'U' ) if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -4_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'ZPOTF2', -info ) return end if ! quick return if possible if( n==0 )return if( upper ) then ! compute the cholesky factorization a = u**h *u. do j = 1, n ! compute u(j,j) and test for non-positive-definiteness. ajj = real( a( j, j ),KIND=dp) - real( stdlib${ii}$_zdotc( j-1, a( 1_${ik}$, j ), 1_${ik}$,a( 1_${ik}$, j ),& 1_${ik}$ ),KIND=dp) if( ajj<=zero.or.stdlib${ii}$_disnan( ajj ) ) then a( j, j ) = ajj go to 30 end if ajj = sqrt( ajj ) a( j, j ) = ajj ! compute elements j+1:n of row j. if( j<n ) then call stdlib${ii}$_zlacgv( j-1, a( 1_${ik}$, j ), 1_${ik}$ ) call stdlib${ii}$_zgemv( 'TRANSPOSE', j-1, n-j, -cone, a( 1_${ik}$, j+1 ),lda, a( 1_${ik}$, j ), & 1_${ik}$, cone, a( j, j+1 ), lda ) call stdlib${ii}$_zlacgv( j-1, a( 1_${ik}$, j ), 1_${ik}$ ) call stdlib${ii}$_zdscal( n-j, one / ajj, a( j, j+1 ), lda ) end if end do else ! compute the cholesky factorization a = l*l**h. do j = 1, n ! compute l(j,j) and test for non-positive-definiteness. ajj = real( a( j, j ),KIND=dp) - real( stdlib${ii}$_zdotc( j-1, a( j, 1_${ik}$ ), lda,a( j, 1_${ik}$ & ), lda ),KIND=dp) if( ajj<=zero.or.stdlib${ii}$_disnan( ajj ) ) then a( j, j ) = ajj go to 30 end if ajj = sqrt( ajj ) a( j, j ) = ajj ! compute elements j+1:n of column j. if( j<n ) then call stdlib${ii}$_zlacgv( j-1, a( j, 1_${ik}$ ), lda ) call stdlib${ii}$_zgemv( 'NO TRANSPOSE', n-j, j-1, -cone, a( j+1, 1_${ik}$ ),lda, a( j, 1_${ik}$ )& , lda, cone, a( j+1, j ), 1_${ik}$ ) call stdlib${ii}$_zlacgv( j-1, a( j, 1_${ik}$ ), lda ) call stdlib${ii}$_zdscal( n-j, one / ajj, a( j+1, j ), 1_${ik}$ ) end if end do end if go to 40 30 continue info = j 40 continue return end subroutine stdlib${ii}$_zpotf2 #:for ck,ct,ci in CMPLX_KINDS_TYPES #:if not ck in ["sp","dp"] pure module subroutine stdlib${ii}$_${ci}$potf2( uplo, n, a, lda, info ) !! ZPOTF2: computes the Cholesky factorization of a complex Hermitian !! positive definite matrix A. !! The factorization has the form !! A = U**H * U , if UPLO = 'U', or !! A = L * L**H, if UPLO = 'L', !! where U is an upper triangular matrix and L is lower triangular. !! This is the unblocked version of the algorithm, calling Level 2 BLAS. ! -- 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) :: uplo integer(${ik}$), intent(out) :: info integer(${ik}$), intent(in) :: lda, n ! Array Arguments complex(${ck}$), intent(inout) :: a(lda,*) ! ===================================================================== ! Local Scalars logical(lk) :: upper integer(${ik}$) :: j real(${ck}$) :: ajj ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ upper = stdlib_lsame( uplo, 'U' ) if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -4_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'ZPOTF2', -info ) return end if ! quick return if possible if( n==0 )return if( upper ) then ! compute the cholesky factorization a = u**h *u. do j = 1, n ! compute u(j,j) and test for non-positive-definiteness. ajj = real( a( j, j ),KIND=${ck}$) - real( stdlib${ii}$_${ci}$dotc( j-1, a( 1_${ik}$, j ), 1_${ik}$,a( 1_${ik}$, j ),& 1_${ik}$ ),KIND=${ck}$) if( ajj<=zero.or.stdlib${ii}$_${c2ri(ci)}$isnan( ajj ) ) then a( j, j ) = ajj go to 30 end if ajj = sqrt( ajj ) a( j, j ) = ajj ! compute elements j+1:n of row j. if( j<n ) then call stdlib${ii}$_${ci}$lacgv( j-1, a( 1_${ik}$, j ), 1_${ik}$ ) call stdlib${ii}$_${ci}$gemv( 'TRANSPOSE', j-1, n-j, -cone, a( 1_${ik}$, j+1 ),lda, a( 1_${ik}$, j ), & 1_${ik}$, cone, a( j, j+1 ), lda ) call stdlib${ii}$_${ci}$lacgv( j-1, a( 1_${ik}$, j ), 1_${ik}$ ) call stdlib${ii}$_${ci}$dscal( n-j, one / ajj, a( j, j+1 ), lda ) end if end do else ! compute the cholesky factorization a = l*l**h. do j = 1, n ! compute l(j,j) and test for non-positive-definiteness. ajj = real( a( j, j ),KIND=${ck}$) - real( stdlib${ii}$_${ci}$dotc( j-1, a( j, 1_${ik}$ ), lda,a( j, 1_${ik}$ & ), lda ),KIND=${ck}$) if( ajj<=zero.or.stdlib${ii}$_${c2ri(ci)}$isnan( ajj ) ) then a( j, j ) = ajj go to 30 end if ajj = sqrt( ajj ) a( j, j ) = ajj ! compute elements j+1:n of column j. if( j<n ) then call stdlib${ii}$_${ci}$lacgv( j-1, a( j, 1_${ik}$ ), lda ) call stdlib${ii}$_${ci}$gemv( 'NO TRANSPOSE', n-j, j-1, -cone, a( j+1, 1_${ik}$ ),lda, a( j, 1_${ik}$ )& , lda, cone, a( j+1, j ), 1_${ik}$ ) call stdlib${ii}$_${ci}$lacgv( j-1, a( j, 1_${ik}$ ), lda ) call stdlib${ii}$_${ci}$dscal( n-j, one / ajj, a( j+1, j ), 1_${ik}$ ) end if end do end if go to 40 30 continue info = j 40 continue return end subroutine stdlib${ii}$_${ci}$potf2 #:endif #:endfor pure module subroutine stdlib${ii}$_spstrf( uplo, n, a, lda, piv, rank, tol, work, info ) !! SPSTRF computes the Cholesky factorization with complete !! pivoting of a real symmetric positive semidefinite matrix A. !! The factorization has the form !! P**T * A * P = U**T * U , if UPLO = 'U', !! P**T * A * P = L * L**T, if UPLO = 'L', !! where U is an upper triangular matrix and L is lower triangular, and !! P is stored as vector PIV. !! This algorithm does not attempt to check that A is positive !! semidefinite. This version of the algorithm calls level 3 BLAS. ! -- 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 real(sp), intent(in) :: tol integer(${ik}$), intent(out) :: info, rank integer(${ik}$), intent(in) :: lda, n character, intent(in) :: uplo ! Array Arguments real(sp), intent(inout) :: a(lda,*) real(sp), intent(out) :: work(2_${ik}$*n) integer(${ik}$), intent(out) :: piv(n) ! ===================================================================== ! Local Scalars real(sp) :: ajj, sstop, stemp integer(${ik}$) :: i, itemp, j, jb, k, nb, pvt logical(lk) :: upper ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ upper = stdlib_lsame( uplo, 'U' ) if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -4_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'SPSTRF', -info ) return end if ! quick return if possible if( n==0 )return ! get block size nb = stdlib${ii}$_ilaenv( 1_${ik}$, 'SPOTRF', uplo, n, -1_${ik}$, -1_${ik}$, -1_${ik}$ ) if( nb<=1_${ik}$ .or. nb>=n ) then ! use unblocked code call stdlib${ii}$_spstf2( uplo, n, a( 1_${ik}$, 1_${ik}$ ), lda, piv, rank, tol, work,info ) go to 200 else ! initialize piv do i = 1, n piv( i ) = i end do ! compute stopping value pvt = 1_${ik}$ ajj = a( pvt, pvt ) do i = 2, n if( a( i, i )>ajj ) then pvt = i ajj = a( pvt, pvt ) end if end do if( ajj<=zero.or.stdlib${ii}$_sisnan( ajj ) ) then rank = 0_${ik}$ info = 1_${ik}$ go to 200 end if ! compute stopping value if not supplied if( tol<zero ) then sstop = n * stdlib${ii}$_slamch( 'EPSILON' ) * ajj else sstop = tol end if if( upper ) then ! compute the cholesky factorization p**t * a * p = u**t * u loop_140: do k = 1, n, nb ! account for last block not being nb wide jb = min( nb, n-k+1 ) ! set relevant part of first half of work to zero, ! holds dot products do i = k, n work( i ) = 0_${ik}$ end do loop_130: do j = k, k + jb - 1 ! find pivot, test for exit, else swap rows and columns ! update dot products, compute possible pivots which are ! stored in the second half of work do i = j, n if( j>k ) then work( i ) = work( i ) + a( j-1, i )**2_${ik}$ end if work( n+i ) = a( i, i ) - work( i ) end do if( j>1_${ik}$ ) then itemp = maxloc( work( (n+j):(2_${ik}$*n) ), 1_${ik}$ ) pvt = itemp + j - 1_${ik}$ ajj = work( n+pvt ) if( ajj<=sstop.or.stdlib${ii}$_sisnan( ajj ) ) then a( j, j ) = ajj go to 190 end if end if if( j/=pvt ) then ! pivot ok, so can now swap pivot rows and columns a( pvt, pvt ) = a( j, j ) call stdlib${ii}$_sswap( j-1, a( 1_${ik}$, j ), 1_${ik}$, a( 1_${ik}$, pvt ), 1_${ik}$ ) if( pvt<n )call stdlib${ii}$_sswap( n-pvt, a( j, pvt+1 ), lda,a( pvt, pvt+1 ),& lda ) call stdlib${ii}$_sswap( pvt-j-1, a( j, j+1 ), lda,a( j+1, pvt ), 1_${ik}$ ) ! swap dot products and piv stemp = work( j ) work( j ) = work( pvt ) work( pvt ) = stemp itemp = piv( pvt ) piv( pvt ) = piv( j ) piv( j ) = itemp end if ajj = sqrt( ajj ) a( j, j ) = ajj ! compute elements j+1:n of row j. if( j<n ) then call stdlib${ii}$_sgemv( 'TRANS', j-k, n-j, -one, a( k, j+1 ),lda, a( k, j ), & 1_${ik}$, one, a( j, j+1 ),lda ) call stdlib${ii}$_sscal( n-j, one / ajj, a( j, j+1 ), lda ) end if end do loop_130 ! update trailing matrix, j already incremented if( k+jb<=n ) then call stdlib${ii}$_ssyrk( 'UPPER', 'TRANS', n-j+1, jb, -one,a( k, j ), lda, one, & a( j, j ), lda ) end if end do loop_140 else ! compute the cholesky factorization p**t * a * p = l * l**t loop_180: do k = 1, n, nb ! account for last block not being nb wide jb = min( nb, n-k+1 ) ! set relevant part of first half of work to zero, ! holds dot products do i = k, n work( i ) = 0_${ik}$ end do loop_170: do j = k, k + jb - 1 ! find pivot, test for exit, else swap rows and columns ! update dot products, compute possible pivots which are ! stored in the second half of work do i = j, n if( j>k ) then work( i ) = work( i ) + a( i, j-1 )**2_${ik}$ end if work( n+i ) = a( i, i ) - work( i ) end do if( j>1_${ik}$ ) then itemp = maxloc( work( (n+j):(2_${ik}$*n) ), 1_${ik}$ ) pvt = itemp + j - 1_${ik}$ ajj = work( n+pvt ) if( ajj<=sstop.or.stdlib${ii}$_sisnan( ajj ) ) then a( j, j ) = ajj go to 190 end if end if if( j/=pvt ) then ! pivot ok, so can now swap pivot rows and columns a( pvt, pvt ) = a( j, j ) call stdlib${ii}$_sswap( j-1, a( j, 1_${ik}$ ), lda, a( pvt, 1_${ik}$ ), lda ) if( pvt<n )call stdlib${ii}$_sswap( n-pvt, a( pvt+1, j ), 1_${ik}$,a( pvt+1, pvt ), & 1_${ik}$ ) call stdlib${ii}$_sswap( pvt-j-1, a( j+1, j ), 1_${ik}$, a( pvt, j+1 ),lda ) ! swap dot products and piv stemp = work( j ) work( j ) = work( pvt ) work( pvt ) = stemp itemp = piv( pvt ) piv( pvt ) = piv( j ) piv( j ) = itemp end if ajj = sqrt( ajj ) a( j, j ) = ajj ! compute elements j+1:n of column j. if( j<n ) then call stdlib${ii}$_sgemv( 'NO TRANS', n-j, j-k, -one,a( j+1, k ), lda, a( j, k & ), lda, one,a( j+1, j ), 1_${ik}$ ) call stdlib${ii}$_sscal( n-j, one / ajj, a( j+1, j ), 1_${ik}$ ) end if end do loop_170 ! update trailing matrix, j already incremented if( k+jb<=n ) then call stdlib${ii}$_ssyrk( 'LOWER', 'NO TRANS', n-j+1, jb, -one,a( j, k ), lda, & one, a( j, j ), lda ) end if end do loop_180 end if end if ! ran to completion, a has full rank rank = n go to 200 190 continue ! rank is the number of steps completed. set info = 1 to signal ! that the factorization cannot be used to solve a system. rank = j - 1_${ik}$ info = 1_${ik}$ 200 continue return end subroutine stdlib${ii}$_spstrf pure module subroutine stdlib${ii}$_dpstrf( uplo, n, a, lda, piv, rank, tol, work, info ) !! DPSTRF computes the Cholesky factorization with complete !! pivoting of a real symmetric positive semidefinite matrix A. !! The factorization has the form !! P**T * A * P = U**T * U , if UPLO = 'U', !! P**T * A * P = L * L**T, if UPLO = 'L', !! where U is an upper triangular matrix and L is lower triangular, and !! P is stored as vector PIV. !! This algorithm does not attempt to check that A is positive !! semidefinite. This version of the algorithm calls level 3 BLAS. ! -- 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 real(dp), intent(in) :: tol integer(${ik}$), intent(out) :: info, rank integer(${ik}$), intent(in) :: lda, n character, intent(in) :: uplo ! Array Arguments real(dp), intent(inout) :: a(lda,*) real(dp), intent(out) :: work(2_${ik}$*n) integer(${ik}$), intent(out) :: piv(n) ! ===================================================================== ! Local Scalars real(dp) :: ajj, dstop, dtemp integer(${ik}$) :: i, itemp, j, jb, k, nb, pvt logical(lk) :: upper ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ upper = stdlib_lsame( uplo, 'U' ) if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -4_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DPSTRF', -info ) return end if ! quick return if possible if( n==0 )return ! get block size nb = stdlib${ii}$_ilaenv( 1_${ik}$, 'DPOTRF', uplo, n, -1_${ik}$, -1_${ik}$, -1_${ik}$ ) if( nb<=1_${ik}$ .or. nb>=n ) then ! use unblocked code call stdlib${ii}$_dpstf2( uplo, n, a( 1_${ik}$, 1_${ik}$ ), lda, piv, rank, tol, work,info ) go to 200 else ! initialize piv do i = 1, n piv( i ) = i end do ! compute stopping value pvt = 1_${ik}$ ajj = a( pvt, pvt ) do i = 2, n if( a( i, i )>ajj ) then pvt = i ajj = a( pvt, pvt ) end if end do if( ajj<=zero.or.stdlib${ii}$_disnan( ajj ) ) then rank = 0_${ik}$ info = 1_${ik}$ go to 200 end if ! compute stopping value if not supplied if( tol<zero ) then dstop = n * stdlib${ii}$_dlamch( 'EPSILON' ) * ajj else dstop = tol end if if( upper ) then ! compute the cholesky factorization p**t * a * p = u**t * u loop_140: do k = 1, n, nb ! account for last block not being nb wide jb = min( nb, n-k+1 ) ! set relevant part of first half of work to zero, ! holds dot products do i = k, n work( i ) = 0_${ik}$ end do loop_130: do j = k, k + jb - 1 ! find pivot, test for exit, else swap rows and columns ! update dot products, compute possible pivots which are ! stored in the second half of work do i = j, n if( j>k ) then work( i ) = work( i ) + a( j-1, i )**2_${ik}$ end if work( n+i ) = a( i, i ) - work( i ) end do if( j>1_${ik}$ ) then itemp = maxloc( work( (n+j):(2_${ik}$*n) ), 1_${ik}$ ) pvt = itemp + j - 1_${ik}$ ajj = work( n+pvt ) if( ajj<=dstop.or.stdlib${ii}$_disnan( ajj ) ) then a( j, j ) = ajj go to 190 end if end if if( j/=pvt ) then ! pivot ok, so can now swap pivot rows and columns a( pvt, pvt ) = a( j, j ) call stdlib${ii}$_dswap( j-1, a( 1_${ik}$, j ), 1_${ik}$, a( 1_${ik}$, pvt ), 1_${ik}$ ) if( pvt<n )call stdlib${ii}$_dswap( n-pvt, a( j, pvt+1 ), lda,a( pvt, pvt+1 ),& lda ) call stdlib${ii}$_dswap( pvt-j-1, a( j, j+1 ), lda,a( j+1, pvt ), 1_${ik}$ ) ! swap dot products and piv dtemp = work( j ) work( j ) = work( pvt ) work( pvt ) = dtemp itemp = piv( pvt ) piv( pvt ) = piv( j ) piv( j ) = itemp end if ajj = sqrt( ajj ) a( j, j ) = ajj ! compute elements j+1:n of row j. if( j<n ) then call stdlib${ii}$_dgemv( 'TRANS', j-k, n-j, -one, a( k, j+1 ),lda, a( k, j ), & 1_${ik}$, one, a( j, j+1 ),lda ) call stdlib${ii}$_dscal( n-j, one / ajj, a( j, j+1 ), lda ) end if end do loop_130 ! update trailing matrix, j already incremented if( k+jb<=n ) then call stdlib${ii}$_dsyrk( 'UPPER', 'TRANS', n-j+1, jb, -one,a( k, j ), lda, one, & a( j, j ), lda ) end if end do loop_140 else ! compute the cholesky factorization p**t * a * p = l * l**t loop_180: do k = 1, n, nb ! account for last block not being nb wide jb = min( nb, n-k+1 ) ! set relevant part of first half of work to zero, ! holds dot products do i = k, n work( i ) = 0_${ik}$ end do loop_170: do j = k, k + jb - 1 ! find pivot, test for exit, else swap rows and columns ! update dot products, compute possible pivots which are ! stored in the second half of work do i = j, n if( j>k ) then work( i ) = work( i ) + a( i, j-1 )**2_${ik}$ end if work( n+i ) = a( i, i ) - work( i ) end do if( j>1_${ik}$ ) then itemp = maxloc( work( (n+j):(2_${ik}$*n) ), 1_${ik}$ ) pvt = itemp + j - 1_${ik}$ ajj = work( n+pvt ) if( ajj<=dstop.or.stdlib${ii}$_disnan( ajj ) ) then a( j, j ) = ajj go to 190 end if end if if( j/=pvt ) then ! pivot ok, so can now swap pivot rows and columns a( pvt, pvt ) = a( j, j ) call stdlib${ii}$_dswap( j-1, a( j, 1_${ik}$ ), lda, a( pvt, 1_${ik}$ ), lda ) if( pvt<n )call stdlib${ii}$_dswap( n-pvt, a( pvt+1, j ), 1_${ik}$,a( pvt+1, pvt ), & 1_${ik}$ ) call stdlib${ii}$_dswap( pvt-j-1, a( j+1, j ), 1_${ik}$, a( pvt, j+1 ),lda ) ! swap dot products and piv dtemp = work( j ) work( j ) = work( pvt ) work( pvt ) = dtemp itemp = piv( pvt ) piv( pvt ) = piv( j ) piv( j ) = itemp end if ajj = sqrt( ajj ) a( j, j ) = ajj ! compute elements j+1:n of column j. if( j<n ) then call stdlib${ii}$_dgemv( 'NO TRANS', n-j, j-k, -one,a( j+1, k ), lda, a( j, k & ), lda, one,a( j+1, j ), 1_${ik}$ ) call stdlib${ii}$_dscal( n-j, one / ajj, a( j+1, j ), 1_${ik}$ ) end if end do loop_170 ! update trailing matrix, j already incremented if( k+jb<=n ) then call stdlib${ii}$_dsyrk( 'LOWER', 'NO TRANS', n-j+1, jb, -one,a( j, k ), lda, & one, a( j, j ), lda ) end if end do loop_180 end if end if ! ran to completion, a has full rank rank = n go to 200 190 continue ! rank is the number of steps completed. set info = 1 to signal ! that the factorization cannot be used to solve a system. rank = j - 1_${ik}$ info = 1_${ik}$ 200 continue return end subroutine stdlib${ii}$_dpstrf #:for rk,rt,ri in REAL_KINDS_TYPES #:if not rk in ["sp","dp"] pure module subroutine stdlib${ii}$_${ri}$pstrf( uplo, n, a, lda, piv, rank, tol, work, info ) !! DPSTRF: computes the Cholesky factorization with complete !! pivoting of a real symmetric positive semidefinite matrix A. !! The factorization has the form !! P**T * A * P = U**T * U , if UPLO = 'U', !! P**T * A * P = L * L**T, if UPLO = 'L', !! where U is an upper triangular matrix and L is lower triangular, and !! P is stored as vector PIV. !! This algorithm does not attempt to check that A is positive !! semidefinite. This version of the algorithm calls level 3 BLAS. ! -- 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 real(${rk}$), intent(in) :: tol integer(${ik}$), intent(out) :: info, rank integer(${ik}$), intent(in) :: lda, n character, intent(in) :: uplo ! Array Arguments real(${rk}$), intent(inout) :: a(lda,*) real(${rk}$), intent(out) :: work(2_${ik}$*n) integer(${ik}$), intent(out) :: piv(n) ! ===================================================================== ! Local Scalars real(${rk}$) :: ajj, dstop, dtemp integer(${ik}$) :: i, itemp, j, jb, k, nb, pvt logical(lk) :: upper ! Intrinsic Functions ! Executable Statements ! test the input parameters. info = 0_${ik}$ upper = stdlib_lsame( uplo, 'U' ) if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then info = -1_${ik}$ else if( n<0_${ik}$ ) then info = -2_${ik}$ else if( lda<max( 1_${ik}$, n ) ) then info = -4_${ik}$ end if if( info/=0_${ik}$ ) then call stdlib${ii}$_xerbla( 'DPSTRF', -info ) return end if ! quick return if possible if( n==0 )return ! get block size nb = stdlib${ii}$_ilaenv( 1_${ik}$, 'DPOTRF', uplo, n, -1_${ik}$, -1_${ik}$, -1_${ik}$ ) if( nb<=1_${ik}$ .or. nb>=n ) then ! use unblocked code call stdlib${ii}$_${ri}$pstf2( uplo, n, a( 1_${ik}$, 1_${ik}$ ), lda, piv, rank, tol, work,info ) go to 200 else ! initialize piv do i = 1, n piv( i ) = i end do ! compute stopping value pvt = 1_${ik}$ ajj = a( pvt, pvt ) do i = 2, n if( a( i, i )>ajj ) then pvt = i ajj = a( pvt, pvt ) end if end do if( ajj<=zero.or.stdlib${ii}$_${ri}$isnan( ajj ) ) then rank = 0_${ik}$ info = 1_${ik}$ go to 200 end if ! compute stopping value if not supplied if( tol<zero ) then dstop = n * stdlib${ii}$_${ri}$lamch( 'EPSILON' ) * ajj else dstop = tol end if if( upper ) then ! compute the cholesky factorization p**t * a * p = u**t * u loop_140: do k = 1, n, nb ! account for last block not being nb wide jb = min( nb, n-k+1 ) ! set relevant part of first half of work to zero, ! holds dot products do i = k, n work( i ) = 0_${ik}$ end do loop_130: do j = k, k + jb - 1 ! find pivot, test for exit, else swap rows and columns ! update dot products, compute possible pivots which are ! stored in the second half of work do i = j, n if( j>k ) then work( i ) = work( i ) + a( j-1, i )**2_${ik}$ end if work( n+i ) = a( i, i ) - work( i ) end do if( j>1_${ik}$ ) then itemp = maxloc( work( (n+j):(2_${ik}$*n) ), 1_${ik}$ ) pvt = itemp + j - 1_${ik}$ ajj = work( n+pvt ) if( ajj<=dstop.or.stdlib${ii}$_${ri}$isnan( ajj ) ) then a( j, j ) = ajj go to 190 end if end if if( j/=pvt ) then ! pivot ok, so can now swap pivot rows and columns a( pvt, pvt ) = a( j, j ) call stdlib${ii}$_${ri}$swap( j-1, a( 1_${ik}$, j ), 1_${ik}$, a( 1_${ik}$, pvt ), 1_${ik}$ ) if( pvt<n )call stdlib${ii}$_${ri}$swap( n-pvt, a( j, pvt+1 ), lda,a( pvt, pvt+1 ),& lda ) call stdlib${ii}$_${ri}$swap( pvt-j-1, a( j, j+1 ), lda,a( j+1, pvt ), 1_${ik}$ ) ! swap dot products and piv dtemp = work( j ) work( j ) = work( pvt ) work( pvt ) = dtemp itemp = piv( pvt ) piv( pvt ) = piv( j ) piv( j ) = itemp end if ajj = sqrt( ajj ) a( j, j ) = ajj ! compute elements j+1:n of row j. if( j<n ) then call stdlib${ii}$_${ri}$gemv( 'TRANS', j-k, n-j, -one, a( k, j+1 ),lda, a( k, j ), & 1_${ik}$, one, a( j, j+1 ),lda ) call stdlib${ii}$_${ri}$scal( n-j, one / ajj, a( j, j+1 ), lda ) end if end do loop_130 ! update trailing matrix, j already incremented if( k+jb<=n ) then call stdlib${ii}$_${ri}$syrk( 'UPPER', 'TRANS', n-j+1, jb, -one,a( k, j ), lda, one, & a( j, j ), lda ) end if end do loop_140 else ! compute the cholesky factorization p**t * a * p = l * l**t loop_180: do k = 1, n, nb ! account for last block not being nb wide jb = min( nb, n-k+1 ) ! set relevant part of first half of work to zero, ! holds dot products do i = k, n work( i ) = 0_${ik}$ end do loop_170: do j = k, k + jb - 1 ! find pivot, test for exit, else swap rows and columns ! update dot products, compute possible pivots which are ! stored in the second half of work do i = j, n if( j>k ) then work( i ) = work( i ) + a( i, j-1 )**2_${ik}$ end if work( n+i ) = a( i, i ) - work( i ) end do if( j>1_${ik}$ ) then itemp = maxloc( work( (n+j):(2_${ik}$*n) ), 1_${ik}$ ) pvt = itemp + j - 1_${ik}$ ajj = work( n+pvt ) if( ajj<=dstop.or.stdlib${ii}$_${ri}$isnan( ajj ) ) then a( j, j ) = ajj go to 190 end if end if if( j/=pvt ) then ! pivot ok, so can now swap pivot rows and columns a( pvt, pvt ) = a( j, j ) call stdlib${ii}$_${ri}$swap( j-1, a( j, 1_${ik}$ ), lda, a( pvt, 1_${ik}$ ), lda ) if( pvt<n )call stdlib${ii}$_${ri}$swap( n-pvt, a( pvt+1, j ), 1_${ik}$,a( pvt+1, pvt ), & 1_${ik}$ ) call stdlib${ii}$_${ri}$swap( pvt-j-1, a( j+1, j ), 1_${ik}$, a( pvt, j+1 ),lda ) ! swap dot products and piv dtemp = work( j ) work( j ) = work( pvt ) work( pvt ) = dtemp itemp = piv( pvt ) piv( pvt ) = piv( j ) piv( j ) = itemp