stdlib_lapack_solve_chol_comp.fypp Source File


Source Code

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