stdlib_lapack_solve_ldl_comp.fypp Source File


Source Code

#:include "common.fypp" 
submodule(stdlib_lapack_solve) stdlib_lapack_solve_ldl_comp
  implicit none


  contains
#:for ik,it,ii in LINALG_INT_KINDS_TYPES

     pure module subroutine stdlib${ii}$_ssycon( uplo, n, a, lda, ipiv, anorm, rcond, work,iwork, info )
     !! SSYCON estimates the reciprocal of the condition number (in the
     !! 1-norm) of a real symmetric matrix A using the factorization
     !! A = U*D*U**T or A = L*D*L**T computed by SSYTRF.
     !! 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(in) :: ipiv(*)
           integer(${ik}$), intent(out) :: iwork(*)
           real(sp), intent(in) :: a(lda,*)
           real(sp), intent(out) :: work(*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: upper
           integer(${ik}$) :: i, kase
           real(sp) :: ainvnm
           ! 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 = -6_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'SSYCON', -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
           ! check that the diagonal matrix d is nonsingular.
           if( upper ) then
              ! upper triangular storage: examine d from bottom to top
              do i = n, 1, -1
                 if( ipiv( i )>0 .and. a( i, i )==zero )return
              end do
           else
              ! lower triangular storage: examine d from top to bottom.
              do i = 1, n
                 if( ipiv( i )>0 .and. a( i, i )==zero )return
              end do
           end if
           ! estimate the 1-norm of the inverse.
           kase = 0_${ik}$
           30 continue
           call stdlib${ii}$_slacn2( n, work( n+1 ), work, iwork, ainvnm, kase, isave )
           if( kase/=0_${ik}$ ) then
              ! multiply by inv(l*d*l**t) or inv(u*d*u**t).
              call stdlib${ii}$_ssytrs( uplo, n, 1_${ik}$, a, lda, ipiv, work, n, info )
              go to 30
           end if
           ! compute the estimate of the reciprocal condition number.
           if( ainvnm/=zero )rcond = ( one / ainvnm ) / anorm
           return
     end subroutine stdlib${ii}$_ssycon

     pure module subroutine stdlib${ii}$_dsycon( uplo, n, a, lda, ipiv, anorm, rcond, work,iwork, info )
     !! DSYCON estimates the reciprocal of the condition number (in the
     !! 1-norm) of a real symmetric matrix A using the factorization
     !! A = U*D*U**T or A = L*D*L**T computed by DSYTRF.
     !! 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(in) :: ipiv(*)
           integer(${ik}$), intent(out) :: iwork(*)
           real(dp), intent(in) :: a(lda,*)
           real(dp), intent(out) :: work(*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: upper
           integer(${ik}$) :: i, kase
           real(dp) :: ainvnm
           ! 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 = -6_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DSYCON', -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
           ! check that the diagonal matrix d is nonsingular.
           if( upper ) then
              ! upper triangular storage: examine d from bottom to top
              do i = n, 1, -1
                 if( ipiv( i )>0 .and. a( i, i )==zero )return
              end do
           else
              ! lower triangular storage: examine d from top to bottom.
              do i = 1, n
                 if( ipiv( i )>0 .and. a( i, i )==zero )return
              end do
           end if
           ! estimate the 1-norm of the inverse.
           kase = 0_${ik}$
           30 continue
           call stdlib${ii}$_dlacn2( n, work( n+1 ), work, iwork, ainvnm, kase, isave )
           if( kase/=0_${ik}$ ) then
              ! multiply by inv(l*d*l**t) or inv(u*d*u**t).
              call stdlib${ii}$_dsytrs( uplo, n, 1_${ik}$, a, lda, ipiv, work, n, info )
              go to 30
           end if
           ! compute the estimate of the reciprocal condition number.
           if( ainvnm/=zero )rcond = ( one / ainvnm ) / anorm
           return
     end subroutine stdlib${ii}$_dsycon

#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
     pure module subroutine stdlib${ii}$_${ri}$sycon( uplo, n, a, lda, ipiv, anorm, rcond, work,iwork, info )
     !! DSYCON: estimates the reciprocal of the condition number (in the
     !! 1-norm) of a real symmetric matrix A using the factorization
     !! A = U*D*U**T or A = L*D*L**T computed by DSYTRF.
     !! 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(in) :: ipiv(*)
           integer(${ik}$), intent(out) :: iwork(*)
           real(${rk}$), intent(in) :: a(lda,*)
           real(${rk}$), intent(out) :: work(*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: upper
           integer(${ik}$) :: i, kase
           real(${rk}$) :: ainvnm
           ! 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 = -6_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DSYCON', -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
           ! check that the diagonal matrix d is nonsingular.
           if( upper ) then
              ! upper triangular storage: examine d from bottom to top
              do i = n, 1, -1
                 if( ipiv( i )>0 .and. a( i, i )==zero )return
              end do
           else
              ! lower triangular storage: examine d from top to bottom.
              do i = 1, n
                 if( ipiv( i )>0 .and. a( i, i )==zero )return
              end do
           end if
           ! estimate the 1-norm of the inverse.
           kase = 0_${ik}$
           30 continue
           call stdlib${ii}$_${ri}$lacn2( n, work( n+1 ), work, iwork, ainvnm, kase, isave )
           if( kase/=0_${ik}$ ) then
              ! multiply by inv(l*d*l**t) or inv(u*d*u**t).
              call stdlib${ii}$_${ri}$sytrs( uplo, n, 1_${ik}$, a, lda, ipiv, work, n, info )
              go to 30
           end if
           ! compute the estimate of the reciprocal condition number.
           if( ainvnm/=zero )rcond = ( one / ainvnm ) / anorm
           return
     end subroutine stdlib${ii}$_${ri}$sycon

#:endif
#:endfor

     pure module subroutine stdlib${ii}$_csycon( uplo, n, a, lda, ipiv, anorm, rcond, work,info )
     !! CSYCON estimates the reciprocal of the condition number (in the
     !! 1-norm) of a complex symmetric matrix A using the factorization
     !! A = U*D*U**T or A = L*D*L**T computed by CSYTRF.
     !! 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(in) :: ipiv(*)
           complex(sp), intent(in) :: a(lda,*)
           complex(sp), intent(out) :: work(*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: upper
           integer(${ik}$) :: i, kase
           real(sp) :: ainvnm
           ! 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 = -6_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'CSYCON', -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
           ! check that the diagonal matrix d is nonsingular.
           if( upper ) then
              ! upper triangular storage: examine d from bottom to top
              do i = n, 1, -1
                 if( ipiv( i )>0 .and. a( i, i )==zero )return
              end do
           else
              ! lower triangular storage: examine d from top to bottom.
              do i = 1, n
                 if( ipiv( i )>0 .and. a( i, i )==zero )return
              end do
           end if
           ! estimate the 1-norm of the inverse.
           kase = 0_${ik}$
           30 continue
           call stdlib${ii}$_clacn2( n, work( n+1 ), work, ainvnm, kase, isave )
           if( kase/=0_${ik}$ ) then
              ! multiply by inv(l*d*l**t) or inv(u*d*u**t).
              call stdlib${ii}$_csytrs( uplo, n, 1_${ik}$, a, lda, ipiv, work, n, info )
              go to 30
           end if
           ! compute the estimate of the reciprocal condition number.
           if( ainvnm/=zero )rcond = ( one / ainvnm ) / anorm
           return
     end subroutine stdlib${ii}$_csycon

     pure module subroutine stdlib${ii}$_zsycon( uplo, n, a, lda, ipiv, anorm, rcond, work,info )
     !! ZSYCON estimates the reciprocal of the condition number (in the
     !! 1-norm) of a complex symmetric matrix A using the factorization
     !! A = U*D*U**T or A = L*D*L**T computed by ZSYTRF.
     !! 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(in) :: ipiv(*)
           complex(dp), intent(in) :: a(lda,*)
           complex(dp), intent(out) :: work(*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: upper
           integer(${ik}$) :: i, kase
           real(dp) :: ainvnm
           ! 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 = -6_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'ZSYCON', -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
           ! check that the diagonal matrix d is nonsingular.
           if( upper ) then
              ! upper triangular storage: examine d from bottom to top
              do i = n, 1, -1
                 if( ipiv( i )>0 .and. a( i, i )==zero )return
              end do
           else
              ! lower triangular storage: examine d from top to bottom.
              do i = 1, n
                 if( ipiv( i )>0 .and. a( i, i )==zero )return
              end do
           end if
           ! estimate the 1-norm of the inverse.
           kase = 0_${ik}$
           30 continue
           call stdlib${ii}$_zlacn2( n, work( n+1 ), work, ainvnm, kase, isave )
           if( kase/=0_${ik}$ ) then
              ! multiply by inv(l*d*l**t) or inv(u*d*u**t).
              call stdlib${ii}$_zsytrs( uplo, n, 1_${ik}$, a, lda, ipiv, work, n, info )
              go to 30
           end if
           ! compute the estimate of the reciprocal condition number.
           if( ainvnm/=zero )rcond = ( one / ainvnm ) / anorm
           return
     end subroutine stdlib${ii}$_zsycon

#:for ck,ct,ci in CMPLX_KINDS_TYPES
#:if not ck in ["sp","dp"]
     pure module subroutine stdlib${ii}$_${ci}$sycon( uplo, n, a, lda, ipiv, anorm, rcond, work,info )
     !! ZSYCON: estimates the reciprocal of the condition number (in the
     !! 1-norm) of a complex symmetric matrix A using the factorization
     !! A = U*D*U**T or A = L*D*L**T computed by ZSYTRF.
     !! 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 
           integer(${ik}$), intent(in) :: ipiv(*)
           complex(${ck}$), intent(in) :: a(lda,*)
           complex(${ck}$), intent(out) :: work(*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: upper
           integer(${ik}$) :: i, kase
           real(${ck}$) :: ainvnm
           ! 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 = -6_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'ZSYCON', -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
           ! check that the diagonal matrix d is nonsingular.
           if( upper ) then
              ! upper triangular storage: examine d from bottom to top
              do i = n, 1, -1
                 if( ipiv( i )>0 .and. a( i, i )==zero )return
              end do
           else
              ! lower triangular storage: examine d from top to bottom.
              do i = 1, n
                 if( ipiv( i )>0 .and. a( i, i )==zero )return
              end do
           end if
           ! estimate the 1-norm of the inverse.
           kase = 0_${ik}$
           30 continue
           call stdlib${ii}$_${ci}$lacn2( n, work( n+1 ), work, ainvnm, kase, isave )
           if( kase/=0_${ik}$ ) then
              ! multiply by inv(l*d*l**t) or inv(u*d*u**t).
              call stdlib${ii}$_${ci}$sytrs( uplo, n, 1_${ik}$, a, lda, ipiv, work, n, info )
              go to 30
           end if
           ! compute the estimate of the reciprocal condition number.
           if( ainvnm/=zero )rcond = ( one / ainvnm ) / anorm
           return
     end subroutine stdlib${ii}$_${ci}$sycon

#:endif
#:endfor



     pure module subroutine stdlib${ii}$_ssytrf( uplo, n, a, lda, ipiv, work, lwork, info )
     !! SSYTRF computes the factorization of a real symmetric matrix A using
     !! the Bunch-Kaufman diagonal pivoting method.  The form of the
     !! factorization is
     !! A = U**T*D*U  or  A = L*D*L**T
     !! where U (or L) is a product of permutation and unit upper (lower)
     !! triangular matrices, and D is symmetric and block diagonal with
     !! 1-by-1 and 2-by-2 diagonal blocks.
     !! This is the blocked 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, lwork, n
           ! Array Arguments 
           integer(${ik}$), intent(out) :: ipiv(*)
           real(sp), intent(inout) :: a(lda,*)
           real(sp), intent(out) :: work(*)
        ! =====================================================================
           ! Local Scalars 
           logical(lk) :: lquery, upper
           integer(${ik}$) :: iinfo, iws, j, k, kb, ldwork, lwkopt, nb, nbmin
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           info = 0_${ik}$
           upper = stdlib_lsame( uplo, 'U' )
           lquery = ( lwork==-1_${ik}$ )
           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( lwork<1_${ik}$ .and. .not.lquery ) then
              info = -7_${ik}$
           end if
           if( info==0_${ik}$ ) then
              ! determine the block size
              nb = stdlib${ii}$_ilaenv( 1_${ik}$, 'SSYTRF', uplo, n, -1_${ik}$, -1_${ik}$, -1_${ik}$ )
              lwkopt = n*nb
              work( 1_${ik}$ ) = lwkopt
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'SSYTRF', -info )
              return
           else if( lquery ) then
              return
           end if
           nbmin = 2_${ik}$
           ldwork = n
           if( nb>1_${ik}$ .and. nb<n ) then
              iws = ldwork*nb
              if( lwork<iws ) then
                 nb = max( lwork / ldwork, 1_${ik}$ )
                 nbmin = max( 2_${ik}$, stdlib${ii}$_ilaenv( 2_${ik}$, 'SSYTRF', uplo, n, -1_${ik}$, -1_${ik}$, -1_${ik}$ ) )
              end if
           else
              iws = 1_${ik}$
           end if
           if( nb<nbmin )nb = n
           if( upper ) then
              ! factorize a as u**t*d*u using the upper triangle of a
              ! k is the main loop index, decreasing from n to 1 in steps of
              ! kb, where kb is the number of columns factorized by stdlib${ii}$_slasyf;
              ! kb is either nb or nb-1, or k for the last block
              k = n
              10 continue
              ! if k < 1, exit from loop
              if( k<1 )go to 40
              if( k>nb ) then
                 ! factorize columns k-kb+1:k of a and use blocked code to
                 ! update columns 1:k-kb
                 call stdlib${ii}$_slasyf( uplo, k, nb, kb, a, lda, ipiv, work, ldwork,iinfo )
              else
                 ! use unblocked code to factorize columns 1:k of a
                 call stdlib${ii}$_ssytf2( uplo, k, a, lda, ipiv, iinfo )
                 kb = k
              end if
              ! set info on the first occurrence of a zero pivot
              if( info==0_${ik}$ .and. iinfo>0_${ik}$ )info = iinfo
              ! decrease k and return to the start of the main loop
              k = k - kb
              go to 10
           else
              ! factorize a as l*d*l**t using the lower triangle of a
              ! k is the main loop index, increasing from 1 to n in steps of
              ! kb, where kb is the number of columns factorized by stdlib${ii}$_slasyf;
              ! kb is either nb or nb-1, or n-k+1 for the last block
              k = 1_${ik}$
              20 continue
              ! if k > n, exit from loop
              if( k>n )go to 40
              if( k<=n-nb ) then
                 ! factorize columns k:k+kb-1 of a and use blocked code to
                 ! update columns k+kb:n
                 call stdlib${ii}$_slasyf( uplo, n-k+1, nb, kb, a( k, k ), lda, ipiv( k ),work, ldwork, &
                           iinfo )
              else
                 ! use unblocked code to factorize columns k:n of a
                 call stdlib${ii}$_ssytf2( uplo, n-k+1, a( k, k ), lda, ipiv( k ), iinfo )
                 kb = n - k + 1_${ik}$
              end if
              ! set info on the first occurrence of a zero pivot
              if( info==0_${ik}$ .and. iinfo>0_${ik}$ )info = iinfo + k - 1_${ik}$
              ! adjust ipiv
              do j = k, k + kb - 1
                 if( ipiv( j )>0_${ik}$ ) then
                    ipiv( j ) = ipiv( j ) + k - 1_${ik}$
                 else
                    ipiv( j ) = ipiv( j ) - k + 1_${ik}$
                 end if
              end do
              ! increase k and return to the start of the main loop
              k = k + kb
              go to 20
           end if
           40 continue
           work( 1_${ik}$ ) = lwkopt
           return
     end subroutine stdlib${ii}$_ssytrf

     pure module subroutine stdlib${ii}$_dsytrf( uplo, n, a, lda, ipiv, work, lwork, info )
     !! DSYTRF computes the factorization of a real symmetric matrix A using
     !! the Bunch-Kaufman diagonal pivoting method.  The form of the
     !! factorization is
     !! A = U**T*D*U  or  A = L*D*L**T
     !! where U (or L) is a product of permutation and unit upper (lower)
     !! triangular matrices, and D is symmetric and block diagonal with
     !! 1-by-1 and 2-by-2 diagonal blocks.
     !! This is the blocked 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, lwork, n
           ! Array Arguments 
           integer(${ik}$), intent(out) :: ipiv(*)
           real(dp), intent(inout) :: a(lda,*)
           real(dp), intent(out) :: work(*)
        ! =====================================================================
           ! Local Scalars 
           logical(lk) :: lquery, upper
           integer(${ik}$) :: iinfo, iws, j, k, kb, ldwork, lwkopt, nb, nbmin
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           info = 0_${ik}$
           upper = stdlib_lsame( uplo, 'U' )
           lquery = ( lwork==-1_${ik}$ )
           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( lwork<1_${ik}$ .and. .not.lquery ) then
              info = -7_${ik}$
           end if
           if( info==0_${ik}$ ) then
              ! determine the block size
              nb = stdlib${ii}$_ilaenv( 1_${ik}$, 'DSYTRF', uplo, n, -1_${ik}$, -1_${ik}$, -1_${ik}$ )
              lwkopt = n*nb
              work( 1_${ik}$ ) = lwkopt
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DSYTRF', -info )
              return
           else if( lquery ) then
              return
           end if
           nbmin = 2_${ik}$
           ldwork = n
           if( nb>1_${ik}$ .and. nb<n ) then
              iws = ldwork*nb
              if( lwork<iws ) then
                 nb = max( lwork / ldwork, 1_${ik}$ )
                 nbmin = max( 2_${ik}$, stdlib${ii}$_ilaenv( 2_${ik}$, 'DSYTRF', uplo, n, -1_${ik}$, -1_${ik}$, -1_${ik}$ ) )
              end if
           else
              iws = 1_${ik}$
           end if
           if( nb<nbmin )nb = n
           if( upper ) then
              ! factorize a as u**t*d*u using the upper triangle of a
              ! k is the main loop index, decreasing from n to 1 in steps of
              ! kb, where kb is the number of columns factorized by stdlib${ii}$_dlasyf;
              ! kb is either nb or nb-1, or k for the last block
              k = n
              10 continue
              ! if k < 1, exit from loop
              if( k<1 )go to 40
              if( k>nb ) then
                 ! factorize columns k-kb+1:k of a and use blocked code to
                 ! update columns 1:k-kb
                 call stdlib${ii}$_dlasyf( uplo, k, nb, kb, a, lda, ipiv, work, ldwork,iinfo )
              else
                 ! use unblocked code to factorize columns 1:k of a
                 call stdlib${ii}$_dsytf2( uplo, k, a, lda, ipiv, iinfo )
                 kb = k
              end if
              ! set info on the first occurrence of a zero pivot
              if( info==0_${ik}$ .and. iinfo>0_${ik}$ )info = iinfo
              ! decrease k and return to the start of the main loop
              k = k - kb
              go to 10
           else
              ! factorize a as l*d*l**t using the lower triangle of a
              ! k is the main loop index, increasing from 1 to n in steps of
              ! kb, where kb is the number of columns factorized by stdlib${ii}$_dlasyf;
              ! kb is either nb or nb-1, or n-k+1 for the last block
              k = 1_${ik}$
              20 continue
              ! if k > n, exit from loop
              if( k>n )go to 40
              if( k<=n-nb ) then
                 ! factorize columns k:k+kb-1 of a and use blocked code to
                 ! update columns k+kb:n
                 call stdlib${ii}$_dlasyf( uplo, n-k+1, nb, kb, a( k, k ), lda, ipiv( k ),work, ldwork, &
                           iinfo )
              else
                 ! use unblocked code to factorize columns k:n of a
                 call stdlib${ii}$_dsytf2( uplo, n-k+1, a( k, k ), lda, ipiv( k ), iinfo )
                 kb = n - k + 1_${ik}$
              end if
              ! set info on the first occurrence of a zero pivot
              if( info==0_${ik}$ .and. iinfo>0_${ik}$ )info = iinfo + k - 1_${ik}$
              ! adjust ipiv
              do j = k, k + kb - 1
                 if( ipiv( j )>0_${ik}$ ) then
                    ipiv( j ) = ipiv( j ) + k - 1_${ik}$
                 else
                    ipiv( j ) = ipiv( j ) - k + 1_${ik}$
                 end if
              end do
              ! increase k and return to the start of the main loop
              k = k + kb
              go to 20
           end if
           40 continue
           work( 1_${ik}$ ) = lwkopt
           return
     end subroutine stdlib${ii}$_dsytrf

#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
     pure module subroutine stdlib${ii}$_${ri}$sytrf( uplo, n, a, lda, ipiv, work, lwork, info )
     !! DSYTRF: computes the factorization of a real symmetric matrix A using
     !! the Bunch-Kaufman diagonal pivoting method.  The form of the
     !! factorization is
     !! A = U**T*D*U  or  A = L*D*L**T
     !! where U (or L) is a product of permutation and unit upper (lower)
     !! triangular matrices, and D is symmetric and block diagonal with
     !! 1-by-1 and 2-by-2 diagonal blocks.
     !! This is the blocked 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, lwork, n
           ! Array Arguments 
           integer(${ik}$), intent(out) :: ipiv(*)
           real(${rk}$), intent(inout) :: a(lda,*)
           real(${rk}$), intent(out) :: work(*)
        ! =====================================================================
           ! Local Scalars 
           logical(lk) :: lquery, upper
           integer(${ik}$) :: iinfo, iws, j, k, kb, ldwork, lwkopt, nb, nbmin
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           info = 0_${ik}$
           upper = stdlib_lsame( uplo, 'U' )
           lquery = ( lwork==-1_${ik}$ )
           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( lwork<1_${ik}$ .and. .not.lquery ) then
              info = -7_${ik}$
           end if
           if( info==0_${ik}$ ) then
              ! determine the block size
              nb = stdlib${ii}$_ilaenv( 1_${ik}$, 'DSYTRF', uplo, n, -1_${ik}$, -1_${ik}$, -1_${ik}$ )
              lwkopt = n*nb
              work( 1_${ik}$ ) = lwkopt
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DSYTRF', -info )
              return
           else if( lquery ) then
              return
           end if
           nbmin = 2_${ik}$
           ldwork = n
           if( nb>1_${ik}$ .and. nb<n ) then
              iws = ldwork*nb
              if( lwork<iws ) then
                 nb = max( lwork / ldwork, 1_${ik}$ )
                 nbmin = max( 2_${ik}$, stdlib${ii}$_ilaenv( 2_${ik}$, 'DSYTRF', uplo, n, -1_${ik}$, -1_${ik}$, -1_${ik}$ ) )
              end if
           else
              iws = 1_${ik}$
           end if
           if( nb<nbmin )nb = n
           if( upper ) then
              ! factorize a as u**t*d*u using the upper triangle of a
              ! k is the main loop index, decreasing from n to 1 in steps of
              ! kb, where kb is the number of columns factorized by stdlib${ii}$_${ri}$lasyf;
              ! kb is either nb or nb-1, or k for the last block
              k = n
              10 continue
              ! if k < 1, exit from loop
              if( k<1 )go to 40
              if( k>nb ) then
                 ! factorize columns k-kb+1:k of a and use blocked code to
                 ! update columns 1:k-kb
                 call stdlib${ii}$_${ri}$lasyf( uplo, k, nb, kb, a, lda, ipiv, work, ldwork,iinfo )
              else
                 ! use unblocked code to factorize columns 1:k of a
                 call stdlib${ii}$_${ri}$sytf2( uplo, k, a, lda, ipiv, iinfo )
                 kb = k
              end if
              ! set info on the first occurrence of a zero pivot
              if( info==0_${ik}$ .and. iinfo>0_${ik}$ )info = iinfo
              ! decrease k and return to the start of the main loop
              k = k - kb
              go to 10
           else
              ! factorize a as l*d*l**t using the lower triangle of a
              ! k is the main loop index, increasing from 1 to n in steps of
              ! kb, where kb is the number of columns factorized by stdlib${ii}$_${ri}$lasyf;
              ! kb is either nb or nb-1, or n-k+1 for the last block
              k = 1_${ik}$
              20 continue
              ! if k > n, exit from loop
              if( k>n )go to 40
              if( k<=n-nb ) then
                 ! factorize columns k:k+kb-1 of a and use blocked code to
                 ! update columns k+kb:n
                 call stdlib${ii}$_${ri}$lasyf( uplo, n-k+1, nb, kb, a( k, k ), lda, ipiv( k ),work, ldwork, &
                           iinfo )
              else
                 ! use unblocked code to factorize columns k:n of a
                 call stdlib${ii}$_${ri}$sytf2( uplo, n-k+1, a( k, k ), lda, ipiv( k ), iinfo )
                 kb = n - k + 1_${ik}$
              end if
              ! set info on the first occurrence of a zero pivot
              if( info==0_${ik}$ .and. iinfo>0_${ik}$ )info = iinfo + k - 1_${ik}$
              ! adjust ipiv
              do j = k, k + kb - 1
                 if( ipiv( j )>0_${ik}$ ) then
                    ipiv( j ) = ipiv( j ) + k - 1_${ik}$
                 else
                    ipiv( j ) = ipiv( j ) - k + 1_${ik}$
                 end if
              end do
              ! increase k and return to the start of the main loop
              k = k + kb
              go to 20
           end if
           40 continue
           work( 1_${ik}$ ) = lwkopt
           return
     end subroutine stdlib${ii}$_${ri}$sytrf

#:endif
#:endfor

     pure module subroutine stdlib${ii}$_csytrf( uplo, n, a, lda, ipiv, work, lwork, info )
     !! CSYTRF computes the factorization of a complex symmetric matrix A
     !! using the Bunch-Kaufman diagonal pivoting method.  The form of the
     !! factorization is
     !! A = U*D*U**T  or  A = L*D*L**T
     !! where U (or L) is a product of permutation and unit upper (lower)
     !! triangular matrices, and D is symmetric and block diagonal with
     !! 1-by-1 and 2-by-2 diagonal blocks.
     !! This is the blocked 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, lwork, n
           ! Array Arguments 
           integer(${ik}$), intent(out) :: ipiv(*)
           complex(sp), intent(inout) :: a(lda,*)
           complex(sp), intent(out) :: work(*)
        ! =====================================================================
           ! Local Scalars 
           logical(lk) :: lquery, upper
           integer(${ik}$) :: iinfo, iws, j, k, kb, ldwork, lwkopt, nb, nbmin
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           info = 0_${ik}$
           upper = stdlib_lsame( uplo, 'U' )
           lquery = ( lwork==-1_${ik}$ )
           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( lwork<1_${ik}$ .and. .not.lquery ) then
              info = -7_${ik}$
           end if
           if( info==0_${ik}$ ) then
              ! determine the block size
              nb = stdlib${ii}$_ilaenv( 1_${ik}$, 'CSYTRF', uplo, n, -1_${ik}$, -1_${ik}$, -1_${ik}$ )
              lwkopt = n*nb
              work( 1_${ik}$ ) = lwkopt
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'CSYTRF', -info )
              return
           else if( lquery ) then
              return
           end if
           nbmin = 2_${ik}$
           ldwork = n
           if( nb>1_${ik}$ .and. nb<n ) then
              iws = ldwork*nb
              if( lwork<iws ) then
                 nb = max( lwork / ldwork, 1_${ik}$ )
                 nbmin = max( 2_${ik}$, stdlib${ii}$_ilaenv( 2_${ik}$, 'CSYTRF', uplo, n, -1_${ik}$, -1_${ik}$, -1_${ik}$ ) )
              end if
           else
              iws = 1_${ik}$
           end if
           if( nb<nbmin )nb = n
           if( upper ) then
              ! factorize a as u*d*u**t using the upper triangle of a
              ! k is the main loop index, decreasing from n to 1 in steps of
              ! kb, where kb is the number of columns factorized by stdlib${ii}$_clasyf;
              ! kb is either nb or nb-1, or k for the last block
              k = n
              10 continue
              ! if k < 1, exit from loop
              if( k<1 )go to 40
              if( k>nb ) then
                 ! factorize columns k-kb+1:k of a and use blocked code to
                 ! update columns 1:k-kb
                 call stdlib${ii}$_clasyf( uplo, k, nb, kb, a, lda, ipiv, work, n, iinfo )
              else
                 ! use unblocked code to factorize columns 1:k of a
                 call stdlib${ii}$_csytf2( uplo, k, a, lda, ipiv, iinfo )
                 kb = k
              end if
              ! set info on the first occurrence of a zero pivot
              if( info==0_${ik}$ .and. iinfo>0_${ik}$ )info = iinfo
              ! decrease k and return to the start of the main loop
              k = k - kb
              go to 10
           else
              ! factorize a as l*d*l**t using the lower triangle of a
              ! k is the main loop index, increasing from 1 to n in steps of
              ! kb, where kb is the number of columns factorized by stdlib${ii}$_clasyf;
              ! kb is either nb or nb-1, or n-k+1 for the last block
              k = 1_${ik}$
              20 continue
              ! if k > n, exit from loop
              if( k>n )go to 40
              if( k<=n-nb ) then
                 ! factorize columns k:k+kb-1 of a and use blocked code to
                 ! update columns k+kb:n
                 call stdlib${ii}$_clasyf( uplo, n-k+1, nb, kb, a( k, k ), lda, ipiv( k ),work, n, &
                           iinfo )
              else
                 ! use unblocked code to factorize columns k:n of a
                 call stdlib${ii}$_csytf2( uplo, n-k+1, a( k, k ), lda, ipiv( k ), iinfo )
                 kb = n - k + 1_${ik}$
              end if
              ! set info on the first occurrence of a zero pivot
              if( info==0_${ik}$ .and. iinfo>0_${ik}$ )info = iinfo + k - 1_${ik}$
              ! adjust ipiv
              do j = k, k + kb - 1
                 if( ipiv( j )>0_${ik}$ ) then
                    ipiv( j ) = ipiv( j ) + k - 1_${ik}$
                 else
                    ipiv( j ) = ipiv( j ) - k + 1_${ik}$
                 end if
              end do
              ! increase k and return to the start of the main loop
              k = k + kb
              go to 20
           end if
           40 continue
           work( 1_${ik}$ ) = lwkopt
           return
     end subroutine stdlib${ii}$_csytrf

     pure module subroutine stdlib${ii}$_zsytrf( uplo, n, a, lda, ipiv, work, lwork, info )
     !! ZSYTRF computes the factorization of a complex symmetric matrix A
     !! using the Bunch-Kaufman diagonal pivoting method.  The form of the
     !! factorization is
     !! A = U*D*U**T  or  A = L*D*L**T
     !! where U (or L) is a product of permutation and unit upper (lower)
     !! triangular matrices, and D is symmetric and block diagonal with
     !! 1-by-1 and 2-by-2 diagonal blocks.
     !! This is the blocked 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, lwork, n
           ! Array Arguments 
           integer(${ik}$), intent(out) :: ipiv(*)
           complex(dp), intent(inout) :: a(lda,*)
           complex(dp), intent(out) :: work(*)
        ! =====================================================================
           ! Local Scalars 
           logical(lk) :: lquery, upper
           integer(${ik}$) :: iinfo, iws, j, k, kb, ldwork, lwkopt, nb, nbmin
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           info = 0_${ik}$
           upper = stdlib_lsame( uplo, 'U' )
           lquery = ( lwork==-1_${ik}$ )
           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( lwork<1_${ik}$ .and. .not.lquery ) then
              info = -7_${ik}$
           end if
           if( info==0_${ik}$ ) then
              ! determine the block size
              nb = stdlib${ii}$_ilaenv( 1_${ik}$, 'ZSYTRF', uplo, n, -1_${ik}$, -1_${ik}$, -1_${ik}$ )
              lwkopt = n*nb
              work( 1_${ik}$ ) = lwkopt
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'ZSYTRF', -info )
              return
           else if( lquery ) then
              return
           end if
           nbmin = 2_${ik}$
           ldwork = n
           if( nb>1_${ik}$ .and. nb<n ) then
              iws = ldwork*nb
              if( lwork<iws ) then
                 nb = max( lwork / ldwork, 1_${ik}$ )
                 nbmin = max( 2_${ik}$, stdlib${ii}$_ilaenv( 2_${ik}$, 'ZSYTRF', uplo, n, -1_${ik}$, -1_${ik}$, -1_${ik}$ ) )
              end if
           else
              iws = 1_${ik}$
           end if
           if( nb<nbmin )nb = n
           if( upper ) then
              ! factorize a as u*d*u**t using the upper triangle of a
              ! k is the main loop index, decreasing from n to 1 in steps of
              ! kb, where kb is the number of columns factorized by stdlib${ii}$_zlasyf;
              ! kb is either nb or nb-1, or k for the last block
              k = n
              10 continue
              ! if k < 1, exit from loop
              if( k<1 )go to 40
              if( k>nb ) then
                 ! factorize columns k-kb+1:k of a and use blocked code to
                 ! update columns 1:k-kb
                 call stdlib${ii}$_zlasyf( uplo, k, nb, kb, a, lda, ipiv, work, n, iinfo )
              else
                 ! use unblocked code to factorize columns 1:k of a
                 call stdlib${ii}$_zsytf2( uplo, k, a, lda, ipiv, iinfo )
                 kb = k
              end if
              ! set info on the first occurrence of a zero pivot
              if( info==0_${ik}$ .and. iinfo>0_${ik}$ )info = iinfo
              ! decrease k and return to the start of the main loop
              k = k - kb
              go to 10
           else
              ! factorize a as l*d*l**t using the lower triangle of a
              ! k is the main loop index, increasing from 1 to n in steps of
              ! kb, where kb is the number of columns factorized by stdlib${ii}$_zlasyf;
              ! kb is either nb or nb-1, or n-k+1 for the last block
              k = 1_${ik}$
              20 continue
              ! if k > n, exit from loop
              if( k>n )go to 40
              if( k<=n-nb ) then
                 ! factorize columns k:k+kb-1 of a and use blocked code to
                 ! update columns k+kb:n
                 call stdlib${ii}$_zlasyf( uplo, n-k+1, nb, kb, a( k, k ), lda, ipiv( k ),work, n, &
                           iinfo )
              else
                 ! use unblocked code to factorize columns k:n of a
                 call stdlib${ii}$_zsytf2( uplo, n-k+1, a( k, k ), lda, ipiv( k ), iinfo )
                 kb = n - k + 1_${ik}$
              end if
              ! set info on the first occurrence of a zero pivot
              if( info==0_${ik}$ .and. iinfo>0_${ik}$ )info = iinfo + k - 1_${ik}$
              ! adjust ipiv
              do j = k, k + kb - 1
                 if( ipiv( j )>0_${ik}$ ) then
                    ipiv( j ) = ipiv( j ) + k - 1_${ik}$
                 else
                    ipiv( j ) = ipiv( j ) - k + 1_${ik}$
                 end if
              end do
              ! increase k and return to the start of the main loop
              k = k + kb
              go to 20
           end if
           40 continue
           work( 1_${ik}$ ) = lwkopt
           return
     end subroutine stdlib${ii}$_zsytrf

#:for ck,ct,ci in CMPLX_KINDS_TYPES
#:if not ck in ["sp","dp"]
     pure module subroutine stdlib${ii}$_${ci}$sytrf( uplo, n, a, lda, ipiv, work, lwork, info )
     !! ZSYTRF: computes the factorization of a complex symmetric matrix A
     !! using the Bunch-Kaufman diagonal pivoting method.  The form of the
     !! factorization is
     !! A = U*D*U**T  or  A = L*D*L**T
     !! where U (or L) is a product of permutation and unit upper (lower)
     !! triangular matrices, and D is symmetric and block diagonal with
     !! 1-by-1 and 2-by-2 diagonal blocks.
     !! This is the blocked 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, lwork, n
           ! Array Arguments 
           integer(${ik}$), intent(out) :: ipiv(*)
           complex(${ck}$), intent(inout) :: a(lda,*)
           complex(${ck}$), intent(out) :: work(*)
        ! =====================================================================
           ! Local Scalars 
           logical(lk) :: lquery, upper
           integer(${ik}$) :: iinfo, iws, j, k, kb, ldwork, lwkopt, nb, nbmin
           ! Intrinsic Functions 
           ! Executable Statements 
           ! test the input parameters.
           info = 0_${ik}$
           upper = stdlib_lsame( uplo, 'U' )
           lquery = ( lwork==-1_${ik}$ )
           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( lwork<1_${ik}$ .and. .not.lquery ) then
              info = -7_${ik}$
           end if
           if( info==0_${ik}$ ) then
              ! determine the block size
              nb = stdlib${ii}$_ilaenv( 1_${ik}$, 'ZSYTRF', uplo, n, -1_${ik}$, -1_${ik}$, -1_${ik}$ )
              lwkopt = n*nb
              work( 1_${ik}$ ) = lwkopt
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'ZSYTRF', -info )
              return
           else if( lquery ) then
              return
           end if
           nbmin = 2_${ik}$
           ldwork = n
           if( nb>1_${ik}$ .and. nb<n ) then
              iws = ldwork*nb
              if( lwork<iws ) then
                 nb = max( lwork / ldwork, 1_${ik}$ )
                 nbmin = max( 2_${ik}$, stdlib${ii}$_ilaenv( 2_${ik}$, 'ZSYTRF', uplo, n, -1_${ik}$, -1_${ik}$, -1_${ik}$ ) )
              end if
           else
              iws = 1_${ik}$
           end if
           if( nb<nbmin )nb = n
           if( upper ) then
              ! factorize a as u*d*u**t using the upper triangle of a
              ! k is the main loop index, decreasing from n to 1 in steps of
              ! kb, where kb is the number of columns factorized by stdlib${ii}$_${ci}$lasyf;
              ! kb is either nb or nb-1, or k for the last block
              k = n
              10 continue
              ! if k < 1, exit from loop
              if( k<1 )go to 40
              if( k>nb ) then
                 ! factorize columns k-kb+1:k of a and use blocked code to
                 ! update columns 1:k-kb
                 call stdlib${ii}$_${ci}$lasyf( uplo, k, nb, kb, a, lda, ipiv, work, n, iinfo )
              else
                 ! use unblocked code to factorize columns 1:k of a
                 call stdlib${ii}$_${ci}$sytf2( uplo, k, a, lda, ipiv, iinfo )
                 kb = k
              end if
              ! set info on the first occurrence of a zero pivot
              if( info==0_${ik}$ .and. iinfo>0_${ik}$ )info = iinfo
              ! decrease k and return to the start of the main loop
              k = k - kb
              go to 10
           else
              ! factorize a as l*d*l**t using the lower triangle of a
              ! k is the main loop index, increasing from 1 to n in steps of
              ! kb, where kb is the number of columns factorized by stdlib${ii}$_${ci}$lasyf;
              ! kb is either nb or nb-1, or n-k+1 for the last block
              k = 1_${ik}$
              20 continue
              ! if k > n, exit from loop
              if( k>n )go to 40
              if( k<=n-nb ) then
                 ! factorize columns k:k+kb-1 of a and use blocked code to
                 ! update columns k+kb:n
                 call stdlib${ii}$_${ci}$lasyf( uplo, n-k+1, nb, kb, a( k, k ), lda, ipiv( k ),work, n, &
                           iinfo )
              else
                 ! use unblocked code to factorize columns k:n of a
                 call stdlib${ii}$_${ci}$sytf2( uplo, n-k+1, a( k, k ), lda, ipiv( k ), iinfo )
                 kb = n - k + 1_${ik}$
              end if
              ! set info on the first occurrence of a zero pivot
              if( info==0_${ik}$ .and. iinfo>0_${ik}$ )info = iinfo + k - 1_${ik}$
              ! adjust ipiv
              do j = k, k + kb - 1
                 if( ipiv( j )>0_${ik}$ ) then
                    ipiv( j ) = ipiv( j ) + k - 1_${ik}$
                 else
                    ipiv( j ) = ipiv( j ) - k + 1_${ik}$
                 end if
              end do
              ! increase k and return to the start of the main loop
              k = k + kb
              go to 20
           end if
           40 continue
           work( 1_${ik}$ ) = lwkopt
           return
     end subroutine stdlib${ii}$_${ci}$sytrf

#:endif
#:endfor



     pure module subroutine stdlib${ii}$_slasyf( uplo, n, nb, kb, a, lda, ipiv, w, ldw, info )
     !! SLASYF computes a partial factorization of a real symmetric matrix A
     !! using the Bunch-Kaufman diagonal pivoting method. The partial
     !! factorization has the form:
     !! A  =  ( I  U12 ) ( A11  0  ) (  I       0    )  if UPLO = 'U', or:
     !! ( 0  U22 ) (  0   D  ) ( U12**T U22**T )
     !! A  =  ( L11  0 ) (  D   0  ) ( L11**T L21**T )  if UPLO = 'L'
     !! ( L21  I ) (  0  A22 ) (  0       I    )
     !! where the order of D is at most NB. The actual order is returned in
     !! the argument KB, and is either NB or NB-1, or N if N <= NB.
     !! SLASYF is an auxiliary routine called by SSYTRF. It uses blocked code
     !! (calling Level 3 BLAS) to update the submatrix A11 (if UPLO = 'U') or
     !! A22 (if UPLO = 'L').
        ! -- 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, kb
           integer(${ik}$), intent(in) :: lda, ldw, n, nb
           ! Array Arguments 
           integer(${ik}$), intent(out) :: ipiv(*)
           real(sp), intent(inout) :: a(lda,*)
           real(sp), intent(out) :: w(ldw,*)
        ! =====================================================================
           ! Parameters 
           real(sp), parameter :: sevten = 17.0e+0_sp
           
           
           ! Local Scalars 
           integer(${ik}$) :: imax, j, jb, jj, jmax, jp, k, kk, kkw, kp, kstep, kw
           real(sp) :: absakk, alpha, colmax, d11, d21, d22, r1, rowmax, t
           ! Intrinsic Functions 
           ! Executable Statements 
           info = 0_${ik}$
           ! initialize alpha for use in choosing pivot block size.
           alpha = ( one+sqrt( sevten ) ) / eight
           if( stdlib_lsame( uplo, 'U' ) ) then
              ! factorize the trailing columns of a using the upper triangle
              ! of a and working backwards, and compute the matrix w = u12*d
              ! for use in updating a11
              ! k is the main loop index, decreasing from n in steps of 1 or 2
              ! kw is the column of w which corresponds to column k of a
              k = n
              10 continue
              kw = nb + k - n
              ! exit from loop
              if( ( k<=n-nb+1 .and. nb<n ) .or. k<1 )go to 30
              ! copy column k of a to column kw of w and update it
              call stdlib${ii}$_scopy( k, a( 1_${ik}$, k ), 1_${ik}$, w( 1_${ik}$, kw ), 1_${ik}$ )
              if( k<n )call stdlib${ii}$_sgemv( 'NO TRANSPOSE', k, n-k, -one, a( 1_${ik}$, k+1 ), lda,w( k, kw+&
                        1_${ik}$ ), ldw, one, w( 1_${ik}$, kw ), 1_${ik}$ )
              kstep = 1_${ik}$
              ! determine rows and columns to be interchanged and whether
              ! a 1-by-1 or 2-by-2 pivot block will be used
              absakk = abs( w( k, kw ) )
              ! imax is the row-index of the largest off-diagonal element in
              ! column k, and colmax is its absolute value.
              ! determine both colmax and imax.
              if( k>1_${ik}$ ) then
                 imax = stdlib${ii}$_isamax( k-1, w( 1_${ik}$, kw ), 1_${ik}$ )
                 colmax = abs( w( imax, kw ) )
              else
                 colmax = zero
              end if
              if( max( absakk, colmax )==zero ) then
                 ! column k is zero or underflow: set info and continue
                 if( info==0_${ik}$ )info = k
                 kp = k
              else
                 if( absakk>=alpha*colmax ) then
                    ! no interchange, use 1-by-1 pivot block
                    kp = k
                 else
                    ! copy column imax to column kw-1 of w and update it
                    call stdlib${ii}$_scopy( imax, a( 1_${ik}$, imax ), 1_${ik}$, w( 1_${ik}$, kw-1 ), 1_${ik}$ )
                    call stdlib${ii}$_scopy( k-imax, a( imax, imax+1 ), lda,w( imax+1, kw-1 ), 1_${ik}$ )
                              
                    if( k<n )call stdlib${ii}$_sgemv( 'NO TRANSPOSE', k, n-k, -one, a( 1_${ik}$, k+1 ),lda, w( &
                              imax, kw+1 ), ldw, one,w( 1_${ik}$, kw-1 ), 1_${ik}$ )
                    ! jmax is the column-index of the largest off-diagonal
                    ! element in row imax, and rowmax is its absolute value
                    jmax = imax + stdlib${ii}$_isamax( k-imax, w( imax+1, kw-1 ), 1_${ik}$ )
                    rowmax = abs( w( jmax, kw-1 ) )
                    if( imax>1_${ik}$ ) then
                       jmax = stdlib${ii}$_isamax( imax-1, w( 1_${ik}$, kw-1 ), 1_${ik}$ )
                       rowmax = max( rowmax, abs( w( jmax, kw-1 ) ) )
                    end if
                    if( absakk>=alpha*colmax*( colmax / rowmax ) ) then
                       ! no interchange, use 1-by-1 pivot block
                       kp = k
                    else if( abs( w( imax, kw-1 ) )>=alpha*rowmax ) then
                       ! interchange rows and columns k and imax, use 1-by-1
                       ! pivot block
                       kp = imax
                       ! copy column kw-1 of w to column kw of w
                       call stdlib${ii}$_scopy( k, w( 1_${ik}$, kw-1 ), 1_${ik}$, w( 1_${ik}$, kw ), 1_${ik}$ )
                    else
                       ! interchange rows and columns k-1 and imax, use 2-by-2
                       ! pivot block
                       kp = imax
                       kstep = 2_${ik}$
                    end if
                 end if
                 ! ============================================================
                 ! kk is the column of a where pivoting step stopped
                 kk = k - kstep + 1_${ik}$
                 ! kkw is the column of w which corresponds to column kk of a
                 kkw = nb + kk - n
                 ! interchange rows and columns kp and kk.
                 ! updated column kp is already stored in column kkw of w.
                 if( kp/=kk ) then
                    ! copy non-updated column kk to column kp of submatrix a
                    ! at step k. no need to copy element into column k
                    ! (or k and k-1 for 2-by-2 pivot) of a, since these columns
                    ! will be later overwritten.
                    a( kp, kp ) = a( kk, kk )
                    call stdlib${ii}$_scopy( kk-1-kp, a( kp+1, kk ), 1_${ik}$, a( kp, kp+1 ),lda )
                    if( kp>1_${ik}$ )call stdlib${ii}$_scopy( kp-1, a( 1_${ik}$, kk ), 1_${ik}$, a( 1_${ik}$, kp ), 1_${ik}$ )
                    ! interchange rows kk and kp in last k+1 to n columns of a
                    ! (columns k (or k and k-1 for 2-by-2 pivot) of a will be
                    ! later overwritten). interchange rows kk and kp
                    ! in last kkw to nb columns of w.
                    if( k<n )call stdlib${ii}$_sswap( n-k, a( kk, k+1 ), lda, a( kp, k+1 ),lda )
                    call stdlib${ii}$_sswap( n-kk+1, w( kk, kkw ), ldw, w( kp, kkw ),ldw )
                 end if
                 if( kstep==1_${ik}$ ) then
                    ! 1-by-1 pivot block d(k): column kw of w now holds
                    ! w(kw) = u(k)*d(k),
                    ! where u(k) is the k-th column of u
                    ! store subdiag. elements of column u(k)
                    ! and 1-by-1 block d(k) in column k of a.
                    ! note: diagonal element u(k,k) is a unit element
                    ! and not stored.
                       ! a(k,k) := d(k,k) = w(k,kw)
                       ! a(1:k-1,k) := u(1:k-1,k) = w(1:k-1,kw)/d(k,k)
                    call stdlib${ii}$_scopy( k, w( 1_${ik}$, kw ), 1_${ik}$, a( 1_${ik}$, k ), 1_${ik}$ )
                    r1 = one / a( k, k )
                    call stdlib${ii}$_sscal( k-1, r1, a( 1_${ik}$, k ), 1_${ik}$ )
                 else
                    ! 2-by-2 pivot block d(k): columns kw and kw-1 of w now hold
                    ! ( w(kw-1) w(kw) ) = ( u(k-1) u(k) )*d(k)
                    ! where u(k) and u(k-1) are the k-th and (k-1)-th columns
                    ! of u
                    ! store u(1:k-2,k-1) and u(1:k-2,k) and 2-by-2
                    ! block d(k-1:k,k-1:k) in columns k-1 and k of a.
                    ! note: 2-by-2 diagonal block u(k-1:k,k-1:k) is a unit
                    ! block and not stored.
                       ! a(k-1:k,k-1:k) := d(k-1:k,k-1:k) = w(k-1:k,kw-1:kw)
                       ! a(1:k-2,k-1:k) := u(1:k-2,k:k-1:k) =
                       ! = w(1:k-2,kw-1:kw) * ( d(k-1:k,k-1:k)**(-1) )
                    if( k>2_${ik}$ ) then
                       ! compose the columns of the inverse of 2-by-2 pivot
                       ! block d in the following way to reduce the number
                       ! of flops when we myltiply panel ( w(kw-1) w(kw) ) by
                       ! this inverse
                       ! d**(-1) = ( d11 d21 )**(-1) =
                                 ! ( d21 d22 )
                       ! = 1/(d11*d22-d21**2) * ( ( d22 ) (-d21 ) ) =
                                              ! ( (-d21 ) ( d11 ) )
                       ! = 1/d21 * 1/((d11/d21)*(d22/d21)-1) *
                         ! * ( ( d22/d21 ) (      -1 ) ) =
                           ! ( (      -1 ) ( d11/d21 ) )
                       ! = 1/d21 * 1/(d22*d11-1) * ( ( d11 ) (  -1 ) ) =
                                                 ! ( ( -1  ) ( d22 ) )
                       ! = 1/d21 * t * ( ( d11 ) (  -1 ) )
                                     ! ( (  -1 ) ( d22 ) )
                       ! = d21 * ( ( d11 ) (  -1 ) )
                               ! ( (  -1 ) ( d22 ) )
                       d21 = w( k-1, kw )
                       d11 = w( k, kw ) / d21
                       d22 = w( k-1, kw-1 ) / d21
                       t = one / ( d11*d22-one )
                       d21 = t / d21
                       ! update elements in columns a(k-1) and a(k) as
                       ! dot products of rows of ( w(kw-1) w(kw) ) and columns
                       ! of d**(-1)
                       do j = 1, k - 2
                          a( j, k-1 ) = d21*( d11*w( j, kw-1 )-w( j, kw ) )
                          a( j, k ) = d21*( d22*w( j, kw )-w( j, kw-1 ) )
                       end do
                    end if
                    ! copy d(k) to a
                    a( k-1, k-1 ) = w( k-1, kw-1 )
                    a( k-1, k ) = w( k-1, kw )
                    a( k, k ) = w( k, kw )
                 end if
              end if
              ! store details of the interchanges in ipiv
              if( kstep==1_${ik}$ ) then
                 ipiv( k ) = kp
              else
                 ipiv( k ) = -kp
                 ipiv( k-1 ) = -kp
              end if
              ! decrease k and return to the start of the main loop
              k = k - kstep
              go to 10
              30 continue
              ! update the upper triangle of a11 (= a(1:k,1:k)) as
              ! a11 := a11 - u12*d*u12**t = a11 - u12*w**t
              ! computing blocks of nb columns at a time
              do j = ( ( k-1 ) / nb )*nb + 1, 1, -nb
                 jb = min( nb, k-j+1 )
                 ! update the upper triangle of the diagonal block
                 do jj = j, j + jb - 1
                    call stdlib${ii}$_sgemv( 'NO TRANSPOSE', jj-j+1, n-k, -one,a( j, k+1 ), lda, w( jj, &
                              kw+1 ), ldw, one,a( j, jj ), 1_${ik}$ )
                 end do
                 ! update the rectangular superdiagonal block
                 call stdlib${ii}$_sgemm( 'NO TRANSPOSE', 'TRANSPOSE', j-1, jb, n-k, -one,a( 1_${ik}$, k+1 ), &
                           lda, w( j, kw+1 ), ldw, one,a( 1_${ik}$, j ), lda )
              end do
              ! put u12 in standard form by partially undoing the interchanges
              ! in columns k+1:n looping backwards from k+1 to n
              j = k + 1_${ik}$
              60 continue
                 ! undo the interchanges (if any) of rows jj and jp at each
                 ! step j
                 ! (here, j is a diagonal index)
                 jj = j
                 jp = ipiv( j )
                 if( jp<0_${ik}$ ) then
                    jp = -jp
                    ! (here, j is a diagonal index)
                    j = j + 1_${ik}$
                 end if
                 ! (note: here, j is used to determine row length. length n-j+1
                 ! of the rows to swap back doesn't include diagonal element)
                 j = j + 1_${ik}$
                 if( jp/=jj .and. j<=n )call stdlib${ii}$_sswap( n-j+1, a( jp, j ), lda, a( jj, j ), &
                           lda )
              if( j<n )go to 60
              ! set kb to the number of columns factorized
              kb = n - k
           else
              ! factorize the leading columns of a using the lower triangle
              ! of a and working forwards, and compute the matrix w = l21*d
              ! for use in updating a22
              ! k is the main loop index, increasing from 1 in steps of 1 or 2
              k = 1_${ik}$
              70 continue
              ! exit from loop
              if( ( k>=nb .and. nb<n ) .or. k>n )go to 90
              ! copy column k of a to column k of w and update it
              call stdlib${ii}$_scopy( n-k+1, a( k, k ), 1_${ik}$, w( k, k ), 1_${ik}$ )
              call stdlib${ii}$_sgemv( 'NO TRANSPOSE', n-k+1, k-1, -one, a( k, 1_${ik}$ ), lda,w( k, 1_${ik}$ ), ldw, &
                        one, w( k, k ), 1_${ik}$ )
              kstep = 1_${ik}$
              ! determine rows and columns to be interchanged and whether
              ! a 1-by-1 or 2-by-2 pivot block will be used
              absakk = abs( w( k, k ) )
              ! imax is the row-index of the largest off-diagonal element in
              ! column k, and colmax is its absolute value.
              ! determine both colmax and imax.
              if( k<n ) then
                 imax = k + stdlib${ii}$_isamax( n-k, w( k+1, k ), 1_${ik}$ )
                 colmax = abs( w( imax, k ) )
              else
                 colmax = zero
              end if
              if( max( absakk, colmax )==zero ) then
                 ! column k is zero or underflow: set info and continue
                 if( info==0_${ik}$ )info = k
                 kp = k
              else
                 if( absakk>=alpha*colmax ) then
                    ! no interchange, use 1-by-1 pivot block
                    kp = k
                 else
                    ! copy column imax to column k+1 of w and update it
                    call stdlib${ii}$_scopy( imax-k, a( imax, k ), lda, w( k, k+1 ), 1_${ik}$ )
                    call stdlib${ii}$_scopy( n-imax+1, a( imax, imax ), 1_${ik}$, w( imax, k+1 ),1_${ik}$ )
                    call stdlib${ii}$_sgemv( 'NO TRANSPOSE', n-k+1, k-1, -one, a( k, 1_${ik}$ ),lda, w( imax, &
                              1_${ik}$ ), ldw, one, w( k, k+1 ), 1_${ik}$ )
                    ! jmax is the column-index of the largest off-diagonal
                    ! element in row imax, and rowmax is its absolute value
                    jmax = k - 1_${ik}$ + stdlib${ii}$_isamax( imax-k, w( k, k+1 ), 1_${ik}$ )
                    rowmax = abs( w( jmax, k+1 ) )
                    if( imax<n ) then
                       jmax = imax + stdlib${ii}$_isamax( n-imax, w( imax+1, k+1 ), 1_${ik}$ )
                       rowmax = max( rowmax, abs( w( jmax, k+1 ) ) )
                    end if
                    if( absakk>=alpha*colmax*( colmax / rowmax ) ) then
                       ! no interchange, use 1-by-1 pivot block
                       kp = k
                    else if( abs( w( imax, k+1 ) )>=alpha*rowmax ) then
                       ! interchange rows and columns k and imax, use 1-by-1
                       ! pivot block
                       kp = imax
                       ! copy column k+1 of w to column k of w
                       call stdlib${ii}$_scopy( n-k+1, w( k, k+1 ), 1_${ik}$, w( k, k ), 1_${ik}$ )
                    else
                       ! interchange rows and columns k+1 and imax, use 2-by-2
                       ! pivot block
                       kp = imax
                       kstep = 2_${ik}$
                    end if
                 end if
                 ! ============================================================
                 ! kk is the column of a where pivoting step stopped
                 kk = k + kstep - 1_${ik}$
                 ! interchange rows and columns kp and kk.
                 ! updated column kp is already stored in column kk of w.
                 if( kp/=kk ) then
                    ! copy non-updated column kk to column kp of submatrix a
                    ! at step k. no need to copy element into column k
                    ! (or k and k+1 for 2-by-2 pivot) of a, since these columns
                    ! will be later overwritten.
                    a( kp, kp ) = a( kk, kk )
                    call stdlib${ii}$_scopy( kp-kk-1, a( kk+1, kk ), 1_${ik}$, a( kp, kk+1 ),lda )
                    if( kp<n )call stdlib${ii}$_scopy( n-kp, a( kp+1, kk ), 1_${ik}$, a( kp+1, kp ), 1_${ik}$ )
                              
                    ! interchange rows kk and kp in first k-1 columns of a
                    ! (columns k (or k and k+1 for 2-by-2 pivot) of a will be
                    ! later overwritten). interchange rows kk and kp
                    ! in first kk columns of w.
                    if( k>1_${ik}$ )call stdlib${ii}$_sswap( k-1, a( kk, 1_${ik}$ ), lda, a( kp, 1_${ik}$ ), lda )
                    call stdlib${ii}$_sswap( kk, w( kk, 1_${ik}$ ), ldw, w( kp, 1_${ik}$ ), ldw )
                 end if
                 if( kstep==1_${ik}$ ) then
                    ! 1-by-1 pivot block d(k): column k of w now holds
                    ! w(k) = l(k)*d(k),
                    ! where l(k) is the k-th column of l
                    ! store subdiag. elements of column l(k)
                    ! and 1-by-1 block d(k) in column k of a.
                    ! (note: diagonal element l(k,k) is a unit element
                    ! and not stored)
                       ! a(k,k) := d(k,k) = w(k,k)
                       ! a(k+1:n,k) := l(k+1:n,k) = w(k+1:n,k)/d(k,k)
                    call stdlib${ii}$_scopy( n-k+1, w( k, k ), 1_${ik}$, a( k, k ), 1_${ik}$ )
                    if( k<n ) then
                       r1 = one / a( k, k )
                       call stdlib${ii}$_sscal( n-k, r1, a( k+1, k ), 1_${ik}$ )
                    end if
                 else
                    ! 2-by-2 pivot block d(k): columns k and k+1 of w now hold
                    ! ( w(k) w(k+1) ) = ( l(k) l(k+1) )*d(k)
                    ! where l(k) and l(k+1) are the k-th and (k+1)-th columns
                    ! of l
                    ! store l(k+2:n,k) and l(k+2:n,k+1) and 2-by-2
                    ! block d(k:k+1,k:k+1) in columns k and k+1 of a.
                    ! (note: 2-by-2 diagonal block l(k:k+1,k:k+1) is a unit
                    ! block and not stored)
                       ! a(k:k+1,k:k+1) := d(k:k+1,k:k+1) = w(k:k+1,k:k+1)
                       ! a(k+2:n,k:k+1) := l(k+2:n,k:k+1) =
                       ! = w(k+2:n,k:k+1) * ( d(k:k+1,k:k+1)**(-1) )
                    if( k<n-1 ) then
                       ! compose the columns of the inverse of 2-by-2 pivot
                       ! block d in the following way to reduce the number
                       ! of flops when we myltiply panel ( w(k) w(k+1) ) by
                       ! this inverse
                       ! d**(-1) = ( d11 d21 )**(-1) =
                                 ! ( d21 d22 )
                       ! = 1/(d11*d22-d21**2) * ( ( d22 ) (-d21 ) ) =
                                              ! ( (-d21 ) ( d11 ) )
                       ! = 1/d21 * 1/((d11/d21)*(d22/d21)-1) *
                         ! * ( ( d22/d21 ) (      -1 ) ) =
                           ! ( (      -1 ) ( d11/d21 ) )
                       ! = 1/d21 * 1/(d22*d11-1) * ( ( d11 ) (  -1 ) ) =
                                                 ! ( ( -1  ) ( d22 ) )
                       ! = 1/d21 * t * ( ( d11 ) (  -1 ) )
                                     ! ( (  -1 ) ( d22 ) )
                       ! = d21 * ( ( d11 ) (  -1 ) )
                               ! ( (  -1 ) ( d22 ) )
                       d21 = w( k+1, k )
                       d11 = w( k+1, k+1 ) / d21
                       d22 = w( k, k ) / d21
                       t = one / ( d11*d22-one )
                       d21 = t / d21
                       ! update elements in columns a(k) and a(k+1) as
                       ! dot products of rows of ( w(k) w(k+1) ) and columns
                       ! of d**(-1)
                       do j = k + 2, n
                          a( j, k ) = d21*( d11*w( j, k )-w( j, k+1 ) )
                          a( j, k+1 ) = d21*( d22*w( j, k+1 )-w( j, k ) )
                       end do
                    end if
                    ! copy d(k) to a
                    a( k, k ) = w( k, k )
                    a( k+1, k ) = w( k+1, k )
                    a( k+1, k+1 ) = w( k+1, k+1 )
                 end if
              end if
              ! store details of the interchanges in ipiv
              if( kstep==1_${ik}$ ) then
                 ipiv( k ) = kp
              else
                 ipiv( k ) = -kp
                 ipiv( k+1 ) = -kp
              end if
              ! increase k and return to the start of the main loop
              k = k + kstep
              go to 70
              90 continue
              ! update the lower triangle of a22 (= a(k:n,k:n)) as
              ! a22 := a22 - l21*d*l21**t = a22 - l21*w**t
              ! computing blocks of nb columns at a time
              do j = k, n, nb
                 jb = min( nb, n-j+1 )
                 ! update the lower triangle of the diagonal block
                 do jj = j, j + jb - 1
                    call stdlib${ii}$_sgemv( 'NO TRANSPOSE', j+jb-jj, k-1, -one,a( jj, 1_${ik}$ ), lda, w( jj, &
                              1_${ik}$ ), ldw, one,a( jj, jj ), 1_${ik}$ )
                 end do
                 ! update the rectangular subdiagonal block
                 if( j+jb<=n )call stdlib${ii}$_sgemm( 'NO TRANSPOSE', 'TRANSPOSE', n-j-jb+1, jb,k-1, -&
                           one, a( j+jb, 1_${ik}$ ), lda, w( j, 1_${ik}$ ), ldw,one, a( j+jb, j ), lda )
              end do
              ! put l21 in standard form by partially undoing the interchanges
              ! of rows in columns 1:k-1 looping backwards from k-1 to 1
              j = k - 1_${ik}$
              120 continue
                 ! undo the interchanges (if any) of rows jj and jp at each
                 ! step j
                 ! (here, j is a diagonal index)
                 jj = j
                 jp = ipiv( j )
                 if( jp<0_${ik}$ ) then
                    jp = -jp
                    ! (here, j is a diagonal index)
                    j = j - 1_${ik}$
                 end if
                 ! (note: here, j is used to determine row length. length j
                 ! of the rows to swap back doesn't include diagonal element)
                 j = j - 1_${ik}$
                 if( jp/=jj .and. j>=1_${ik}$ )call stdlib${ii}$_sswap( j, a( jp, 1_${ik}$ ), lda, a( jj, 1_${ik}$ ), lda )
                           
              if( j>1 )go to 120
              ! set kb to the number of columns factorized
              kb = k - 1_${ik}$
           end if
           return
     end subroutine stdlib${ii}$_slasyf

     pure module subroutine stdlib${ii}$_dlasyf( uplo, n, nb, kb, a, lda, ipiv, w, ldw, info )
     !! DLASYF computes a partial factorization of a real symmetric matrix A
     !! using the Bunch-Kaufman diagonal pivoting method. The partial
     !! factorization has the form:
     !! A  =  ( I  U12 ) ( A11  0  ) (  I       0    )  if UPLO = 'U', or:
     !! ( 0  U22 ) (  0   D  ) ( U12**T U22**T )
     !! A  =  ( L11  0 ) (  D   0  ) ( L11**T L21**T )  if UPLO = 'L'
     !! ( L21  I ) (  0  A22 ) (  0       I    )
     !! where the order of D is at most NB. The actual order is returned in
     !! the argument KB, and is either NB or NB-1, or N if N <= NB.
     !! DLASYF is an auxiliary routine called by DSYTRF. It uses blocked code
     !! (calling Level 3 BLAS) to update the submatrix A11 (if UPLO = 'U') or
     !! A22 (if UPLO = 'L').
        ! -- 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, kb
           integer(${ik}$), intent(in) :: lda, ldw, n, nb
           ! Array Arguments 
           integer(${ik}$), intent(out) :: ipiv(*)
           real(dp), intent(inout) :: a(lda,*)
           real(dp), intent(out) :: w(ldw,*)
        ! =====================================================================
           ! Parameters 
           real(dp), parameter :: sevten = 17.0e+0_dp
           
           
           ! Local Scalars 
           integer(${ik}$) :: imax, j, jb, jj, jmax, jp, k, kk, kkw, kp, kstep, kw
           real(dp) :: absakk, alpha, colmax, d11, d21, d22, r1, rowmax, t
           ! Intrinsic Functions 
           ! Executable Statements 
           info = 0_${ik}$
           ! initialize alpha for use in choosing pivot block size.
           alpha = ( one+sqrt( sevten ) ) / eight
           if( stdlib_lsame( uplo, 'U' ) ) then
              ! factorize the trailing columns of a using the upper triangle
              ! of a and working backwards, and compute the matrix w = u12*d
              ! for use in updating a11
              ! k is the main loop index, decreasing from n in steps of 1 or 2
              ! kw is the column of w which corresponds to column k of a
              k = n
              10 continue
              kw = nb + k - n
              ! exit from loop
              if( ( k<=n-nb+1 .and. nb<n ) .or. k<1 )go to 30
              ! copy column k of a to column kw of w and update it
              call stdlib${ii}$_dcopy( k, a( 1_${ik}$, k ), 1_${ik}$, w( 1_${ik}$, kw ), 1_${ik}$ )
              if( k<n )call stdlib${ii}$_dgemv( 'NO TRANSPOSE', k, n-k, -one, a( 1_${ik}$, k+1 ), lda,w( k, kw+&
                        1_${ik}$ ), ldw, one, w( 1_${ik}$, kw ), 1_${ik}$ )
              kstep = 1_${ik}$
              ! determine rows and columns to be interchanged and whether
              ! a 1-by-1 or 2-by-2 pivot block will be used
              absakk = abs( w( k, kw ) )
              ! imax is the row-index of the largest off-diagonal element in
              ! column k, and colmax is its absolute value.
              ! determine both colmax and imax.
              if( k>1_${ik}$ ) then
                 imax = stdlib${ii}$_idamax( k-1, w( 1_${ik}$, kw ), 1_${ik}$ )
                 colmax = abs( w( imax, kw ) )
              else
                 colmax = zero
              end if
              if( max( absakk, colmax )==zero ) then
                 ! column k is zero or underflow: set info and continue
                 if( info==0_${ik}$ )info = k
                 kp = k
              else
                 if( absakk>=alpha*colmax ) then
                    ! no interchange, use 1-by-1 pivot block
                    kp = k
                 else
                    ! copy column imax to column kw-1 of w and update it
                    call stdlib${ii}$_dcopy( imax, a( 1_${ik}$, imax ), 1_${ik}$, w( 1_${ik}$, kw-1 ), 1_${ik}$ )
                    call stdlib${ii}$_dcopy( k-imax, a( imax, imax+1 ), lda,w( imax+1, kw-1 ), 1_${ik}$ )
                              
                    if( k<n )call stdlib${ii}$_dgemv( 'NO TRANSPOSE', k, n-k, -one, a( 1_${ik}$, k+1 ),lda, w( &
                              imax, kw+1 ), ldw, one,w( 1_${ik}$, kw-1 ), 1_${ik}$ )
                    ! jmax is the column-index of the largest off-diagonal
                    ! element in row imax, and rowmax is its absolute value
                    jmax = imax + stdlib${ii}$_idamax( k-imax, w( imax+1, kw-1 ), 1_${ik}$ )
                    rowmax = abs( w( jmax, kw-1 ) )
                    if( imax>1_${ik}$ ) then
                       jmax = stdlib${ii}$_idamax( imax-1, w( 1_${ik}$, kw-1 ), 1_${ik}$ )
                       rowmax = max( rowmax, abs( w( jmax, kw-1 ) ) )
                    end if
                    if( absakk>=alpha*colmax*( colmax / rowmax ) ) then
                       ! no interchange, use 1-by-1 pivot block
                       kp = k
                    else if( abs( w( imax, kw-1 ) )>=alpha*rowmax ) then
                       ! interchange rows and columns k and imax, use 1-by-1
                       ! pivot block
                       kp = imax
                       ! copy column kw-1 of w to column kw of w
                       call stdlib${ii}$_dcopy( k, w( 1_${ik}$, kw-1 ), 1_${ik}$, w( 1_${ik}$, kw ), 1_${ik}$ )
                    else
                       ! interchange rows and columns k-1 and imax, use 2-by-2
                       ! pivot block
                       kp = imax
                       kstep = 2_${ik}$
                    end if
                 end if
                 ! ============================================================
                 ! kk is the column of a where pivoting step stopped
                 kk = k - kstep + 1_${ik}$
                 ! kkw is the column of w which corresponds to column kk of a
                 kkw = nb + kk - n
                 ! interchange rows and columns kp and kk.
                 ! updated column kp is already stored in column kkw of w.
                 if( kp/=kk ) then
                    ! copy non-updated column kk to column kp of submatrix a
                    ! at step k. no need to copy element into column k
                    ! (or k and k-1 for 2-by-2 pivot) of a, since these columns
                    ! will be later overwritten.
                    a( kp, kp ) = a( kk, kk )
                    call stdlib${ii}$_dcopy( kk-1-kp, a( kp+1, kk ), 1_${ik}$, a( kp, kp+1 ),lda )
                    if( kp>1_${ik}$ )call stdlib${ii}$_dcopy( kp-1, a( 1_${ik}$, kk ), 1_${ik}$, a( 1_${ik}$, kp ), 1_${ik}$ )
                    ! interchange rows kk and kp in last k+1 to n columns of a
                    ! (columns k (or k and k-1 for 2-by-2 pivot) of a will be
                    ! later overwritten). interchange rows kk and kp
                    ! in last kkw to nb columns of w.
                    if( k<n )call stdlib${ii}$_dswap( n-k, a( kk, k+1 ), lda, a( kp, k+1 ),lda )
                    call stdlib${ii}$_dswap( n-kk+1, w( kk, kkw ), ldw, w( kp, kkw ),ldw )
                 end if
                 if( kstep==1_${ik}$ ) then
                    ! 1-by-1 pivot block d(k): column kw of w now holds
                    ! w(kw) = u(k)*d(k),
                    ! where u(k) is the k-th column of u
                    ! store subdiag. elements of column u(k)
                    ! and 1-by-1 block d(k) in column k of a.
                    ! note: diagonal element u(k,k) is a unit element
                    ! and not stored.
                       ! a(k,k) := d(k,k) = w(k,kw)
                       ! a(1:k-1,k) := u(1:k-1,k) = w(1:k-1,kw)/d(k,k)
                    call stdlib${ii}$_dcopy( k, w( 1_${ik}$, kw ), 1_${ik}$, a( 1_${ik}$, k ), 1_${ik}$ )
                    r1 = one / a( k, k )
                    call stdlib${ii}$_dscal( k-1, r1, a( 1_${ik}$, k ), 1_${ik}$ )
                 else
                    ! 2-by-2 pivot block d(k): columns kw and kw-1 of w now hold
                    ! ( w(kw-1) w(kw) ) = ( u(k-1) u(k) )*d(k)
                    ! where u(k) and u(k-1) are the k-th and (k-1)-th columns
                    ! of u
                    ! store u(1:k-2,k-1) and u(1:k-2,k) and 2-by-2
                    ! block d(k-1:k,k-1:k) in columns k-1 and k of a.
                    ! note: 2-by-2 diagonal block u(k-1:k,k-1:k) is a unit
                    ! block and not stored.
                       ! a(k-1:k,k-1:k) := d(k-1:k,k-1:k) = w(k-1:k,kw-1:kw)
                       ! a(1:k-2,k-1:k) := u(1:k-2,k:k-1:k) =
                       ! = w(1:k-2,kw-1:kw) * ( d(k-1:k,k-1:k)**(-1) )
                    if( k>2_${ik}$ ) then
                       ! compose the columns of the inverse of 2-by-2 pivot
                       ! block d in the following way to reduce the number
                       ! of flops when we myltiply panel ( w(kw-1) w(kw) ) by
                       ! this inverse
                       ! d**(-1) = ( d11 d21 )**(-1) =
                                 ! ( d21 d22 )
                       ! = 1/(d11*d22-d21**2) * ( ( d22 ) (-d21 ) ) =
                                              ! ( (-d21 ) ( d11 ) )
                       ! = 1/d21 * 1/((d11/d21)*(d22/d21)-1) *
                         ! * ( ( d22/d21 ) (      -1 ) ) =
                           ! ( (      -1 ) ( d11/d21 ) )
                       ! = 1/d21 * 1/(d22*d11-1) * ( ( d11 ) (  -1 ) ) =
                                                 ! ( ( -1  ) ( d22 ) )
                       ! = 1/d21 * t * ( ( d11 ) (  -1 ) )
                                     ! ( (  -1 ) ( d22 ) )
                       ! = d21 * ( ( d11 ) (  -1 ) )
                               ! ( (  -1 ) ( d22 ) )
                       d21 = w( k-1, kw )
                       d11 = w( k, kw ) / d21
                       d22 = w( k-1, kw-1 ) / d21
                       t = one / ( d11*d22-one )
                       d21 = t / d21
                       ! update elements in columns a(k-1) and a(k) as
                       ! dot products of rows of ( w(kw-1) w(kw) ) and columns
                       ! of d**(-1)
                       do j = 1, k - 2
                          a( j, k-1 ) = d21*( d11*w( j, kw-1 )-w( j, kw ) )
                          a( j, k ) = d21*( d22*w( j, kw )-w( j, kw-1 ) )
                       end do
                    end if
                    ! copy d(k) to a
                    a( k-1, k-1 ) = w( k-1, kw-1 )
                    a( k-1, k ) = w( k-1, kw )
                    a( k, k ) = w( k, kw )
                 end if
              end if
              ! store details of the interchanges in ipiv
              if( kstep==1_${ik}$ ) then
                 ipiv( k ) = kp
              else
                 ipiv( k ) = -kp
                 ipiv( k-1 ) = -kp
              end if
              ! decrease k and return to the start of the main loop
              k = k - kstep
              go to 10
              30 continue
              ! update the upper triangle of a11 (= a(1:k,1:k)) as
              ! a11 := a11 - u12*d*u12**t = a11 - u12*w**t
              ! computing blocks of nb columns at a time
              do j = ( ( k-1 ) / nb )*nb + 1, 1, -nb
                 jb = min( nb, k-j+1 )
                 ! update the upper triangle of the diagonal block
                 do jj = j, j + jb - 1
                    call stdlib${ii}$_dgemv( 'NO TRANSPOSE', jj-j+1, n-k, -one,a( j, k+1 ), lda, w( jj, &
                              kw+1 ), ldw, one,a( j, jj ), 1_${ik}$ )
                 end do
                 ! update the rectangular superdiagonal block
                 call stdlib${ii}$_dgemm( 'NO TRANSPOSE', 'TRANSPOSE', j-1, jb, n-k, -one,a( 1_${ik}$, k+1 ), &
                           lda, w( j, kw+1 ), ldw, one,a( 1_${ik}$, j ), lda )
              end do
              ! put u12 in standard form by partially undoing the interchanges
              ! in columns k+1:n looping backwards from k+1 to n
              j = k + 1_${ik}$
              60 continue
                 ! undo the interchanges (if any) of rows jj and jp at each
                 ! step j
                 ! (here, j is a diagonal index)
                 jj = j
                 jp = ipiv( j )
                 if( jp<0_${ik}$ ) then
                    jp = -jp
                    ! (here, j is a diagonal index)
                    j = j + 1_${ik}$
                 end if
                 ! (note: here, j is used to determine row length. length n-j+1
                 ! of the rows to swap back doesn't include diagonal element)
                 j = j + 1_${ik}$
                 if( jp/=jj .and. j<=n )call stdlib${ii}$_dswap( n-j+1, a( jp, j ), lda, a( jj, j ), &
                           lda )
              if( j<n )go to 60
              ! set kb to the number of columns factorized
              kb = n - k
           else
              ! factorize the leading columns of a using the lower triangle
              ! of a and working forwards, and compute the matrix w = l21*d
              ! for use in updating a22
              ! k is the main loop index, increasing from 1 in steps of 1 or 2
              k = 1_${ik}$
              70 continue
              ! exit from loop
              if( ( k>=nb .and. nb<n ) .or. k>n )go to 90
              ! copy column k of a to column k of w and update it
              call stdlib${ii}$_dcopy( n-k+1, a( k, k ), 1_${ik}$, w( k, k ), 1_${ik}$ )
              call stdlib${ii}$_dgemv( 'NO TRANSPOSE', n-k+1, k-1, -one, a( k, 1_${ik}$ ), lda,w( k, 1_${ik}$ ), ldw, &
                        one, w( k, k ), 1_${ik}$ )
              kstep = 1_${ik}$
              ! determine rows and columns to be interchanged and whether
              ! a 1-by-1 or 2-by-2 pivot block will be used
              absakk = abs( w( k, k ) )
              ! imax is the row-index of the largest off-diagonal element in
              ! column k, and colmax is its absolute value.
              ! determine both colmax and imax.
              if( k<n ) then
                 imax = k + stdlib${ii}$_idamax( n-k, w( k+1, k ), 1_${ik}$ )
                 colmax = abs( w( imax, k ) )
              else
                 colmax = zero
              end if
              if( max( absakk, colmax )==zero ) then
                 ! column k is zero or underflow: set info and continue
                 if( info==0_${ik}$ )info = k
                 kp = k
              else
                 if( absakk>=alpha*colmax ) then
                    ! no interchange, use 1-by-1 pivot block
                    kp = k
                 else
                    ! copy column imax to column k+1 of w and update it
                    call stdlib${ii}$_dcopy( imax-k, a( imax, k ), lda, w( k, k+1 ), 1_${ik}$ )
                    call stdlib${ii}$_dcopy( n-imax+1, a( imax, imax ), 1_${ik}$, w( imax, k+1 ),1_${ik}$ )
                    call stdlib${ii}$_dgemv( 'NO TRANSPOSE', n-k+1, k-1, -one, a( k, 1_${ik}$ ),lda, w( imax, &
                              1_${ik}$ ), ldw, one, w( k, k+1 ), 1_${ik}$ )
                    ! jmax is the column-index of the largest off-diagonal
                    ! element in row imax, and rowmax is its absolute value
                    jmax = k - 1_${ik}$ + stdlib${ii}$_idamax( imax-k, w( k, k+1 ), 1_${ik}$ )
                    rowmax = abs( w( jmax, k+1 ) )
                    if( imax<n ) then
                       jmax = imax + stdlib${ii}$_idamax( n-imax, w( imax+1, k+1 ), 1_${ik}$ )
                       rowmax = max( rowmax, abs( w( jmax, k+1 ) ) )
                    end if
                    if( absakk>=alpha*colmax*( colmax / rowmax ) ) then
                       ! no interchange, use 1-by-1 pivot block
                       kp = k
                    else if( abs( w( imax, k+1 ) )>=alpha*rowmax ) then
                       ! interchange rows and columns k and imax, use 1-by-1
                       ! pivot block
                       kp = imax
                       ! copy column k+1 of w to column k of w
                       call stdlib${ii}$_dcopy( n-k+1, w( k, k+1 ), 1_${ik}$, w( k, k ), 1_${ik}$ )
                    else
                       ! interchange rows and columns k+1 and imax, use 2-by-2
                       ! pivot block
                       kp = imax
                       kstep = 2_${ik}$
                    end if
                 end if
                 ! ============================================================
                 ! kk is the column of a where pivoting step stopped
                 kk = k + kstep - 1_${ik}$
                 ! interchange rows and columns kp and kk.
                 ! updated column kp is already stored in column kk of w.
                 if( kp/=kk ) then
                    ! copy non-updated column kk to column kp of submatrix a
                    ! at step k. no need to copy element into column k
                    ! (or k and k+1 for 2-by-2 pivot) of a, since these columns
                    ! will be later overwritten.
                    a( kp, kp ) = a( kk, kk )
                    call stdlib${ii}$_dcopy( kp-kk-1, a( kk+1, kk ), 1_${ik}$, a( kp, kk+1 ),lda )
                    if( kp<n )call stdlib${ii}$_dcopy( n-kp, a( kp+1, kk ), 1_${ik}$, a( kp+1, kp ), 1_${ik}$ )
                              
                    ! interchange rows kk and kp in first k-1 columns of a
                    ! (columns k (or k and k+1 for 2-by-2 pivot) of a will be
                    ! later overwritten). interchange rows kk and kp
                    ! in first kk columns of w.
                    if( k>1_${ik}$ )call stdlib${ii}$_dswap( k-1, a( kk, 1_${ik}$ ), lda, a( kp, 1_${ik}$ ), lda )
                    call stdlib${ii}$_dswap( kk, w( kk, 1_${ik}$ ), ldw, w( kp, 1_${ik}$ ), ldw )
                 end if
                 if( kstep==1_${ik}$ ) then
                    ! 1-by-1 pivot block d(k): column k of w now holds
                    ! w(k) = l(k)*d(k),
                    ! where l(k) is the k-th column of l
                    ! store subdiag. elements of column l(k)
                    ! and 1-by-1 block d(k) in column k of a.
                    ! (note: diagonal element l(k,k) is a unit element
                    ! and not stored)
                       ! a(k,k) := d(k,k) = w(k,k)
                       ! a(k+1:n,k) := l(k+1:n,k) = w(k+1:n,k)/d(k,k)
                    call stdlib${ii}$_dcopy( n-k+1, w( k, k ), 1_${ik}$, a( k, k ), 1_${ik}$ )
                    if( k<n ) then
                       r1 = one / a( k, k )
                       call stdlib${ii}$_dscal( n-k, r1, a( k+1, k ), 1_${ik}$ )
                    end if
                 else
                    ! 2-by-2 pivot block d(k): columns k and k+1 of w now hold
                    ! ( w(k) w(k+1) ) = ( l(k) l(k+1) )*d(k)
                    ! where l(k) and l(k+1) are the k-th and (k+1)-th columns
                    ! of l
                    ! store l(k+2:n,k) and l(k+2:n,k+1) and 2-by-2
                    ! block d(k:k+1,k:k+1) in columns k and k+1 of a.
                    ! (note: 2-by-2 diagonal block l(k:k+1,k:k+1) is a unit
                    ! block and not stored)
                       ! a(k:k+1,k:k+1) := d(k:k+1,k:k+1) = w(k:k+1,k:k+1)
                       ! a(k+2:n,k:k+1) := l(k+2:n,k:k+1) =
                       ! = w(k+2:n,k:k+1) * ( d(k:k+1,k:k+1)**(-1) )
                    if( k<n-1 ) then
                       ! compose the columns of the inverse of 2-by-2 pivot
                       ! block d in the following way to reduce the number
                       ! of flops when we myltiply panel ( w(k) w(k+1) ) by
                       ! this inverse
                       ! d**(-1) = ( d11 d21 )**(-1) =
                                 ! ( d21 d22 )
                       ! = 1/(d11*d22-d21**2) * ( ( d22 ) (-d21 ) ) =
                                              ! ( (-d21 ) ( d11 ) )
                       ! = 1/d21 * 1/((d11/d21)*(d22/d21)-1) *
                         ! * ( ( d22/d21 ) (      -1 ) ) =
                           ! ( (      -1 ) ( d11/d21 ) )
                       ! = 1/d21 * 1/(d22*d11-1) * ( ( d11 ) (  -1 ) ) =
                                                 ! ( ( -1  ) ( d22 ) )
                       ! = 1/d21 * t * ( ( d11 ) (  -1 ) )
                                     ! ( (  -1 ) ( d22 ) )
                       ! = d21 * ( ( d11 ) (  -1 ) )
                               ! ( (  -1 ) ( d22 ) )
                       d21 = w( k+1, k )
                       d11 = w( k+1, k+1 ) / d21
                       d22 = w( k, k ) / d21
                       t = one / ( d11*d22-one )
                       d21 = t / d21
                       ! update elements in columns a(k) and a(k+1) as
                       ! dot products of rows of ( w(k) w(k+1) ) and columns
                       ! of d**(-1)
                       do j = k + 2, n
                          a( j, k ) = d21*( d11*w( j, k )-w( j, k+1 ) )
                          a( j, k+1 ) = d21*( d22*w( j, k+1 )-w( j, k ) )
                       end do
                    end if
                    ! copy d(k) to a
                    a( k, k ) = w( k, k )
                    a( k+1, k ) = w( k+1, k )
                    a( k+1, k+1 ) = w( k+1, k+1 )
                 end if
              end if
              ! store details of the interchanges in ipiv
              if( kstep==1_${ik}$ ) then
                 ipiv( k ) = kp
              else
                 ipiv( k ) = -kp
                 ipiv( k+1 ) = -kp
              end if
              ! increase k and return to the start of the main loop
              k = k + kstep
              go to 70
              90 continue
              ! update the lower triangle of a22 (= a(k:n,k:n)) as
              ! a22 := a22 - l21*d*l21**t = a22 - l21*w**t
              ! computing blocks of nb columns at a time
              do j = k, n, nb
                 jb = min( nb, n-j+1 )
                 ! update the lower triangle of the diagonal block
                 do jj = j, j + jb - 1
                    call stdlib${ii}$_dgemv( 'NO TRANSPOSE', j+jb-jj, k-1, -one,a( jj, 1_${ik}$ ), lda, w( jj, &
                              1_${ik}$ ), ldw, one,a( jj, jj ), 1_${ik}$ )
                 end do
                 ! update the rectangular subdiagonal block
                 if( j+jb<=n )call stdlib${ii}$_dgemm( 'NO TRANSPOSE', 'TRANSPOSE', n-j-jb+1, jb,k-1, -&
                           one, a( j+jb, 1_${ik}$ ), lda, w( j, 1_${ik}$ ), ldw,one, a( j+jb, j ), lda )
              end do
              ! put l21 in standard form by partially undoing the interchanges
              ! of rows in columns 1:k-1 looping backwards from k-1 to 1
              j = k - 1_${ik}$
              120 continue
                 ! undo the interchanges (if any) of rows jj and jp at each
                 ! step j
                 ! (here, j is a diagonal index)
                 jj = j
                 jp = ipiv( j )
                 if( jp<0_${ik}$ ) then
                    jp = -jp
                    ! (here, j is a diagonal index)
                    j = j - 1_${ik}$
                 end if
                 ! (note: here, j is used to determine row length. length j
                 ! of the rows to swap back doesn't include diagonal element)
                 j = j - 1_${ik}$
                 if( jp/=jj .and. j>=1_${ik}$ )call stdlib${ii}$_dswap( j, a( jp, 1_${ik}$ ), lda, a( jj, 1_${ik}$ ), lda )
                           
              if( j>1 )go to 120
              ! set kb to the number of columns factorized
              kb = k - 1_${ik}$
           end if
           return
     end subroutine stdlib${ii}$_dlasyf

#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
     pure module subroutine stdlib${ii}$_${ri}$lasyf( uplo, n, nb, kb, a, lda, ipiv, w, ldw, info )
     !! DLASYF: computes a partial factorization of a real symmetric matrix A
     !! using the Bunch-Kaufman diagonal pivoting method. The partial
     !! factorization has the form:
     !! A  =  ( I  U12 ) ( A11  0  ) (  I       0    )  if UPLO = 'U', or:
     !! ( 0  U22 ) (  0   D  ) ( U12**T U22**T )
     !! A  =  ( L11  0 ) (  D   0  ) ( L11**T L21**T )  if UPLO = 'L'
     !! ( L21  I ) (  0  A22 ) (  0       I    )
     !! where the order of D is at most NB. The actual order is returned in
     !! the argument KB, and is either NB or NB-1, or N if N <= NB.
     !! DLASYF is an auxiliary routine called by DSYTRF. It uses blocked code
     !! (calling Level 3 BLAS) to update the submatrix A11 (if UPLO = 'U') or
     !! A22 (if UPLO = 'L').
        ! -- 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, kb
           integer(${ik}$), intent(in) :: lda, ldw, n, nb
           ! Array Arguments 
           integer(${ik}$), intent(out) :: ipiv(*)
           real(${rk}$), intent(inout) :: a(lda,*)
           real(${rk}$), intent(out) :: w(ldw,*)
        ! =====================================================================
           ! Parameters 
           real(${rk}$), parameter :: sevten = 17.0e+0_${rk}$
           
           
           ! Local Scalars 
           integer(${ik}$) :: imax, j, jb, jj, jmax, jp, k, kk, kkw, kp, kstep, kw
           real(${rk}$) :: absakk, alpha, colmax, d11, d21, d22, r1, rowmax, t
           ! Intrinsic Functions 
           ! Executable Statements 
           info = 0_${ik}$
           ! initialize alpha for use in choosing pivot block size.
           alpha = ( one+sqrt( sevten ) ) / eight
           if( stdlib_lsame( uplo, 'U' ) ) then
              ! factorize the trailing columns of a using the upper triangle
              ! of a and working backwards, and compute the matrix w = u12*d
              ! for use in updating a11
              ! k is the main loop index, decreasing from n in steps of 1 or 2
              ! kw is the column of w which corresponds to column k of a
              k = n
              10 continue
              kw = nb + k - n
              ! exit from loop
              if( ( k<=n-nb+1 .and. nb<n ) .or. k<1 )go to 30
              ! copy column k of a to column kw of w and update it
              call stdlib${ii}$_${ri}$copy( k, a( 1_${ik}$, k ), 1_${ik}$, w( 1_${ik}$, kw ), 1_${ik}$ )
              if( k<n )call stdlib${ii}$_${ri}$gemv( 'NO TRANSPOSE', k, n-k, -one, a( 1_${ik}$, k+1 ), lda,w( k, kw+&
                        1_${ik}$ ), ldw, one, w( 1_${ik}$, kw ), 1_${ik}$ )
              kstep = 1_${ik}$
              ! determine rows and columns to be interchanged and whether
              ! a 1-by-1 or 2-by-2 pivot block will be used
              absakk = abs( w( k, kw ) )
              ! imax is the row-index of the largest off-diagonal element in
              ! column k, and colmax is its absolute value.
              ! determine both colmax and imax.
              if( k>1_${ik}$ ) then
                 imax = stdlib${ii}$_i${ri}$amax( k-1, w( 1_${ik}$, kw ), 1_${ik}$ )
                 colmax = abs( w( imax, kw ) )
              else
                 colmax = zero
              end if
              if( max( absakk, colmax )==zero ) then
                 ! column k is zero or underflow: set info and continue
                 if( info==0_${ik}$ )info = k
                 kp = k
              else
                 if( absakk>=alpha*colmax ) then
                    ! no interchange, use 1-by-1 pivot block
                    kp = k
                 else
                    ! copy column imax to column kw-1 of w and update it
                    call stdlib${ii}$_${ri}$copy( imax, a( 1_${ik}$, imax ), 1_${ik}$, w( 1_${ik}$, kw-1 ), 1_${ik}$ )
                    call stdlib${ii}$_${ri}$copy( k-imax, a( imax, imax+1 ), lda,w( imax+1, kw-1 ), 1_${ik}$ )
                              
                    if( k<n )call stdlib${ii}$_${ri}$gemv( 'NO TRANSPOSE', k, n-k, -one, a( 1_${ik}$, k+1 ),lda, w( &
                              imax, kw+1 ), ldw, one,w( 1_${ik}$, kw-1 ), 1_${ik}$ )
                    ! jmax is the column-index of the largest off-diagonal
                    ! element in row imax, and rowmax is its absolute value
                    jmax = imax + stdlib${ii}$_i${ri}$amax( k-imax, w( imax+1, kw-1 ), 1_${ik}$ )
                    rowmax = abs( w( jmax, kw-1 ) )
                    if( imax>1_${ik}$ ) then
                       jmax = stdlib${ii}$_i${ri}$amax( imax-1, w( 1_${ik}$, kw-1 ), 1_${ik}$ )
                       rowmax = max( rowmax, abs( w( jmax, kw-1 ) ) )
                    end if
                    if( absakk>=alpha*colmax*( colmax / rowmax ) ) then
                       ! no interchange, use 1-by-1 pivot block
                       kp = k
                    else if( abs( w( imax, kw-1 ) )>=alpha*rowmax ) then
                       ! interchange rows and columns k and imax, use 1-by-1
                       ! pivot block
                       kp = imax
                       ! copy column kw-1 of w to column kw of w
                       call stdlib${ii}$_${ri}$copy( k, w( 1_${ik}$, kw-1 ), 1_${ik}$, w( 1_${ik}$, kw ), 1_${ik}$ )
                    else
                       ! interchange rows and columns k-1 and imax, use 2-by-2
                       ! pivot block
                       kp = imax
                       kstep = 2_${ik}$
                    end if
                 end if
                 ! ============================================================
                 ! kk is the column of a where pivoting step stopped
                 kk = k - kstep + 1_${ik}$
                 ! kkw is the column of w which corresponds to column kk of a
                 kkw = nb + kk - n
                 ! interchange rows and columns kp and kk.
                 ! updated column kp is already stored in column kkw of w.
                 if( kp/=kk ) then
                    ! copy non-updated column kk to column kp of submatrix a
                    ! at step k. no need to copy element into column k
                    ! (or k and k-1 for 2-by-2 pivot) of a, since these columns
                    ! will be later overwritten.
                    a( kp, kp ) = a( kk, kk )
                    call stdlib${ii}$_${ri}$copy( kk-1-kp, a( kp+1, kk ), 1_${ik}$, a( kp, kp+1 ),lda )
                    if( kp>1_${ik}$ )call stdlib${ii}$_${ri}$copy( kp-1, a( 1_${ik}$, kk ), 1_${ik}$, a( 1_${ik}$, kp ), 1_${ik}$ )
                    ! interchange rows kk and kp in last k+1 to n columns of a
                    ! (columns k (or k and k-1 for 2-by-2 pivot) of a will be
                    ! later overwritten). interchange rows kk and kp
                    ! in last kkw to nb columns of w.
                    if( k<n )call stdlib${ii}$_${ri}$swap( n-k, a( kk, k+1 ), lda, a( kp, k+1 ),lda )
                    call stdlib${ii}$_${ri}$swap( n-kk+1, w( kk, kkw ), ldw, w( kp, kkw ),ldw )
                 end if
                 if( kstep==1_${ik}$ ) then
                    ! 1-by-1 pivot block d(k): column kw of w now holds
                    ! w(kw) = u(k)*d(k),
                    ! where u(k) is the k-th column of u
                    ! store subdiag. elements of column u(k)
                    ! and 1-by-1 block d(k) in column k of a.
                    ! note: diagonal element u(k,k) is a unit element
                    ! and not stored.
                       ! a(k,k) := d(k,k) = w(k,kw)
                       ! a(1:k-1,k) := u(1:k-1,k) = w(1:k-1,kw)/d(k,k)
                    call stdlib${ii}$_${ri}$copy( k, w( 1_${ik}$, kw ), 1_${ik}$, a( 1_${ik}$, k ), 1_${ik}$ )
                    r1 = one / a( k, k )
                    call stdlib${ii}$_${ri}$scal( k-1, r1, a( 1_${ik}$, k ), 1_${ik}$ )
                 else
                    ! 2-by-2 pivot block d(k): columns kw and kw-1 of w now hold
                    ! ( w(kw-1) w(kw) ) = ( u(k-1) u(k) )*d(k)
                    ! where u(k) and u(k-1) are the k-th and (k-1)-th columns
                    ! of u
                    ! store u(1:k-2,k-1) and u(1:k-2,k) and 2-by-2
                    ! block d(k-1:k,k-1:k) in columns k-1 and k of a.
                    ! note: 2-by-2 diagonal block u(k-1:k,k-1:k) is a unit
                    ! block and not stored.
                       ! a(k-1:k,k-1:k) := d(k-1:k,k-1:k) = w(k-1:k,kw-1:kw)
                       ! a(1:k-2,k-1:k) := u(1:k-2,k:k-1:k) =
                       ! = w(1:k-2,kw-1:kw) * ( d(k-1:k,k-1:k)**(-1) )
                    if( k>2_${ik}$ ) then
                       ! compose the columns of the inverse of 2-by-2 pivot
                       ! block d in the following way to reduce the number
                       ! of flops when we myltiply panel ( w(kw-1) w(kw) ) by
                       ! this inverse
                       ! d**(-1) = ( d11 d21 )**(-1) =
                                 ! ( d21 d22 )
                       ! = 1/(d11*d22-d21**2) * ( ( d22 ) (-d21 ) ) =
                                              ! ( (-d21 ) ( d11 ) )
                       ! = 1/d21 * 1/((d11/d21)*(d22/d21)-1) *
                         ! * ( ( d22/d21 ) (      -1 ) ) =
                           ! ( (      -1 ) ( d11/d21 ) )
                       ! = 1/d21 * 1/(d22*d11-1) * ( ( d11 ) (  -1 ) ) =
                                                 ! ( ( -1  ) ( d22 ) )
                       ! = 1/d21 * t * ( ( d11 ) (  -1 ) )
                                     ! ( (  -1 ) ( d22 ) )
                       ! = d21 * ( ( d11 ) (  -1 ) )
                               ! ( (  -1 ) ( d22 ) )
                       d21 = w( k-1, kw )
                       d11 = w( k, kw ) / d21
                       d22 = w( k-1, kw-1 ) / d21
                       t = one / ( d11*d22-one )
                       d21 = t / d21
                       ! update elements in columns a(k-1) and a(k) as
                       ! dot products of rows of ( w(kw-1) w(kw) ) and columns
                       ! of d**(-1)
                       do j = 1, k - 2
                          a( j, k-1 ) = d21*( d11*w( j, kw-1 )-w( j, kw ) )
                          a( j, k ) = d21*( d22*w( j, kw )-w( j, kw-1 ) )
                       end do
                    end if
                    ! copy d(k) to a
                    a( k-1, k-1 ) = w( k-1, kw-1 )
                    a( k-1, k ) = w( k-1, kw )
                    a( k, k ) = w( k, kw )
                 end if
              end if
              ! store details of the interchanges in ipiv
              if( kstep==1_${ik}$ ) then
                 ipiv( k ) = kp
              else
                 ipiv( k ) = -kp
                 ipiv( k-1 ) = -kp
              end if
              ! decrease k and return to the start of the main loop
              k = k - kstep
              go to 10
              30 continue
              ! update the upper triangle of a11 (= a(1:k,1:k)) as
              ! a11 := a11 - u12*d*u12**t = a11 - u12*w**t
              ! computing blocks of nb columns at a time
              do j = ( ( k-1 ) / nb )*nb + 1, 1, -nb
                 jb = min( nb, k-j+1 )
                 ! update the upper triangle of the diagonal block
                 do jj = j, j + jb - 1
                    call stdlib${ii}$_${ri}$gemv( 'NO TRANSPOSE', jj-j+1, n-k, -one,a( j, k+1 ), lda, w( jj, &
                              kw+1 ), ldw, one,a( j, jj ), 1_${ik}$ )
                 end do
                 ! update the rectangular superdiagonal block
                 call stdlib${ii}$_${ri}$gemm( 'NO TRANSPOSE', 'TRANSPOSE', j-1, jb, n-k, -one,a( 1_${ik}$, k+1 ), &
                           lda, w( j, kw+1 ), ldw, one,a( 1_${ik}$, j ), lda )
              end do
              ! put u12 in standard form by partially undoing the interchanges
              ! in columns k+1:n looping backwards from k+1 to n
              j = k + 1_${ik}$
              60 continue
                 ! undo the interchanges (if any) of rows jj and jp at each
                 ! step j
                 ! (here, j is a diagonal index)
                 jj = j
                 jp = ipiv( j )
                 if( jp<0_${ik}$ ) then
                    jp = -jp
                    ! (here, j is a diagonal index)
                    j = j + 1_${ik}$
                 end if
                 ! (note: here, j is used to determine row length. length n-j+1
                 ! of the rows to swap back doesn't include diagonal element)
                 j = j + 1_${ik}$
                 if( jp/=jj .and. j<=n )call stdlib${ii}$_${ri}$swap( n-j+1, a( jp, j ), lda, a( jj, j ), &
                           lda )
              if( j<n )go to 60
              ! set kb to the number of columns factorized
              kb = n - k
           else
              ! factorize the leading columns of a using the lower triangle
              ! of a and working forwards, and compute the matrix w = l21*d
              ! for use in updating a22
              ! k is the main loop index, increasing from 1 in steps of 1 or 2
              k = 1_${ik}$
              70 continue
              ! exit from loop
              if( ( k>=nb .and. nb<n ) .or. k>n )go to 90
              ! copy column k of a to column k of w and update it
              call stdlib${ii}$_${ri}$copy( n-k+1, a( k, k ), 1_${ik}$, w( k, k ), 1_${ik}$ )
              call stdlib${ii}$_${ri}$gemv( 'NO TRANSPOSE', n-k+1, k-1, -one, a( k, 1_${ik}$ ), lda,w( k, 1_${ik}$ ), ldw, &
                        one, w( k, k ), 1_${ik}$ )
              kstep = 1_${ik}$
              ! determine rows and columns to be interchanged and whether
              ! a 1-by-1 or 2-by-2 pivot block will be used
              absakk = abs( w( k, k ) )
              ! imax is the row-index of the largest off-diagonal element in
              ! column k, and colmax is its absolute value.
              ! determine both colmax and imax.
              if( k<n ) then
                 imax = k + stdlib${ii}$_i${ri}$amax( n-k, w( k+1, k ), 1_${ik}$ )
                 colmax = abs( w( imax, k ) )
              else
                 colmax = zero
              end if
              if( max( absakk, colmax )==zero ) then
                 ! column k is zero or underflow: set info and continue
                 if( info==0_${ik}$ )info = k
                 kp = k
              else
                 if( absakk>=alpha*colmax ) then
                    ! no interchange, use 1-by-1 pivot block
                    kp = k
                 else
                    ! copy column imax to column k+1 of w and update it
                    call stdlib${ii}$_${ri}$copy( imax-k, a( imax, k ), lda, w( k, k+1 ), 1_${ik}$ )
                    call stdlib${ii}$_${ri}$copy( n-imax+1, a( imax, imax ), 1_${ik}$, w( imax, k+1 ),1_${ik}$ )
                    call stdlib${ii}$_${ri}$gemv( 'NO TRANSPOSE', n-k+1, k-1, -one, a( k, 1_${ik}$ ),lda, w( imax, &
                              1_${ik}$ ), ldw, one, w( k, k+1 ), 1_${ik}$ )
                    ! jmax is the column-index of the largest off-diagonal
                    ! element in row imax, and rowmax is its absolute value
                    jmax = k - 1_${ik}$ + stdlib${ii}$_i${ri}$amax( imax-k, w( k, k+1 ), 1_${ik}$ )
                    rowmax = abs( w( jmax, k+1 ) )
                    if( imax<n ) then
                       jmax = imax + stdlib${ii}$_i${ri}$amax( n-imax, w( imax+1, k+1 ), 1_${ik}$ )
                       rowmax = max( rowmax, abs( w( jmax, k+1 ) ) )
                    end if
                    if( absakk>=alpha*colmax*( colmax / rowmax ) ) then
                       ! no interchange, use 1-by-1 pivot block
                       kp = k
                    else if( abs( w( imax, k+1 ) )>=alpha*rowmax ) then
                       ! interchange rows and columns k and imax, use 1-by-1
                       ! pivot block
                       kp = imax
                       ! copy column k+1 of w to column k of w
                       call stdlib${ii}$_${ri}$copy( n-k+1, w( k, k+1 ), 1_${ik}$, w( k, k ), 1_${ik}$ )
                    else
                       ! interchange rows and columns k+1 and imax, use 2-by-2
                       ! pivot block
                       kp = imax
                       kstep = 2_${ik}$
                    end if
                 end if
                 ! ============================================================
                 ! kk is the column of a where pivoting step stopped
                 kk = k + kstep - 1_${ik}$
                 ! interchange rows and columns kp and kk.
                 ! updated column kp is already stored in column kk of w.
                 if( kp/=kk ) then
                    ! copy non-updated column kk to column kp of submatrix a
                    ! at step k. no need to copy element into column k
                    ! (or k and k+1 for 2-by-2 pivot) of a, since these columns
                    ! will be later overwritten.
                    a( kp, kp ) = a( kk, kk )
                    call stdlib${ii}$_${ri}$copy( kp-kk-1, a( kk+1, kk ), 1_${ik}$, a( kp, kk+1 ),lda )
                    if( kp<n )call stdlib${ii}$_${ri}$copy( n-kp, a( kp+1, kk ), 1_${ik}$, a( kp+1, kp ), 1_${ik}$ )
                              
                    ! interchange rows kk and kp in first k-1 columns of a
                    ! (columns k (or k and k+1 for 2-by-2 pivot) of a will be
                    ! later overwritten). interchange rows kk and kp
                    ! in first kk columns of w.
                    if( k>1_${ik}$ )call stdlib${ii}$_${ri}$swap( k-1, a( kk, 1_${ik}$ ), lda, a( kp, 1_${ik}$ ), lda )
                    call stdlib${ii}$_${ri}$swap( kk, w( kk, 1_${ik}$ ), ldw, w( kp, 1_${ik}$ ), ldw )
                 end if
                 if( kstep==1_${ik}$ ) then
                    ! 1-by-1 pivot block d(k): column k of w now holds
                    ! w(k) = l(k)*d(k),
                    ! where l(k) is the k-th column of l
                    ! store subdiag. elements of column l(k)
                    ! and 1-by-1 block d(k) in column k of a.
                    ! (note: diagonal element l(k,k) is a unit element
                    ! and not stored)
                       ! a(k,k) := d(k,k) = w(k,k)
                       ! a(k+1:n,k) := l(k+1:n,k) = w(k+1:n,k)/d(k,k)
                    call stdlib${ii}$_${ri}$copy( n-k+1, w( k, k ), 1_${ik}$, a( k, k ), 1_${ik}$ )
                    if( k<n ) then
                       r1 = one / a( k, k )
                       call stdlib${ii}$_${ri}$scal( n-k, r1, a( k+1, k ), 1_${ik}$ )
                    end if
                 else
                    ! 2-by-2 pivot block d(k): columns k and k+1 of w now hold
                    ! ( w(k) w(k+1) ) = ( l(k) l(k+1) )*d(k)
                    ! where l(k) and l(k+1) are the k-th and (k+1)-th columns
                    ! of l
                    ! store l(k+2:n,k) and l(k+2:n,k+1) and 2-by-2
                    ! block d(k:k+1,k:k+1) in columns k and k+1 of a.
                    ! (note: 2-by-2 diagonal block l(k:k+1,k:k+1) is a unit
                    ! block and not stored)
                       ! a(k:k+1,k:k+1) := d(k:k+1,k:k+1) = w(k:k+1,k:k+1)
                       ! a(k+2:n,k:k+1) := l(k+2:n,k:k+1) =
                       ! = w(k+2:n,k:k+1) * ( d(k:k+1,k:k+1)**(-1) )
                    if( k<n-1 ) then
                       ! compose the columns of the inverse of 2-by-2 pivot
                       ! block d in the following way to reduce the number
                       ! of flops when we myltiply panel ( w(k) w(k+1) ) by
                       ! this inverse
                       ! d**(-1) = ( d11 d21 )**(-1) =
                                 ! ( d21 d22 )
                       ! = 1/(d11*d22-d21**2) * ( ( d22 ) (-d21 ) ) =
                                              ! ( (-d21 ) ( d11 ) )
                       ! = 1/d21 * 1/((d11/d21)*(d22/d21)-1) *
                         ! * ( ( d22/d21 ) (      -1 ) ) =
                           ! ( (      -1 ) ( d11/d21 ) )
                       ! = 1/d21 * 1/(d22*d11-1) * ( ( d11 ) (  -1 ) ) =
                                                 ! ( ( -1  ) ( d22 ) )
                       ! = 1/d21 * t * ( ( d11 ) (  -1 ) )
                                     ! ( (  -1 ) ( d22 ) )
                       ! = d21 * ( ( d11 ) (  -1 ) )
                               ! ( (  -1 ) ( d22 ) )
                       d21 = w( k+1, k )
                       d11 = w( k+1, k+1 ) / d21
                       d22 = w( k, k ) / d21
                       t = one / ( d11*d22-one )
                       d21 = t / d21
                       ! update elements in columns a(k) and a(k+1) as
                       ! dot products of rows of ( w(k) w(k+1) ) and columns
                       ! of d**(-1)
                       do j = k + 2, n
                          a( j, k ) = d21*( d11*w( j, k )-w( j, k+1 ) )
                          a( j, k+1 ) = d21*( d22*w( j, k+1 )-w( j, k ) )
                       end do
                    end if
                    ! copy d(k) to a
                    a( k, k ) = w( k, k )
                    a( k+1, k ) = w( k+1, k )
                    a( k+1, k+1 ) = w( k+1, k+1 )
                 end if
              end if
              ! store details of the interchanges in ipiv
              if( kstep==1_${ik}$ ) then
                 ipiv( k ) = kp
              else
                 ipiv( k ) = -kp
                 ipiv( k+1 ) = -kp
              end if
              ! increase k and return to the start of the main loop
              k = k + kstep
              go to 70
              90 continue
              ! update the lower triangle of a22 (= a(k:n,k:n)) as
              ! a22 := a22 - l21*d*l21**t = a22 - l21*w**t
              ! computing blocks of nb columns at a time
              do j = k, n, nb
                 jb = min( nb, n-j+1 )
                 ! update the lower triangle of the diagonal block
                 do jj = j, j + jb - 1
                    call stdlib${ii}$_${ri}$gemv( 'NO TRANSPOSE', j+jb-jj, k-1, -one,a( jj, 1_${ik}$ ), lda, w( jj, &
                              1_${ik}$ ), ldw, one,a( jj, jj ), 1_${ik}$ )
                 end do
                 ! update the rectangular subdiagonal block
                 if( j+jb<=n )call stdlib${ii}$_${ri}$gemm( 'NO TRANSPOSE', 'TRANSPOSE', n-j-jb+1, jb,k-1, -&
                           one, a( j+jb, 1_${ik}$ ), lda, w( j, 1_${ik}$ ), ldw,one, a( j+jb, j ), lda )
              end do
              ! put l21 in standard form by partially undoing the interchanges
              ! of rows in columns 1:k-1 looping backwards from k-1 to 1
              j = k - 1_${ik}$
              120 continue
                 ! undo the interchanges (if any) of rows jj and jp at each
                 ! step j
                 ! (here, j is a diagonal index)
                 jj = j
                 jp = ipiv( j )
                 if( jp<0_${ik}$ ) then
                    jp = -jp
                    ! (here, j is a diagonal index)
                    j = j - 1_${ik}$
                 end if
                 ! (note: here, j is used to determine row length. length j
                 ! of the rows to swap back doesn't include diagonal element)
                 j = j - 1_${ik}$
                 if( jp/=jj .and. j>=1_${ik}$ )call stdlib${ii}$_${ri}$swap( j, a( jp, 1_${ik}$ ), lda, a( jj, 1_${ik}$ ), lda )
                           
              if( j>1 )go to 120
              ! set kb to the number of columns factorized
              kb = k - 1_${ik}$
           end if
           return
     end subroutine stdlib${ii}$_${ri}$lasyf

#:endif
#:endfor

     pure module subroutine stdlib${ii}$_clasyf( uplo, n, nb, kb, a, lda, ipiv, w, ldw, info )
     !! CLASYF computes a partial factorization of a complex symmetric matrix
     !! A using the Bunch-Kaufman diagonal pivoting method. The partial
     !! factorization has the form:
     !! A  =  ( I  U12 ) ( A11  0  ) (  I       0    )  if UPLO = 'U', or:
     !! ( 0  U22 ) (  0   D  ) ( U12**T U22**T )
     !! A  =  ( L11  0 ) ( D    0  ) ( L11**T L21**T )  if UPLO = 'L'
     !! ( L21  I ) ( 0   A22 ) (  0       I    )
     !! where the order of D is at most NB. The actual order is returned in
     !! the argument KB, and is either NB or NB-1, or N if N <= NB.
     !! Note that U**T denotes the transpose of U.
     !! CLASYF is an auxiliary routine called by CSYTRF. It uses blocked code
     !! (calling Level 3 BLAS) to update the submatrix A11 (if UPLO = 'U') or
     !! A22 (if UPLO = 'L').
        ! -- 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, kb
           integer(${ik}$), intent(in) :: lda, ldw, n, nb
           ! Array Arguments 
           integer(${ik}$), intent(out) :: ipiv(*)
           complex(sp), intent(inout) :: a(lda,*)
           complex(sp), intent(out) :: w(ldw,*)
        ! =====================================================================
           ! Parameters 
           real(sp), parameter :: sevten = 17.0e+0_sp
           
           
           
           ! Local Scalars 
           integer(${ik}$) :: imax, j, jb, jj, jmax, jp, k, kk, kkw, kp, kstep, kw
           real(sp) :: absakk, alpha, colmax, rowmax
           complex(sp) :: d11, d21, d22, r1, t, z
           ! Intrinsic Functions 
           ! Statement Functions 
           real(sp) :: cabs1
           ! Statement Function Definitions 
           cabs1( z ) = abs( real( z,KIND=sp) ) + abs( aimag( z ) )
           ! Executable Statements 
           info = 0_${ik}$
           ! initialize alpha for use in choosing pivot block size.
           alpha = ( one+sqrt( sevten ) ) / eight
           if( stdlib_lsame( uplo, 'U' ) ) then
              ! factorize the trailing columns of a using the upper triangle
              ! of a and working backwards, and compute the matrix w = u12*d
              ! for use in updating a11
              ! k is the main loop index, decreasing from n in steps of 1 or 2
              ! kw is the column of w which corresponds to column k of a
              k = n
              10 continue
              kw = nb + k - n
              ! exit from loop
              if( ( k<=n-nb+1 .and. nb<n ) .or. k<1 )go to 30
              ! copy column k of a to column kw of w and update it
              call stdlib${ii}$_ccopy( k, a( 1_${ik}$, k ), 1_${ik}$, w( 1_${ik}$, kw ), 1_${ik}$ )
              if( k<n )call stdlib${ii}$_cgemv( 'NO TRANSPOSE', k, n-k, -cone, a( 1_${ik}$, k+1 ), lda,w( k, &
                        kw+1 ), ldw, cone, w( 1_${ik}$, kw ), 1_${ik}$ )
              kstep = 1_${ik}$
              ! determine rows and columns to be interchanged and whether
              ! a 1-by-1 or 2-by-2 pivot block will be used
              absakk = cabs1( w( k, kw ) )
              ! imax is the row-index of the largest off-diagonal element in
              ! column k, and colmax is its absolute value.
              ! determine both colmax and imax.
              if( k>1_${ik}$ ) then
                 imax = stdlib${ii}$_icamax( k-1, w( 1_${ik}$, kw ), 1_${ik}$ )
                 colmax = cabs1( w( imax, kw ) )
              else
                 colmax = zero
              end if
              if( max( absakk, colmax )==zero ) then
                 ! column k is zero or underflow: set info and continue
                 if( info==0_${ik}$ )info = k
                 kp = k
              else
                 if( absakk>=alpha*colmax ) then
                    ! no interchange, use 1-by-1 pivot block
                    kp = k
                 else
                    ! copy column imax to column kw-1 of w and update it
                    call stdlib${ii}$_ccopy( imax, a( 1_${ik}$, imax ), 1_${ik}$, w( 1_${ik}$, kw-1 ), 1_${ik}$ )
                    call stdlib${ii}$_ccopy( k-imax, a( imax, imax+1 ), lda,w( imax+1, kw-1 ), 1_${ik}$ )
                              
                    if( k<n )call stdlib${ii}$_cgemv( 'NO TRANSPOSE', k, n-k, -cone,a( 1_${ik}$, k+1 ), lda, w(&
                               imax, kw+1 ), ldw,cone, w( 1_${ik}$, kw-1 ), 1_${ik}$ )
                    ! jmax is the column-index of the largest off-diagonal
                    ! element in row imax, and rowmax is its absolute value
                    jmax = imax + stdlib${ii}$_icamax( k-imax, w( imax+1, kw-1 ), 1_${ik}$ )
                    rowmax = cabs1( w( jmax, kw-1 ) )
                    if( imax>1_${ik}$ ) then
                       jmax = stdlib${ii}$_icamax( imax-1, w( 1_${ik}$, kw-1 ), 1_${ik}$ )
                       rowmax = max( rowmax, cabs1( w( jmax, kw-1 ) ) )
                    end if
                    if( absakk>=alpha*colmax*( colmax / rowmax ) ) then
                       ! no interchange, use 1-by-1 pivot block
                       kp = k
                    else if( cabs1( w( imax, kw-1 ) )>=alpha*rowmax ) then
                       ! interchange rows and columns k and imax, use 1-by-1
                       ! pivot block
                       kp = imax
                       ! copy column kw-1 of w to column kw of w
                       call stdlib${ii}$_ccopy( k, w( 1_${ik}$, kw-1 ), 1_${ik}$, w( 1_${ik}$, kw ), 1_${ik}$ )
                    else
                       ! interchange rows and columns k-1 and imax, use 2-by-2
                       ! pivot block
                       kp = imax
                       kstep = 2_${ik}$
                    end if
                 end if
                 ! ============================================================
                 ! kk is the column of a where pivoting step stopped
                 kk = k - kstep + 1_${ik}$
                 ! kkw is the column of w which corresponds to column kk of a
                 kkw = nb + kk - n
                 ! interchange rows and columns kp and kk.
                 ! updated column kp is already stored in column kkw of w.
                 if( kp/=kk ) then
                    ! copy non-updated column kk to column kp of submatrix a
                    ! at step k. no need to copy element into column k
                    ! (or k and k-1 for 2-by-2 pivot) of a, since these columns
                    ! will be later overwritten.
                    a( kp, kp ) = a( kk, kk )
                    call stdlib${ii}$_ccopy( kk-1-kp, a( kp+1, kk ), 1_${ik}$, a( kp, kp+1 ),lda )
                    if( kp>1_${ik}$ )call stdlib${ii}$_ccopy( kp-1, a( 1_${ik}$, kk ), 1_${ik}$, a( 1_${ik}$, kp ), 1_${ik}$ )
                    ! interchange rows kk and kp in last k+1 to n columns of a
                    ! (columns k (or k and k-1 for 2-by-2 pivot) of a will be
                    ! later overwritten). interchange rows kk and kp
                    ! in last kkw to nb columns of w.
                    if( k<n )call stdlib${ii}$_cswap( n-k, a( kk, k+1 ), lda, a( kp, k+1 ),lda )
                    call stdlib${ii}$_cswap( n-kk+1, w( kk, kkw ), ldw, w( kp, kkw ),ldw )
                 end if
                 if( kstep==1_${ik}$ ) then
                    ! 1-by-1 pivot block d(k): column kw of w now holds
                    ! w(kw) = u(k)*d(k),
                    ! where u(k) is the k-th column of u
                    ! store subdiag. elements of column u(k)
                    ! and 1-by-1 block d(k) in column k of a.
                    ! note: diagonal element u(k,k) is a unit element
                    ! and not stored.
                       ! a(k,k) := d(k,k) = w(k,kw)
                       ! a(1:k-1,k) := u(1:k-1,k) = w(1:k-1,kw)/d(k,k)
                    call stdlib${ii}$_ccopy( k, w( 1_${ik}$, kw ), 1_${ik}$, a( 1_${ik}$, k ), 1_${ik}$ )
                    r1 = cone / a( k, k )
                    call stdlib${ii}$_cscal( k-1, r1, a( 1_${ik}$, k ), 1_${ik}$ )
                 else
                    ! 2-by-2 pivot block d(k): columns kw and kw-1 of w now hold
                    ! ( w(kw-1) w(kw) ) = ( u(k-1) u(k) )*d(k)
                    ! where u(k) and u(k-1) are the k-th and (k-1)-th columns
                    ! of u
                    ! store u(1:k-2,k-1) and u(1:k-2,k) and 2-by-2
                    ! block d(k-1:k,k-1:k) in columns k-1 and k of a.
                    ! note: 2-by-2 diagonal block u(k-1:k,k-1:k) is a unit
                    ! block and not stored.
                       ! a(k-1:k,k-1:k) := d(k-1:k,k-1:k) = w(k-1:k,kw-1:kw)
                       ! a(1:k-2,k-1:k) := u(1:k-2,k:k-1:k) =
                       ! = w(1:k-2,kw-1:kw) * ( d(k-1:k,k-1:k)**(-1) )
                    if( k>2_${ik}$ ) then
                       ! compose the columns of the inverse of 2-by-2 pivot
                       ! block d in the following way to reduce the number
                       ! of flops when we myltiply panel ( w(kw-1) w(kw) ) by
                       ! this inverse
                       ! d**(-1) = ( d11 d21 )**(-1) =
                                 ! ( d21 d22 )
                       ! = 1/(d11*d22-d21**2) * ( ( d22 ) (-d21 ) ) =
                                              ! ( (-d21 ) ( d11 ) )
                       ! = 1/d21 * 1/((d11/d21)*(d22/d21)-1) *
                         ! * ( ( d22/d21 ) (      -1 ) ) =
                           ! ( (      -1 ) ( d11/d21 ) )
                       ! = 1/d21 * 1/(d22*d11-1) * ( ( d11 ) (  -1 ) ) =
                                                 ! ( ( -1  ) ( d22 ) )
                       ! = 1/d21 * t * ( ( d11 ) (  -1 ) )
                                     ! ( (  -1 ) ( d22 ) )
                       ! = d21 * ( ( d11 ) (  -1 ) )
                               ! ( (  -1 ) ( d22 ) )
                       d21 = w( k-1, kw )
                       d11 = w( k, kw ) / d21
                       d22 = w( k-1, kw-1 ) / d21
                       t = cone / ( d11*d22-cone )
                       ! update elements in columns a(k-1) and a(k) as
                       ! dot products of rows of ( w(kw-1) w(kw) ) and columns
                       ! of d**(-1)
                       d21 = t / d21
                       do j = 1, k - 2
                          a( j, k-1 ) = d21*( d11*w( j, kw-1 )-w( j, kw ) )
                          a( j, k ) = d21*( d22*w( j, kw )-w( j, kw-1 ) )
                       end do
                    end if
                    ! copy d(k) to a
                    a( k-1, k-1 ) = w( k-1, kw-1 )
                    a( k-1, k ) = w( k-1, kw )
                    a( k, k ) = w( k, kw )
                 end if
              end if
              ! store details of the interchanges in ipiv
              if( kstep==1_${ik}$ ) then
                 ipiv( k ) = kp
              else
                 ipiv( k ) = -kp
                 ipiv( k-1 ) = -kp
              end if
              ! decrease k and return to the start of the main loop
              k = k - kstep
              go to 10
              30 continue
              ! update the upper triangle of a11 (= a(1:k,1:k)) as
              ! a11 := a11 - u12*d*u12**t = a11 - u12*w**t
              ! computing blocks of nb columns at a time
              do j = ( ( k-1 ) / nb )*nb + 1, 1, -nb
                 jb = min( nb, k-j+1 )
                 ! update the upper triangle of the diagonal block
                 do jj = j, j + jb - 1
                    call stdlib${ii}$_cgemv( 'NO TRANSPOSE', jj-j+1, n-k, -cone,a( j, k+1 ), lda, w( jj,&
                               kw+1 ), ldw, cone,a( j, jj ), 1_${ik}$ )
                 end do
                 ! update the rectangular superdiagonal block
                 call stdlib${ii}$_cgemm( 'NO TRANSPOSE', 'TRANSPOSE', j-1, jb, n-k,-cone, a( 1_${ik}$, k+1 ), &
                           lda, w( j, kw+1 ), ldw,cone, a( 1_${ik}$, j ), lda )
              end do
              ! put u12 in standard form by partially undoing the interchanges
              ! in columns k+1:n looping backwards from k+1 to n
              j = k + 1_${ik}$
              60 continue
                 ! undo the interchanges (if any) of rows jj and jp at each
                 ! step j
                 ! (here, j is a diagonal index)
                 jj = j
                 jp = ipiv( j )
                 if( jp<0_${ik}$ ) then
                    jp = -jp
                    ! (here, j is a diagonal index)
                    j = j + 1_${ik}$
                 end if
                 ! (note: here, j is used to determine row length. length n-j+1
                 ! of the rows to swap back doesn't include diagonal element)
                 j = j + 1_${ik}$
                 if( jp/=jj .and. j<=n )call stdlib${ii}$_cswap( n-j+1, a( jp, j ), lda, a( jj, j ), &
                           lda )
              if( j<n )go to 60
              ! set kb to the number of columns factorized
              kb = n - k
           else
              ! factorize the leading columns of a using the lower triangle
              ! of a and working forwards, and compute the matrix w = l21*d
              ! for use in updating a22
              ! k is the main loop index, increasing from 1 in steps of 1 or 2
              k = 1_${ik}$
              70 continue
              ! exit from loop
              if( ( k>=nb .and. nb<n ) .or. k>n )go to 90
              ! copy column k of a to column k of w and update it
              call stdlib${ii}$_ccopy( n-k+1, a( k, k ), 1_${ik}$, w( k, k ), 1_${ik}$ )
              call stdlib${ii}$_cgemv( 'NO TRANSPOSE', n-k+1, k-1, -cone, a( k, 1_${ik}$ ), lda,w( k, 1_${ik}$ ), ldw,&
                         cone, w( k, k ), 1_${ik}$ )
              kstep = 1_${ik}$
              ! determine rows and columns to be interchanged and whether
              ! a 1-by-1 or 2-by-2 pivot block will be used
              absakk = cabs1( w( k, k ) )
              ! imax is the row-index of the largest off-diagonal element in
              ! column k, and colmax is its absolute value.
              ! determine both colmax and imax.
              if( k<n ) then
                 imax = k + stdlib${ii}$_icamax( n-k, w( k+1, k ), 1_${ik}$ )
                 colmax = cabs1( w( imax, k ) )
              else
                 colmax = zero
              end if
              if( max( absakk, colmax )==zero ) then
                 ! column k is zero or underflow: set info and continue
                 if( info==0_${ik}$ )info = k
                 kp = k
              else
                 if( absakk>=alpha*colmax ) then
                    ! no interchange, use 1-by-1 pivot block
                    kp = k
                 else
                    ! copy column imax to column k+1 of w and update it
                    call stdlib${ii}$_ccopy( imax-k, a( imax, k ), lda, w( k, k+1 ), 1_${ik}$ )
                    call stdlib${ii}$_ccopy( n-imax+1, a( imax, imax ), 1_${ik}$, w( imax, k+1 ),1_${ik}$ )
                    call stdlib${ii}$_cgemv( 'NO TRANSPOSE', n-k+1, k-1, -cone, a( k, 1_${ik}$ ),lda, w( imax, &
                              1_${ik}$ ), ldw, cone, w( k, k+1 ),1_${ik}$ )
                    ! jmax is the column-index of the largest off-diagonal
                    ! element in row imax, and rowmax is its absolute value
                    jmax = k - 1_${ik}$ + stdlib${ii}$_icamax( imax-k, w( k, k+1 ), 1_${ik}$ )
                    rowmax = cabs1( w( jmax, k+1 ) )
                    if( imax<n ) then
                       jmax = imax + stdlib${ii}$_icamax( n-imax, w( imax+1, k+1 ), 1_${ik}$ )
                       rowmax = max( rowmax, cabs1( w( jmax, k+1 ) ) )
                    end if
                    if( absakk>=alpha*colmax*( colmax / rowmax ) ) then
                       ! no interchange, use 1-by-1 pivot block
                       kp = k
                    else if( cabs1( w( imax, k+1 ) )>=alpha*rowmax ) then
                       ! interchange rows and columns k and imax, use 1-by-1
                       ! pivot block
                       kp = imax
                       ! copy column k+1 of w to column k of w
                       call stdlib${ii}$_ccopy( n-k+1, w( k, k+1 ), 1_${ik}$, w( k, k ), 1_${ik}$ )
                    else
                       ! interchange rows and columns k+1 and imax, use 2-by-2
                       ! pivot block
                       kp = imax
                       kstep = 2_${ik}$
                    end if
                 end if
                 ! ============================================================
                 ! kk is the column of a where pivoting step stopped
                 kk = k + kstep - 1_${ik}$
                 ! interchange rows and columns kp and kk.
                 ! updated column kp is already stored in column kk of w.
                 if( kp/=kk ) then
                    ! copy non-updated column kk to column kp of submatrix a
                    ! at step k. no need to copy element into column k
                    ! (or k and k+1 for 2-by-2 pivot) of a, since these columns
                    ! will be later overwritten.
                    a( kp, kp ) = a( kk, kk )
                    call stdlib${ii}$_ccopy( kp-kk-1, a( kk+1, kk ), 1_${ik}$, a( kp, kk+1 ),lda )
                    if( kp<n )call stdlib${ii}$_ccopy( n-kp, a( kp+1, kk ), 1_${ik}$, a( kp+1, kp ), 1_${ik}$ )
                              
                    ! interchange rows kk and kp in first k-1 columns of a
                    ! (columns k (or k and k+1 for 2-by-2 pivot) of a will be
                    ! later overwritten). interchange rows kk and kp
                    ! in first kk columns of w.
                    if( k>1_${ik}$ )call stdlib${ii}$_cswap( k-1, a( kk, 1_${ik}$ ), lda, a( kp, 1_${ik}$ ), lda )
                    call stdlib${ii}$_cswap( kk, w( kk, 1_${ik}$ ), ldw, w( kp, 1_${ik}$ ), ldw )
                 end if
                 if( kstep==1_${ik}$ ) then
                    ! 1-by-1 pivot block d(k): column k of w now holds
                    ! w(k) = l(k)*d(k),
                    ! where l(k) is the k-th column of l
                    ! store subdiag. elements of column l(k)
                    ! and 1-by-1 block d(k) in column k of a.
                    ! (note: diagonal element l(k,k) is a unit element
                    ! and not stored)
                       ! a(k,k) := d(k,k) = w(k,k)
                       ! a(k+1:n,k) := l(k+1:n,k) = w(k+1:n,k)/d(k,k)
                    call stdlib${ii}$_ccopy( n-k+1, w( k, k ), 1_${ik}$, a( k, k ), 1_${ik}$ )
                    if( k<n ) then
                       r1 = cone / a( k, k )
                       call stdlib${ii}$_cscal( n-k, r1, a( k+1, k ), 1_${ik}$ )
                    end if
                 else
                    ! 2-by-2 pivot block d(k): columns k and k+1 of w now hold
                    ! ( w(k) w(k+1) ) = ( l(k) l(k+1) )*d(k)
                    ! where l(k) and l(k+1) are the k-th and (k+1)-th columns
                    ! of l
                    ! store l(k+2:n,k) and l(k+2:n,k+1) and 2-by-2
                    ! block d(k:k+1,k:k+1) in columns k and k+1 of a.
                    ! (note: 2-by-2 diagonal block l(k:k+1,k:k+1) is a unit
                    ! block and not stored)
                       ! a(k:k+1,k:k+1) := d(k:k+1,k:k+1) = w(k:k+1,k:k+1)
                       ! a(k+2:n,k:k+1) := l(k+2:n,k:k+1) =
                       ! = w(k+2:n,k:k+1) * ( d(k:k+1,k:k+1)**(-1) )
                    if( k<n-1 ) then
                       ! compose the columns of the inverse of 2-by-2 pivot
                       ! block d in the following way to reduce the number
                       ! of flops when we myltiply panel ( w(k) w(k+1) ) by
                       ! this inverse
                       ! d**(-1) = ( d11 d21 )**(-1) =
                                 ! ( d21 d22 )
                       ! = 1/(d11*d22-d21**2) * ( ( d22 ) (-d21 ) ) =
                                              ! ( (-d21 ) ( d11 ) )
                       ! = 1/d21 * 1/((d11/d21)*(d22/d21)-1) *
                         ! * ( ( d22/d21 ) (      -1 ) ) =
                           ! ( (      -1 ) ( d11/d21 ) )
                       ! = 1/d21 * 1/(d22*d11-1) * ( ( d11 ) (  -1 ) ) =
                                                 ! ( ( -1  ) ( d22 ) )
                       ! = 1/d21 * t * ( ( d11 ) (  -1 ) )
                                     ! ( (  -1 ) ( d22 ) )
                       ! = d21 * ( ( d11 ) (  -1 ) )
                               ! ( (  -1 ) ( d22 ) )
                       d21 = w( k+1, k )
                       d11 = w( k+1, k+1 ) / d21
                       d22 = w( k, k ) / d21
                       t = cone / ( d11*d22-cone )
                       d21 = t / d21
                       ! update elements in columns a(k) and a(k+1) as
                       ! dot products of rows of ( w(k) w(k+1) ) and columns
                       ! of d**(-1)
                       do j = k + 2, n
                          a( j, k ) = d21*( d11*w( j, k )-w( j, k+1 ) )
                          a( j, k+1 ) = d21*( d22*w( j, k+1 )-w( j, k ) )
                       end do
                    end if
                    ! copy d(k) to a
                    a( k, k ) = w( k, k )
                    a( k+1, k ) = w( k+1, k )
                    a( k+1, k+1 ) = w( k+1, k+1 )
                 end if
              end if
              ! store details of the interchanges in ipiv
              if( kstep==1_${ik}$ ) then
                 ipiv( k ) = kp
              else
                 ipiv( k ) = -kp
                 ipiv( k+1 ) = -kp
              end if
              ! increase k and return to the start of the main loop
              k = k + kstep
              go to 70
              90 continue
              ! update the lower triangle of a22 (= a(k:n,k:n)) as
              ! a22 := a22 - l21*d*l21**t = a22 - l21*w**t
              ! computing blocks of nb columns at a time
              do j = k, n, nb
                 jb = min( nb, n-j+1 )
                 ! update the lower triangle of the diagonal block
                 do jj = j, j + jb - 1
                    call stdlib${ii}$_cgemv( 'NO TRANSPOSE', j+jb-jj, k-1, -cone,a( jj, 1_${ik}$ ), lda, w( jj,&
                               1_${ik}$ ), ldw, cone,a( jj, jj ), 1_${ik}$ )
                 end do
                 ! update the rectangular subdiagonal block
                 if( j+jb<=n )call stdlib${ii}$_cgemm( 'NO TRANSPOSE', 'TRANSPOSE', n-j-jb+1, jb,k-1, -&
                           cone, a( j+jb, 1_${ik}$ ), lda, w( j, 1_${ik}$ ),ldw, cone, a( j+jb, j ), lda )
              end do
              ! put l21 in standard form by partially undoing the interchanges
              ! of rows in columns 1:k-1 looping backwards from k-1 to 1
              j = k - 1_${ik}$
              120 continue
                 ! undo the interchanges (if any) of rows jj and jp at each
                 ! step j
                 ! (here, j is a diagonal index)
                 jj = j
                 jp = ipiv( j )
                 if( jp<0_${ik}$ ) then
                    jp = -jp
                    ! (here, j is a diagonal index)
                    j = j - 1_${ik}$
                 end if
                 ! (note: here, j is used to determine row length. length j
                 ! of the rows to swap back doesn't include diagonal element)
                 j = j - 1_${ik}$
                 if( jp/=jj .and. j>=1_${ik}$ )call stdlib${ii}$_cswap( j, a( jp, 1_${ik}$ ), lda, a( jj, 1_${ik}$ ), lda )
                           
              if( j>1 )go to 120
              ! set kb to the number of columns factorized
              kb = k - 1_${ik}$
           end if
           return
     end subroutine stdlib${ii}$_clasyf

     pure module subroutine stdlib${ii}$_zlasyf( uplo, n, nb, kb, a, lda, ipiv, w, ldw, info )
     !! ZLASYF computes a partial factorization of a complex symmetric matrix
     !! A using the Bunch-Kaufman diagonal pivoting method. The partial
     !! factorization has the form:
     !! A  =  ( I  U12 ) ( A11  0  ) (  I       0    )  if UPLO = 'U', or:
     !! ( 0  U22 ) (  0   D  ) ( U12**T U22**T )
     !! A  =  ( L11  0 ) ( D    0  ) ( L11**T L21**T )  if UPLO = 'L'
     !! ( L21  I ) ( 0   A22 ) (  0       I    )
     !! where the order of D is at most NB. The actual order is returned in
     !! the argument KB, and is either NB or NB-1, or N if N <= NB.
     !! Note that U**T denotes the transpose of U.
     !! ZLASYF is an auxiliary routine called by ZSYTRF. It uses blocked code
     !! (calling Level 3 BLAS) to update the submatrix A11 (if UPLO = 'U') or
     !! A22 (if UPLO = 'L').
        ! -- 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, kb
           integer(${ik}$), intent(in) :: lda, ldw, n, nb
           ! Array Arguments 
           integer(${ik}$), intent(out) :: ipiv(*)
           complex(dp), intent(inout) :: a(lda,*)
           complex(dp), intent(out) :: w(ldw,*)
        ! =====================================================================
           ! Parameters 
           real(dp), parameter :: sevten = 17.0e+0_dp
           
           
           
           ! Local Scalars 
           integer(${ik}$) :: imax, j