stdlib_lapack_solve_ldl_comp2.fypp Source File


Source Code

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


  contains
#:for ik,it,ii in LINALG_INT_KINDS_TYPES

     pure module subroutine stdlib${ii}$_ssptrs( uplo, n, nrhs, ap, ipiv, b, ldb, info )
     !! SSPTRS solves a system of linear equations A*X = B with a real
     !! symmetric matrix A stored in packed format using the factorization
     !! A = U*D*U**T or A = L*D*L**T computed by SSPTRF.
        ! -- 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) :: ldb, n, nrhs
           ! Array Arguments 
           integer(${ik}$), intent(in) :: ipiv(*)
           real(sp), intent(in) :: ap(*)
           real(sp), intent(inout) :: b(ldb,*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: upper
           integer(${ik}$) :: j, k, kc, kp
           real(sp) :: ak, akm1, akm1k, bk, bkm1, denom
           ! Intrinsic Functions 
           ! Executable Statements 
           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( nrhs<0_${ik}$ ) then
              info = -3_${ik}$
           else if( ldb<max( 1_${ik}$, n ) ) then
              info = -7_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'SSPTRS', -info )
              return
           end if
           ! quick return if possible
           if( n==0 .or. nrhs==0 )return
           if( upper ) then
              ! solve a*x = b, where a = u*d*u**t.
              ! first solve u*d*x = b, overwriting b with x.
              ! k is the main loop index, decreasing from n to 1 in steps of
              ! 1 or 2, depending on the size of the diagonal blocks.
              k = n
              kc = n*( n+1 ) / 2_${ik}$ + 1_${ik}$
              10 continue
              ! if k < 1, exit from loop.
              if( k<1 )go to 30
              kc = kc - k
              if( ipiv( k )>0_${ik}$ ) then
                 ! 1 x 1 diagonal block
                 ! interchange rows k and ipiv(k).
                 kp = ipiv( k )
                 if( kp/=k )call stdlib${ii}$_sswap( nrhs, b( k, 1_${ik}$ ), ldb, b( kp, 1_${ik}$ ), ldb )
                 ! multiply by inv(u(k)), where u(k) is the transformation
                 ! stored in column k of a.
                 call stdlib${ii}$_sger( k-1, nrhs, -one, ap( kc ), 1_${ik}$, b( k, 1_${ik}$ ), ldb,b( 1_${ik}$, 1_${ik}$ ), ldb )
                           
                 ! multiply by the inverse of the diagonal block.
                 call stdlib${ii}$_sscal( nrhs, one / ap( kc+k-1 ), b( k, 1_${ik}$ ), ldb )
                 k = k - 1_${ik}$
              else
                 ! 2 x 2 diagonal block
                 ! interchange rows k-1 and -ipiv(k).
                 kp = -ipiv( k )
                 if( kp/=k-1 )call stdlib${ii}$_sswap( nrhs, b( k-1, 1_${ik}$ ), ldb, b( kp, 1_${ik}$ ), ldb )
                 ! multiply by inv(u(k)), where u(k) is the transformation
                 ! stored in columns k-1 and k of a.
                 call stdlib${ii}$_sger( k-2, nrhs, -one, ap( kc ), 1_${ik}$, b( k, 1_${ik}$ ), ldb,b( 1_${ik}$, 1_${ik}$ ), ldb )
                           
                 call stdlib${ii}$_sger( k-2, nrhs, -one, ap( kc-( k-1 ) ), 1_${ik}$,b( k-1, 1_${ik}$ ), ldb, b( 1_${ik}$, 1_${ik}$ &
                           ), ldb )
                 ! multiply by the inverse of the diagonal block.
                 akm1k = ap( kc+k-2 )
                 akm1 = ap( kc-1 ) / akm1k
                 ak = ap( kc+k-1 ) / akm1k
                 denom = akm1*ak - one
                 do j = 1, nrhs
                    bkm1 = b( k-1, j ) / akm1k
                    bk = b( k, j ) / akm1k
                    b( k-1, j ) = ( ak*bkm1-bk ) / denom
                    b( k, j ) = ( akm1*bk-bkm1 ) / denom
                 end do
                 kc = kc - k + 1_${ik}$
                 k = k - 2_${ik}$
              end if
              go to 10
              30 continue
              ! next solve u**t*x = b, overwriting b with x.
              ! k is the main loop index, increasing from 1 to n in steps of
              ! 1 or 2, depending on the size of the diagonal blocks.
              k = 1_${ik}$
              kc = 1_${ik}$
              40 continue
              ! if k > n, exit from loop.
              if( k>n )go to 50
              if( ipiv( k )>0_${ik}$ ) then
                 ! 1 x 1 diagonal block
                 ! multiply by inv(u**t(k)), where u(k) is the transformation
                 ! stored in column k of a.
                 call stdlib${ii}$_sgemv( 'TRANSPOSE', k-1, nrhs, -one, b, ldb, ap( kc ),1_${ik}$, one, b( k, &
                           1_${ik}$ ), ldb )
                 ! interchange rows k and ipiv(k).
                 kp = ipiv( k )
                 if( kp/=k )call stdlib${ii}$_sswap( nrhs, b( k, 1_${ik}$ ), ldb, b( kp, 1_${ik}$ ), ldb )
                 kc = kc + k
                 k = k + 1_${ik}$
              else
                 ! 2 x 2 diagonal block
                 ! multiply by inv(u**t(k+1)), where u(k+1) is the transformation
                 ! stored in columns k and k+1 of a.
                 call stdlib${ii}$_sgemv( 'TRANSPOSE', k-1, nrhs, -one, b, ldb, ap( kc ),1_${ik}$, one, b( k, &
                           1_${ik}$ ), ldb )
                 call stdlib${ii}$_sgemv( 'TRANSPOSE', k-1, nrhs, -one, b, ldb,ap( kc+k ), 1_${ik}$, one, b( k+&
                           1_${ik}$, 1_${ik}$ ), ldb )
                 ! interchange rows k and -ipiv(k).
                 kp = -ipiv( k )
                 if( kp/=k )call stdlib${ii}$_sswap( nrhs, b( k, 1_${ik}$ ), ldb, b( kp, 1_${ik}$ ), ldb )
                 kc = kc + 2_${ik}$*k + 1_${ik}$
                 k = k + 2_${ik}$
              end if
              go to 40
              50 continue
           else
              ! solve a*x = b, where a = l*d*l**t.
              ! first solve l*d*x = b, overwriting b with x.
              ! k is the main loop index, increasing from 1 to n in steps of
              ! 1 or 2, depending on the size of the diagonal blocks.
              k = 1_${ik}$
              kc = 1_${ik}$
              60 continue
              ! if k > n, exit from loop.
              if( k>n )go to 80
              if( ipiv( k )>0_${ik}$ ) then
                 ! 1 x 1 diagonal block
                 ! interchange rows k and ipiv(k).
                 kp = ipiv( k )
                 if( kp/=k )call stdlib${ii}$_sswap( nrhs, b( k, 1_${ik}$ ), ldb, b( kp, 1_${ik}$ ), ldb )
                 ! multiply by inv(l(k)), where l(k) is the transformation
                 ! stored in column k of a.
                 if( k<n )call stdlib${ii}$_sger( n-k, nrhs, -one, ap( kc+1 ), 1_${ik}$, b( k, 1_${ik}$ ),ldb, b( k+1,&
                            1_${ik}$ ), ldb )
                 ! multiply by the inverse of the diagonal block.
                 call stdlib${ii}$_sscal( nrhs, one / ap( kc ), b( k, 1_${ik}$ ), ldb )
                 kc = kc + n - k + 1_${ik}$
                 k = k + 1_${ik}$
              else
                 ! 2 x 2 diagonal block
                 ! interchange rows k+1 and -ipiv(k).
                 kp = -ipiv( k )
                 if( kp/=k+1 )call stdlib${ii}$_sswap( nrhs, b( k+1, 1_${ik}$ ), ldb, b( kp, 1_${ik}$ ), ldb )
                 ! multiply by inv(l(k)), where l(k) is the transformation
                 ! stored in columns k and k+1 of a.
                 if( k<n-1 ) then
                    call stdlib${ii}$_sger( n-k-1, nrhs, -one, ap( kc+2 ), 1_${ik}$, b( k, 1_${ik}$ ),ldb, b( k+2, 1_${ik}$ )&
                              , ldb )
                    call stdlib${ii}$_sger( n-k-1, nrhs, -one, ap( kc+n-k+2 ), 1_${ik}$,b( k+1, 1_${ik}$ ), ldb, b( k+&
                              2_${ik}$, 1_${ik}$ ), ldb )
                 end if
                 ! multiply by the inverse of the diagonal block.
                 akm1k = ap( kc+1 )
                 akm1 = ap( kc ) / akm1k
                 ak = ap( kc+n-k+1 ) / akm1k
                 denom = akm1*ak - one
                 do j = 1, nrhs
                    bkm1 = b( k, j ) / akm1k
                    bk = b( k+1, j ) / akm1k
                    b( k, j ) = ( ak*bkm1-bk ) / denom
                    b( k+1, j ) = ( akm1*bk-bkm1 ) / denom
                 end do
                 kc = kc + 2_${ik}$*( n-k ) + 1_${ik}$
                 k = k + 2_${ik}$
              end if
              go to 60
              80 continue
              ! next solve l**t*x = b, overwriting b with x.
              ! k is the main loop index, decreasing from n to 1 in steps of
              ! 1 or 2, depending on the size of the diagonal blocks.
              k = n
              kc = n*( n+1 ) / 2_${ik}$ + 1_${ik}$
              90 continue
              ! if k < 1, exit from loop.
              if( k<1 )go to 100
              kc = kc - ( n-k+1 )
              if( ipiv( k )>0_${ik}$ ) then
                 ! 1 x 1 diagonal block
                 ! multiply by inv(l**t(k)), where l(k) is the transformation
                 ! stored in column k of a.
                 if( k<n )call stdlib${ii}$_sgemv( 'TRANSPOSE', n-k, nrhs, -one, b( k+1, 1_${ik}$ ),ldb, ap( &
                           kc+1 ), 1_${ik}$, one, b( k, 1_${ik}$ ), ldb )
                 ! interchange rows k and ipiv(k).
                 kp = ipiv( k )
                 if( kp/=k )call stdlib${ii}$_sswap( nrhs, b( k, 1_${ik}$ ), ldb, b( kp, 1_${ik}$ ), ldb )
                 k = k - 1_${ik}$
              else
                 ! 2 x 2 diagonal block
                 ! multiply by inv(l**t(k-1)), where l(k-1) is the transformation
                 ! stored in columns k-1 and k of a.
                 if( k<n ) then
                    call stdlib${ii}$_sgemv( 'TRANSPOSE', n-k, nrhs, -one, b( k+1, 1_${ik}$ ),ldb, ap( kc+1 ), &
                              1_${ik}$, one, b( k, 1_${ik}$ ), ldb )
                    call stdlib${ii}$_sgemv( 'TRANSPOSE', n-k, nrhs, -one, b( k+1, 1_${ik}$ ),ldb, ap( kc-( n-&
                              k ) ), 1_${ik}$, one, b( k-1, 1_${ik}$ ),ldb )
                 end if
                 ! interchange rows k and -ipiv(k).
                 kp = -ipiv( k )
                 if( kp/=k )call stdlib${ii}$_sswap( nrhs, b( k, 1_${ik}$ ), ldb, b( kp, 1_${ik}$ ), ldb )
                 kc = kc - ( n-k+2 )
                 k = k - 2_${ik}$
              end if
              go to 90
              100 continue
           end if
           return
     end subroutine stdlib${ii}$_ssptrs

     pure module subroutine stdlib${ii}$_dsptrs( uplo, n, nrhs, ap, ipiv, b, ldb, info )
     !! DSPTRS solves a system of linear equations A*X = B with a real
     !! symmetric matrix A stored in packed format using the factorization
     !! A = U*D*U**T or A = L*D*L**T computed by DSPTRF.
        ! -- 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) :: ldb, n, nrhs
           ! Array Arguments 
           integer(${ik}$), intent(in) :: ipiv(*)
           real(dp), intent(in) :: ap(*)
           real(dp), intent(inout) :: b(ldb,*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: upper
           integer(${ik}$) :: j, k, kc, kp
           real(dp) :: ak, akm1, akm1k, bk, bkm1, denom
           ! Intrinsic Functions 
           ! Executable Statements 
           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( nrhs<0_${ik}$ ) then
              info = -3_${ik}$
           else if( ldb<max( 1_${ik}$, n ) ) then
              info = -7_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DSPTRS', -info )
              return
           end if
           ! quick return if possible
           if( n==0 .or. nrhs==0 )return
           if( upper ) then
              ! solve a*x = b, where a = u*d*u**t.
              ! first solve u*d*x = b, overwriting b with x.
              ! k is the main loop index, decreasing from n to 1 in steps of
              ! 1 or 2, depending on the size of the diagonal blocks.
              k = n
              kc = n*( n+1 ) / 2_${ik}$ + 1_${ik}$
              10 continue
              ! if k < 1, exit from loop.
              if( k<1 )go to 30
              kc = kc - k
              if( ipiv( k )>0_${ik}$ ) then
                 ! 1 x 1 diagonal block
                 ! interchange rows k and ipiv(k).
                 kp = ipiv( k )
                 if( kp/=k )call stdlib${ii}$_dswap( nrhs, b( k, 1_${ik}$ ), ldb, b( kp, 1_${ik}$ ), ldb )
                 ! multiply by inv(u(k)), where u(k) is the transformation
                 ! stored in column k of a.
                 call stdlib${ii}$_dger( k-1, nrhs, -one, ap( kc ), 1_${ik}$, b( k, 1_${ik}$ ), ldb,b( 1_${ik}$, 1_${ik}$ ), ldb )
                           
                 ! multiply by the inverse of the diagonal block.
                 call stdlib${ii}$_dscal( nrhs, one / ap( kc+k-1 ), b( k, 1_${ik}$ ), ldb )
                 k = k - 1_${ik}$
              else
                 ! 2 x 2 diagonal block
                 ! interchange rows k-1 and -ipiv(k).
                 kp = -ipiv( k )
                 if( kp/=k-1 )call stdlib${ii}$_dswap( nrhs, b( k-1, 1_${ik}$ ), ldb, b( kp, 1_${ik}$ ), ldb )
                 ! multiply by inv(u(k)), where u(k) is the transformation
                 ! stored in columns k-1 and k of a.
                 call stdlib${ii}$_dger( k-2, nrhs, -one, ap( kc ), 1_${ik}$, b( k, 1_${ik}$ ), ldb,b( 1_${ik}$, 1_${ik}$ ), ldb )
                           
                 call stdlib${ii}$_dger( k-2, nrhs, -one, ap( kc-( k-1 ) ), 1_${ik}$,b( k-1, 1_${ik}$ ), ldb, b( 1_${ik}$, 1_${ik}$ &
                           ), ldb )
                 ! multiply by the inverse of the diagonal block.
                 akm1k = ap( kc+k-2 )
                 akm1 = ap( kc-1 ) / akm1k
                 ak = ap( kc+k-1 ) / akm1k
                 denom = akm1*ak - one
                 do j = 1, nrhs
                    bkm1 = b( k-1, j ) / akm1k
                    bk = b( k, j ) / akm1k
                    b( k-1, j ) = ( ak*bkm1-bk ) / denom
                    b( k, j ) = ( akm1*bk-bkm1 ) / denom
                 end do
                 kc = kc - k + 1_${ik}$
                 k = k - 2_${ik}$
              end if
              go to 10
              30 continue
              ! next solve u**t*x = b, overwriting b with x.
              ! k is the main loop index, increasing from 1 to n in steps of
              ! 1 or 2, depending on the size of the diagonal blocks.
              k = 1_${ik}$
              kc = 1_${ik}$
              40 continue
              ! if k > n, exit from loop.
              if( k>n )go to 50
              if( ipiv( k )>0_${ik}$ ) then
                 ! 1 x 1 diagonal block
                 ! multiply by inv(u**t(k)), where u(k) is the transformation
                 ! stored in column k of a.
                 call stdlib${ii}$_dgemv( 'TRANSPOSE', k-1, nrhs, -one, b, ldb, ap( kc ),1_${ik}$, one, b( k, &
                           1_${ik}$ ), ldb )
                 ! interchange rows k and ipiv(k).
                 kp = ipiv( k )
                 if( kp/=k )call stdlib${ii}$_dswap( nrhs, b( k, 1_${ik}$ ), ldb, b( kp, 1_${ik}$ ), ldb )
                 kc = kc + k
                 k = k + 1_${ik}$
              else
                 ! 2 x 2 diagonal block
                 ! multiply by inv(u**t(k+1)), where u(k+1) is the transformation
                 ! stored in columns k and k+1 of a.
                 call stdlib${ii}$_dgemv( 'TRANSPOSE', k-1, nrhs, -one, b, ldb, ap( kc ),1_${ik}$, one, b( k, &
                           1_${ik}$ ), ldb )
                 call stdlib${ii}$_dgemv( 'TRANSPOSE', k-1, nrhs, -one, b, ldb,ap( kc+k ), 1_${ik}$, one, b( k+&
                           1_${ik}$, 1_${ik}$ ), ldb )
                 ! interchange rows k and -ipiv(k).
                 kp = -ipiv( k )
                 if( kp/=k )call stdlib${ii}$_dswap( nrhs, b( k, 1_${ik}$ ), ldb, b( kp, 1_${ik}$ ), ldb )
                 kc = kc + 2_${ik}$*k + 1_${ik}$
                 k = k + 2_${ik}$
              end if
              go to 40
              50 continue
           else
              ! solve a*x = b, where a = l*d*l**t.
              ! first solve l*d*x = b, overwriting b with x.
              ! k is the main loop index, increasing from 1 to n in steps of
              ! 1 or 2, depending on the size of the diagonal blocks.
              k = 1_${ik}$
              kc = 1_${ik}$
              60 continue
              ! if k > n, exit from loop.
              if( k>n )go to 80
              if( ipiv( k )>0_${ik}$ ) then
                 ! 1 x 1 diagonal block
                 ! interchange rows k and ipiv(k).
                 kp = ipiv( k )
                 if( kp/=k )call stdlib${ii}$_dswap( nrhs, b( k, 1_${ik}$ ), ldb, b( kp, 1_${ik}$ ), ldb )
                 ! multiply by inv(l(k)), where l(k) is the transformation
                 ! stored in column k of a.
                 if( k<n )call stdlib${ii}$_dger( n-k, nrhs, -one, ap( kc+1 ), 1_${ik}$, b( k, 1_${ik}$ ),ldb, b( k+1,&
                            1_${ik}$ ), ldb )
                 ! multiply by the inverse of the diagonal block.
                 call stdlib${ii}$_dscal( nrhs, one / ap( kc ), b( k, 1_${ik}$ ), ldb )
                 kc = kc + n - k + 1_${ik}$
                 k = k + 1_${ik}$
              else
                 ! 2 x 2 diagonal block
                 ! interchange rows k+1 and -ipiv(k).
                 kp = -ipiv( k )
                 if( kp/=k+1 )call stdlib${ii}$_dswap( nrhs, b( k+1, 1_${ik}$ ), ldb, b( kp, 1_${ik}$ ), ldb )
                 ! multiply by inv(l(k)), where l(k) is the transformation
                 ! stored in columns k and k+1 of a.
                 if( k<n-1 ) then
                    call stdlib${ii}$_dger( n-k-1, nrhs, -one, ap( kc+2 ), 1_${ik}$, b( k, 1_${ik}$ ),ldb, b( k+2, 1_${ik}$ )&
                              , ldb )
                    call stdlib${ii}$_dger( n-k-1, nrhs, -one, ap( kc+n-k+2 ), 1_${ik}$,b( k+1, 1_${ik}$ ), ldb, b( k+&
                              2_${ik}$, 1_${ik}$ ), ldb )
                 end if
                 ! multiply by the inverse of the diagonal block.
                 akm1k = ap( kc+1 )
                 akm1 = ap( kc ) / akm1k
                 ak = ap( kc+n-k+1 ) / akm1k
                 denom = akm1*ak - one
                 do j = 1, nrhs
                    bkm1 = b( k, j ) / akm1k
                    bk = b( k+1, j ) / akm1k
                    b( k, j ) = ( ak*bkm1-bk ) / denom
                    b( k+1, j ) = ( akm1*bk-bkm1 ) / denom
                 end do
                 kc = kc + 2_${ik}$*( n-k ) + 1_${ik}$
                 k = k + 2_${ik}$
              end if
              go to 60
              80 continue
              ! next solve l**t*x = b, overwriting b with x.
              ! k is the main loop index, decreasing from n to 1 in steps of
              ! 1 or 2, depending on the size of the diagonal blocks.
              k = n
              kc = n*( n+1 ) / 2_${ik}$ + 1_${ik}$
              90 continue
              ! if k < 1, exit from loop.
              if( k<1 )go to 100
              kc = kc - ( n-k+1 )
              if( ipiv( k )>0_${ik}$ ) then
                 ! 1 x 1 diagonal block
                 ! multiply by inv(l**t(k)), where l(k) is the transformation
                 ! stored in column k of a.
                 if( k<n )call stdlib${ii}$_dgemv( 'TRANSPOSE', n-k, nrhs, -one, b( k+1, 1_${ik}$ ),ldb, ap( &
                           kc+1 ), 1_${ik}$, one, b( k, 1_${ik}$ ), ldb )
                 ! interchange rows k and ipiv(k).
                 kp = ipiv( k )
                 if( kp/=k )call stdlib${ii}$_dswap( nrhs, b( k, 1_${ik}$ ), ldb, b( kp, 1_${ik}$ ), ldb )
                 k = k - 1_${ik}$
              else
                 ! 2 x 2 diagonal block
                 ! multiply by inv(l**t(k-1)), where l(k-1) is the transformation
                 ! stored in columns k-1 and k of a.
                 if( k<n ) then
                    call stdlib${ii}$_dgemv( 'TRANSPOSE', n-k, nrhs, -one, b( k+1, 1_${ik}$ ),ldb, ap( kc+1 ), &
                              1_${ik}$, one, b( k, 1_${ik}$ ), ldb )
                    call stdlib${ii}$_dgemv( 'TRANSPOSE', n-k, nrhs, -one, b( k+1, 1_${ik}$ ),ldb, ap( kc-( n-&
                              k ) ), 1_${ik}$, one, b( k-1, 1_${ik}$ ),ldb )
                 end if
                 ! interchange rows k and -ipiv(k).
                 kp = -ipiv( k )
                 if( kp/=k )call stdlib${ii}$_dswap( nrhs, b( k, 1_${ik}$ ), ldb, b( kp, 1_${ik}$ ), ldb )
                 kc = kc - ( n-k+2 )
                 k = k - 2_${ik}$
              end if
              go to 90
              100 continue
           end if
           return
     end subroutine stdlib${ii}$_dsptrs

#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
     pure module subroutine stdlib${ii}$_${ri}$sptrs( uplo, n, nrhs, ap, ipiv, b, ldb, info )
     !! DSPTRS: solves a system of linear equations A*X = B with a real
     !! symmetric matrix A stored in packed format using the factorization
     !! A = U*D*U**T or A = L*D*L**T computed by DSPTRF.
        ! -- 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) :: ldb, n, nrhs
           ! Array Arguments 
           integer(${ik}$), intent(in) :: ipiv(*)
           real(${rk}$), intent(in) :: ap(*)
           real(${rk}$), intent(inout) :: b(ldb,*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: upper
           integer(${ik}$) :: j, k, kc, kp
           real(${rk}$) :: ak, akm1, akm1k, bk, bkm1, denom
           ! Intrinsic Functions 
           ! Executable Statements 
           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( nrhs<0_${ik}$ ) then
              info = -3_${ik}$
           else if( ldb<max( 1_${ik}$, n ) ) then
              info = -7_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DSPTRS', -info )
              return
           end if
           ! quick return if possible
           if( n==0 .or. nrhs==0 )return
           if( upper ) then
              ! solve a*x = b, where a = u*d*u**t.
              ! first solve u*d*x = b, overwriting b with x.
              ! k is the main loop index, decreasing from n to 1 in steps of
              ! 1 or 2, depending on the size of the diagonal blocks.
              k = n
              kc = n*( n+1 ) / 2_${ik}$ + 1_${ik}$
              10 continue
              ! if k < 1, exit from loop.
              if( k<1 )go to 30
              kc = kc - k
              if( ipiv( k )>0_${ik}$ ) then
                 ! 1 x 1 diagonal block
                 ! interchange rows k and ipiv(k).
                 kp = ipiv( k )
                 if( kp/=k )call stdlib${ii}$_${ri}$swap( nrhs, b( k, 1_${ik}$ ), ldb, b( kp, 1_${ik}$ ), ldb )
                 ! multiply by inv(u(k)), where u(k) is the transformation
                 ! stored in column k of a.
                 call stdlib${ii}$_${ri}$ger( k-1, nrhs, -one, ap( kc ), 1_${ik}$, b( k, 1_${ik}$ ), ldb,b( 1_${ik}$, 1_${ik}$ ), ldb )
                           
                 ! multiply by the inverse of the diagonal block.
                 call stdlib${ii}$_${ri}$scal( nrhs, one / ap( kc+k-1 ), b( k, 1_${ik}$ ), ldb )
                 k = k - 1_${ik}$
              else
                 ! 2 x 2 diagonal block
                 ! interchange rows k-1 and -ipiv(k).
                 kp = -ipiv( k )
                 if( kp/=k-1 )call stdlib${ii}$_${ri}$swap( nrhs, b( k-1, 1_${ik}$ ), ldb, b( kp, 1_${ik}$ ), ldb )
                 ! multiply by inv(u(k)), where u(k) is the transformation
                 ! stored in columns k-1 and k of a.
                 call stdlib${ii}$_${ri}$ger( k-2, nrhs, -one, ap( kc ), 1_${ik}$, b( k, 1_${ik}$ ), ldb,b( 1_${ik}$, 1_${ik}$ ), ldb )
                           
                 call stdlib${ii}$_${ri}$ger( k-2, nrhs, -one, ap( kc-( k-1 ) ), 1_${ik}$,b( k-1, 1_${ik}$ ), ldb, b( 1_${ik}$, 1_${ik}$ &
                           ), ldb )
                 ! multiply by the inverse of the diagonal block.
                 akm1k = ap( kc+k-2 )
                 akm1 = ap( kc-1 ) / akm1k
                 ak = ap( kc+k-1 ) / akm1k
                 denom = akm1*ak - one
                 do j = 1, nrhs
                    bkm1 = b( k-1, j ) / akm1k
                    bk = b( k, j ) / akm1k
                    b( k-1, j ) = ( ak*bkm1-bk ) / denom
                    b( k, j ) = ( akm1*bk-bkm1 ) / denom
                 end do
                 kc = kc - k + 1_${ik}$
                 k = k - 2_${ik}$
              end if
              go to 10
              30 continue
              ! next solve u**t*x = b, overwriting b with x.
              ! k is the main loop index, increasing from 1 to n in steps of
              ! 1 or 2, depending on the size of the diagonal blocks.
              k = 1_${ik}$
              kc = 1_${ik}$
              40 continue
              ! if k > n, exit from loop.
              if( k>n )go to 50
              if( ipiv( k )>0_${ik}$ ) then
                 ! 1 x 1 diagonal block
                 ! multiply by inv(u**t(k)), where u(k) is the transformation
                 ! stored in column k of a.
                 call stdlib${ii}$_${ri}$gemv( 'TRANSPOSE', k-1, nrhs, -one, b, ldb, ap( kc ),1_${ik}$, one, b( k, &
                           1_${ik}$ ), ldb )
                 ! interchange rows k and ipiv(k).
                 kp = ipiv( k )
                 if( kp/=k )call stdlib${ii}$_${ri}$swap( nrhs, b( k, 1_${ik}$ ), ldb, b( kp, 1_${ik}$ ), ldb )
                 kc = kc + k
                 k = k + 1_${ik}$
              else
                 ! 2 x 2 diagonal block
                 ! multiply by inv(u**t(k+1)), where u(k+1) is the transformation
                 ! stored in columns k and k+1 of a.
                 call stdlib${ii}$_${ri}$gemv( 'TRANSPOSE', k-1, nrhs, -one, b, ldb, ap( kc ),1_${ik}$, one, b( k, &
                           1_${ik}$ ), ldb )
                 call stdlib${ii}$_${ri}$gemv( 'TRANSPOSE', k-1, nrhs, -one, b, ldb,ap( kc+k ), 1_${ik}$, one, b( k+&
                           1_${ik}$, 1_${ik}$ ), ldb )
                 ! interchange rows k and -ipiv(k).
                 kp = -ipiv( k )
                 if( kp/=k )call stdlib${ii}$_${ri}$swap( nrhs, b( k, 1_${ik}$ ), ldb, b( kp, 1_${ik}$ ), ldb )
                 kc = kc + 2_${ik}$*k + 1_${ik}$
                 k = k + 2_${ik}$
              end if
              go to 40
              50 continue
           else
              ! solve a*x = b, where a = l*d*l**t.
              ! first solve l*d*x = b, overwriting b with x.
              ! k is the main loop index, increasing from 1 to n in steps of
              ! 1 or 2, depending on the size of the diagonal blocks.
              k = 1_${ik}$
              kc = 1_${ik}$
              60 continue
              ! if k > n, exit from loop.
              if( k>n )go to 80
              if( ipiv( k )>0_${ik}$ ) then
                 ! 1 x 1 diagonal block
                 ! interchange rows k and ipiv(k).
                 kp = ipiv( k )
                 if( kp/=k )call stdlib${ii}$_${ri}$swap( nrhs, b( k, 1_${ik}$ ), ldb, b( kp, 1_${ik}$ ), ldb )
                 ! multiply by inv(l(k)), where l(k) is the transformation
                 ! stored in column k of a.
                 if( k<n )call stdlib${ii}$_${ri}$ger( n-k, nrhs, -one, ap( kc+1 ), 1_${ik}$, b( k, 1_${ik}$ ),ldb, b( k+1,&
                            1_${ik}$ ), ldb )
                 ! multiply by the inverse of the diagonal block.
                 call stdlib${ii}$_${ri}$scal( nrhs, one / ap( kc ), b( k, 1_${ik}$ ), ldb )
                 kc = kc + n - k + 1_${ik}$
                 k = k + 1_${ik}$
              else
                 ! 2 x 2 diagonal block
                 ! interchange rows k+1 and -ipiv(k).
                 kp = -ipiv( k )
                 if( kp/=k+1 )call stdlib${ii}$_${ri}$swap( nrhs, b( k+1, 1_${ik}$ ), ldb, b( kp, 1_${ik}$ ), ldb )
                 ! multiply by inv(l(k)), where l(k) is the transformation
                 ! stored in columns k and k+1 of a.
                 if( k<n-1 ) then
                    call stdlib${ii}$_${ri}$ger( n-k-1, nrhs, -one, ap( kc+2 ), 1_${ik}$, b( k, 1_${ik}$ ),ldb, b( k+2, 1_${ik}$ )&
                              , ldb )
                    call stdlib${ii}$_${ri}$ger( n-k-1, nrhs, -one, ap( kc+n-k+2 ), 1_${ik}$,b( k+1, 1_${ik}$ ), ldb, b( k+&
                              2_${ik}$, 1_${ik}$ ), ldb )
                 end if
                 ! multiply by the inverse of the diagonal block.
                 akm1k = ap( kc+1 )
                 akm1 = ap( kc ) / akm1k
                 ak = ap( kc+n-k+1 ) / akm1k
                 denom = akm1*ak - one
                 do j = 1, nrhs
                    bkm1 = b( k, j ) / akm1k
                    bk = b( k+1, j ) / akm1k
                    b( k, j ) = ( ak*bkm1-bk ) / denom
                    b( k+1, j ) = ( akm1*bk-bkm1 ) / denom
                 end do
                 kc = kc + 2_${ik}$*( n-k ) + 1_${ik}$
                 k = k + 2_${ik}$
              end if
              go to 60
              80 continue
              ! next solve l**t*x = b, overwriting b with x.
              ! k is the main loop index, decreasing from n to 1 in steps of
              ! 1 or 2, depending on the size of the diagonal blocks.
              k = n
              kc = n*( n+1 ) / 2_${ik}$ + 1_${ik}$
              90 continue
              ! if k < 1, exit from loop.
              if( k<1 )go to 100
              kc = kc - ( n-k+1 )
              if( ipiv( k )>0_${ik}$ ) then
                 ! 1 x 1 diagonal block
                 ! multiply by inv(l**t(k)), where l(k) is the transformation
                 ! stored in column k of a.
                 if( k<n )call stdlib${ii}$_${ri}$gemv( 'TRANSPOSE', n-k, nrhs, -one, b( k+1, 1_${ik}$ ),ldb, ap( &
                           kc+1 ), 1_${ik}$, one, b( k, 1_${ik}$ ), ldb )
                 ! interchange rows k and ipiv(k).
                 kp = ipiv( k )
                 if( kp/=k )call stdlib${ii}$_${ri}$swap( nrhs, b( k, 1_${ik}$ ), ldb, b( kp, 1_${ik}$ ), ldb )
                 k = k - 1_${ik}$
              else
                 ! 2 x 2 diagonal block
                 ! multiply by inv(l**t(k-1)), where l(k-1) is the transformation
                 ! stored in columns k-1 and k of a.
                 if( k<n ) then
                    call stdlib${ii}$_${ri}$gemv( 'TRANSPOSE', n-k, nrhs, -one, b( k+1, 1_${ik}$ ),ldb, ap( kc+1 ), &
                              1_${ik}$, one, b( k, 1_${ik}$ ), ldb )
                    call stdlib${ii}$_${ri}$gemv( 'TRANSPOSE', n-k, nrhs, -one, b( k+1, 1_${ik}$ ),ldb, ap( kc-( n-&
                              k ) ), 1_${ik}$, one, b( k-1, 1_${ik}$ ),ldb )
                 end if
                 ! interchange rows k and -ipiv(k).
                 kp = -ipiv( k )
                 if( kp/=k )call stdlib${ii}$_${ri}$swap( nrhs, b( k, 1_${ik}$ ), ldb, b( kp, 1_${ik}$ ), ldb )
                 kc = kc - ( n-k+2 )
                 k = k - 2_${ik}$
              end if
              go to 90
              100 continue
           end if
           return
     end subroutine stdlib${ii}$_${ri}$sptrs

#:endif
#:endfor

     pure module subroutine stdlib${ii}$_csptrs( uplo, n, nrhs, ap, ipiv, b, ldb, info )
     !! CSPTRS solves a system of linear equations A*X = B with a complex
     !! symmetric matrix A stored in packed format using the factorization
     !! A = U*D*U**T or A = L*D*L**T computed by CSPTRF.
        ! -- 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) :: ldb, n, nrhs
           ! Array Arguments 
           integer(${ik}$), intent(in) :: ipiv(*)
           complex(sp), intent(in) :: ap(*)
           complex(sp), intent(inout) :: b(ldb,*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: upper
           integer(${ik}$) :: j, k, kc, kp
           complex(sp) :: ak, akm1, akm1k, bk, bkm1, denom
           ! Intrinsic Functions 
           ! Executable Statements 
           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( nrhs<0_${ik}$ ) then
              info = -3_${ik}$
           else if( ldb<max( 1_${ik}$, n ) ) then
              info = -7_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'CSPTRS', -info )
              return
           end if
           ! quick return if possible
           if( n==0 .or. nrhs==0 )return
           if( upper ) then
              ! solve a*x = b, where a = u*d*u**t.
              ! first solve u*d*x = b, overwriting b with x.
              ! k is the main loop index, decreasing from n to 1 in steps of
              ! 1 or 2, depending on the size of the diagonal blocks.
              k = n
              kc = n*( n+1 ) / 2_${ik}$ + 1_${ik}$
              10 continue
              ! if k < 1, exit from loop.
              if( k<1 )go to 30
              kc = kc - k
              if( ipiv( k )>0_${ik}$ ) then
                 ! 1 x 1 diagonal block
                 ! interchange rows k and ipiv(k).
                 kp = ipiv( k )
                 if( kp/=k )call stdlib${ii}$_cswap( nrhs, b( k, 1_${ik}$ ), ldb, b( kp, 1_${ik}$ ), ldb )
                 ! multiply by inv(u(k)), where u(k) is the transformation
                 ! stored in column k of a.
                 call stdlib${ii}$_cgeru( k-1, nrhs, -cone, ap( kc ), 1_${ik}$, b( k, 1_${ik}$ ), ldb,b( 1_${ik}$, 1_${ik}$ ), ldb )
                           
                 ! multiply by the inverse of the diagonal block.
                 call stdlib${ii}$_cscal( nrhs, cone / ap( kc+k-1 ), b( k, 1_${ik}$ ), ldb )
                 k = k - 1_${ik}$
              else
                 ! 2 x 2 diagonal block
                 ! interchange rows k-1 and -ipiv(k).
                 kp = -ipiv( k )
                 if( kp/=k-1 )call stdlib${ii}$_cswap( nrhs, b( k-1, 1_${ik}$ ), ldb, b( kp, 1_${ik}$ ), ldb )
                 ! multiply by inv(u(k)), where u(k) is the transformation
                 ! stored in columns k-1 and k of a.
                 call stdlib${ii}$_cgeru( k-2, nrhs, -cone, ap( kc ), 1_${ik}$, b( k, 1_${ik}$ ), ldb,b( 1_${ik}$, 1_${ik}$ ), ldb )
                           
                 call stdlib${ii}$_cgeru( k-2, nrhs, -cone, ap( kc-( k-1 ) ), 1_${ik}$,b( k-1, 1_${ik}$ ), ldb, b( 1_${ik}$, &
                           1_${ik}$ ), ldb )
                 ! multiply by the inverse of the diagonal block.
                 akm1k = ap( kc+k-2 )
                 akm1 = ap( kc-1 ) / akm1k
                 ak = ap( kc+k-1 ) / akm1k
                 denom = akm1*ak - cone
                 do j = 1, nrhs
                    bkm1 = b( k-1, j ) / akm1k
                    bk = b( k, j ) / akm1k
                    b( k-1, j ) = ( ak*bkm1-bk ) / denom
                    b( k, j ) = ( akm1*bk-bkm1 ) / denom
                 end do
                 kc = kc - k + 1_${ik}$
                 k = k - 2_${ik}$
              end if
              go to 10
              30 continue
              ! next solve u**t*x = b, overwriting b with x.
              ! k is the main loop index, increasing from 1 to n in steps of
              ! 1 or 2, depending on the size of the diagonal blocks.
              k = 1_${ik}$
              kc = 1_${ik}$
              40 continue
              ! if k > n, exit from loop.
              if( k>n )go to 50
              if( ipiv( k )>0_${ik}$ ) then
                 ! 1 x 1 diagonal block
                 ! multiply by inv(u**t(k)), where u(k) is the transformation
                 ! stored in column k of a.
                 call stdlib${ii}$_cgemv( 'TRANSPOSE', k-1, nrhs, -cone, b, ldb, ap( kc ),1_${ik}$, cone, b( k,&
                            1_${ik}$ ), ldb )
                 ! interchange rows k and ipiv(k).
                 kp = ipiv( k )
                 if( kp/=k )call stdlib${ii}$_cswap( nrhs, b( k, 1_${ik}$ ), ldb, b( kp, 1_${ik}$ ), ldb )
                 kc = kc + k
                 k = k + 1_${ik}$
              else
                 ! 2 x 2 diagonal block
                 ! multiply by inv(u**t(k+1)), where u(k+1) is the transformation
                 ! stored in columns k and k+1 of a.
                 call stdlib${ii}$_cgemv( 'TRANSPOSE', k-1, nrhs, -cone, b, ldb, ap( kc ),1_${ik}$, cone, b( k,&
                            1_${ik}$ ), ldb )
                 call stdlib${ii}$_cgemv( 'TRANSPOSE', k-1, nrhs, -cone, b, ldb,ap( kc+k ), 1_${ik}$, cone, b( &
                           k+1, 1_${ik}$ ), ldb )
                 ! interchange rows k and -ipiv(k).
                 kp = -ipiv( k )
                 if( kp/=k )call stdlib${ii}$_cswap( nrhs, b( k, 1_${ik}$ ), ldb, b( kp, 1_${ik}$ ), ldb )
                 kc = kc + 2_${ik}$*k + 1_${ik}$
                 k = k + 2_${ik}$
              end if
              go to 40
              50 continue
           else
              ! solve a*x = b, where a = l*d*l**t.
              ! first solve l*d*x = b, overwriting b with x.
              ! k is the main loop index, increasing from 1 to n in steps of
              ! 1 or 2, depending on the size of the diagonal blocks.
              k = 1_${ik}$
              kc = 1_${ik}$
              60 continue
              ! if k > n, exit from loop.
              if( k>n )go to 80
              if( ipiv( k )>0_${ik}$ ) then
                 ! 1 x 1 diagonal block
                 ! interchange rows k and ipiv(k).
                 kp = ipiv( k )
                 if( kp/=k )call stdlib${ii}$_cswap( nrhs, b( k, 1_${ik}$ ), ldb, b( kp, 1_${ik}$ ), ldb )
                 ! multiply by inv(l(k)), where l(k) is the transformation
                 ! stored in column k of a.
                 if( k<n )call stdlib${ii}$_cgeru( n-k, nrhs, -cone, ap( kc+1 ), 1_${ik}$, b( k, 1_${ik}$ ),ldb, b( k+&
                           1_${ik}$, 1_${ik}$ ), ldb )
                 ! multiply by the inverse of the diagonal block.
                 call stdlib${ii}$_cscal( nrhs, cone / ap( kc ), b( k, 1_${ik}$ ), ldb )
                 kc = kc + n - k + 1_${ik}$
                 k = k + 1_${ik}$
              else
                 ! 2 x 2 diagonal block
                 ! interchange rows k+1 and -ipiv(k).
                 kp = -ipiv( k )
                 if( kp/=k+1 )call stdlib${ii}$_cswap( nrhs, b( k+1, 1_${ik}$ ), ldb, b( kp, 1_${ik}$ ), ldb )
                 ! multiply by inv(l(k)), where l(k) is the transformation
                 ! stored in columns k and k+1 of a.
                 if( k<n-1 ) then
                    call stdlib${ii}$_cgeru( n-k-1, nrhs, -cone, ap( kc+2 ), 1_${ik}$, b( k, 1_${ik}$ ),ldb, b( k+2, &
                              1_${ik}$ ), ldb )
                    call stdlib${ii}$_cgeru( n-k-1, nrhs, -cone, ap( kc+n-k+2 ), 1_${ik}$,b( k+1, 1_${ik}$ ), ldb, b( &
                              k+2, 1_${ik}$ ), ldb )
                 end if
                 ! multiply by the inverse of the diagonal block.
                 akm1k = ap( kc+1 )
                 akm1 = ap( kc ) / akm1k
                 ak = ap( kc+n-k+1 ) / akm1k
                 denom = akm1*ak - cone
                 do j = 1, nrhs
                    bkm1 = b( k, j ) / akm1k
                    bk = b( k+1, j ) / akm1k
                    b( k, j ) = ( ak*bkm1-bk ) / denom
                    b( k+1, j ) = ( akm1*bk-bkm1 ) / denom
                 end do
                 kc = kc + 2_${ik}$*( n-k ) + 1_${ik}$
                 k = k + 2_${ik}$
              end if
              go to 60
              80 continue
              ! next solve l**t*x = b, overwriting b with x.
              ! k is the main loop index, decreasing from n to 1 in steps of
              ! 1 or 2, depending on the size of the diagonal blocks.
              k = n
              kc = n*( n+1 ) / 2_${ik}$ + 1_${ik}$
              90 continue
              ! if k < 1, exit from loop.
              if( k<1 )go to 100
              kc = kc - ( n-k+1 )
              if( ipiv( k )>0_${ik}$ ) then
                 ! 1 x 1 diagonal block
                 ! multiply by inv(l**t(k)), where l(k) is the transformation
                 ! stored in column k of a.
                 if( k<n )call stdlib${ii}$_cgemv( 'TRANSPOSE', n-k, nrhs, -cone, b( k+1, 1_${ik}$ ),ldb, ap( &
                           kc+1 ), 1_${ik}$, cone, b( k, 1_${ik}$ ), ldb )
                 ! interchange rows k and ipiv(k).
                 kp = ipiv( k )
                 if( kp/=k )call stdlib${ii}$_cswap( nrhs, b( k, 1_${ik}$ ), ldb, b( kp, 1_${ik}$ ), ldb )
                 k = k - 1_${ik}$
              else
                 ! 2 x 2 diagonal block
                 ! multiply by inv(l**t(k-1)), where l(k-1) is the transformation
                 ! stored in columns k-1 and k of a.
                 if( k<n ) then
                    call stdlib${ii}$_cgemv( 'TRANSPOSE', n-k, nrhs, -cone, b( k+1, 1_${ik}$ ),ldb, ap( kc+1 ),&
                               1_${ik}$, cone, b( k, 1_${ik}$ ), ldb )
                    call stdlib${ii}$_cgemv( 'TRANSPOSE', n-k, nrhs, -cone, b( k+1, 1_${ik}$ ),ldb, ap( kc-( n-&
                              k ) ), 1_${ik}$, cone, b( k-1, 1_${ik}$ ),ldb )
                 end if
                 ! interchange rows k and -ipiv(k).
                 kp = -ipiv( k )
                 if( kp/=k )call stdlib${ii}$_cswap( nrhs, b( k, 1_${ik}$ ), ldb, b( kp, 1_${ik}$ ), ldb )
                 kc = kc - ( n-k+2 )
                 k = k - 2_${ik}$
              end if
              go to 90
              100 continue
           end if
           return
     end subroutine stdlib${ii}$_csptrs

     pure module subroutine stdlib${ii}$_zsptrs( uplo, n, nrhs, ap, ipiv, b, ldb, info )
     !! ZSPTRS solves a system of linear equations A*X = B with a complex
     !! symmetric matrix A stored in packed format using the factorization
     !! A = U*D*U**T or A = L*D*L**T computed by ZSPTRF.
        ! -- 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) :: ldb, n, nrhs
           ! Array Arguments 
           integer(${ik}$), intent(in) :: ipiv(*)
           complex(dp), intent(in) :: ap(*)
           complex(dp), intent(inout) :: b(ldb,*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: upper
           integer(${ik}$) :: j, k, kc, kp
           complex(dp) :: ak, akm1, akm1k, bk, bkm1, denom
           ! Intrinsic Functions 
           ! Executable Statements 
           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( nrhs<0_${ik}$ ) then
              info = -3_${ik}$
           else if( ldb<max( 1_${ik}$, n ) ) then
              info = -7_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'ZSPTRS', -info )
              return
           end if
           ! quick return if possible
           if( n==0 .or. nrhs==0 )return
           if( upper ) then
              ! solve a*x = b, where a = u*d*u**t.
              ! first solve u*d*x = b, overwriting b with x.
              ! k is the main loop index, decreasing from n to 1 in steps of
              ! 1 or 2, depending on the size of the diagonal blocks.
              k = n
              kc = n*( n+1 ) / 2_${ik}$ + 1_${ik}$
              10 continue
              ! if k < 1, exit from loop.
              if( k<1 )go to 30
              kc = kc - k
              if( ipiv( k )>0_${ik}$ ) then
                 ! 1 x 1 diagonal block
                 ! interchange rows k and ipiv(k).
                 kp = ipiv( k )
                 if( kp/=k )call stdlib${ii}$_zswap( nrhs, b( k, 1_${ik}$ ), ldb, b( kp, 1_${ik}$ ), ldb )
                 ! multiply by inv(u(k)), where u(k) is the transformation
                 ! stored in column k of a.
                 call stdlib${ii}$_zgeru( k-1, nrhs, -cone, ap( kc ), 1_${ik}$, b( k, 1_${ik}$ ), ldb,b( 1_${ik}$, 1_${ik}$ ), ldb )
                           
                 ! multiply by the inverse of the diagonal block.
                 call stdlib${ii}$_zscal( nrhs, cone / ap( kc+k-1 ), b( k, 1_${ik}$ ), ldb )
                 k = k - 1_${ik}$
              else
                 ! 2 x 2 diagonal block
                 ! interchange rows k-1 and -ipiv(k).
                 kp = -ipiv( k )
                 if( kp/=k-1 )call stdlib${ii}$_zswap( nrhs, b( k-1, 1_${ik}$ ), ldb, b( kp, 1_${ik}$ ), ldb )
                 ! multiply by inv(u(k)), where u(k) is the transformation
                 ! stored in columns k-1 and k of a.
                 call stdlib${ii}$_zgeru( k-2, nrhs, -cone, ap( kc ), 1_${ik}$, b( k, 1_${ik}$ ), ldb,b( 1_${ik}$, 1_${ik}$ ), ldb )
                           
                 call stdlib${ii}$_zgeru( k-2, nrhs, -cone, ap( kc-( k-1 ) ), 1_${ik}$,b( k-1, 1_${ik}$ ), ldb, b( 1_${ik}$, &
                           1_${ik}$ ), ldb )
                 ! multiply by the inverse of the diagonal block.
                 akm1k = ap( kc+k-2 )
                 akm1 = ap( kc-1 ) / akm1k
                 ak = ap( kc+k-1 ) / akm1k
                 denom = akm1*ak - cone
                 do j = 1, nrhs
                    bkm1 = b( k-1, j ) / akm1k
                    bk = b( k, j ) / akm1k
                    b( k-1, j ) = ( ak*bkm1-bk ) / denom
                    b( k, j ) = ( akm1*bk-bkm1 ) / denom
                 end do
                 kc = kc - k + 1_${ik}$
                 k = k - 2_${ik}$
              end if
              go to 10
              30 continue
              ! next solve u**t*x = b, overwriting b with x.
              ! k is the main loop index, increasing from 1 to n in steps of
              ! 1 or 2, depending on the size of the diagonal blocks.
              k = 1_${ik}$
              kc = 1_${ik}$
              40 continue
              ! if k > n, exit from loop.
              if( k>n )go to 50
              if( ipiv( k )>0_${ik}$ ) then
                 ! 1 x 1 diagonal block
                 ! multiply by inv(u**t(k)), where u(k) is the transformation
                 ! stored in column k of a.
                 call stdlib${ii}$_zgemv( 'TRANSPOSE', k-1, nrhs, -cone, b, ldb, ap( kc ),1_${ik}$, cone, b( k,&
                            1_${ik}$ ), ldb )
                 ! interchange rows k and ipiv(k).
                 kp = ipiv( k )
                 if( kp/=k )call stdlib${ii}$_zswap( nrhs, b( k, 1_${ik}$ ), ldb, b( kp, 1_${ik}$ ), ldb )
                 kc = kc + k
                 k = k + 1_${ik}$
              else
                 ! 2 x 2 diagonal block
                 ! multiply by inv(u**t(k+1)), where u(k+1) is the transformation
                 ! stored in columns k and k+1 of a.
                 call stdlib${ii}$_zgemv( 'TRANSPOSE', k-1, nrhs, -cone, b, ldb, ap( kc ),1_${ik}$, cone, b( k,&
                            1_${ik}$ ), ldb )
                 call stdlib${ii}$_zgemv( 'TRANSPOSE', k-1, nrhs, -cone, b, ldb,ap( kc+k ), 1_${ik}$, cone, b( &
                           k+1, 1_${ik}$ ), ldb )
                 ! interchange rows k and -ipiv(k).
                 kp = -ipiv( k )
                 if( kp/=k )call stdlib${ii}$_zswap( nrhs, b( k, 1_${ik}$ ), ldb, b( kp, 1_${ik}$ ), ldb )
                 kc = kc + 2_${ik}$*k + 1_${ik}$
                 k = k + 2_${ik}$
              end if
              go to 40
              50 continue
           else
              ! solve a*x = b, where a = l*d*l**t.
              ! first solve l*d*x = b, overwriting b with x.
              ! k is the main loop index, increasing from 1 to n in steps of
              ! 1 or 2, depending on the size of the diagonal blocks.
              k = 1_${ik}$
              kc = 1_${ik}$
              60 continue
              ! if k > n, exit from loop.
              if( k>n )go to 80
              if( ipiv( k )>0_${ik}$ ) then
                 ! 1 x 1 diagonal block
                 ! interchange rows k and ipiv(k).
                 kp = ipiv( k )
                 if( kp/=k )call stdlib${ii}$_zswap( nrhs, b( k, 1_${ik}$ ), ldb, b( kp, 1_${ik}$ ), ldb )
                 ! multiply by inv(l(k)), where l(k) is the transformation
                 ! stored in column k of a.
                 if( k<n )call stdlib${ii}$_zgeru( n-k, nrhs, -cone, ap( kc+1 ), 1_${ik}$, b( k, 1_${ik}$ ),ldb, b( k+&
                           1_${ik}$, 1_${ik}$ ), ldb )
                 ! multiply by the inverse of the diagonal block.
                 call stdlib${ii}$_zscal( nrhs, cone / ap( kc ), b( k, 1_${ik}$ ), ldb )
                 kc = kc + n - k + 1_${ik}$
                 k = k + 1_${ik}$
              else
                 ! 2 x 2 diagonal block
                 ! interchange rows k+1 and -ipiv(k).
                 kp = -ipiv( k )
                 if( kp/=k+1 )call stdlib${ii}$_zswap( nrhs, b( k+1, 1_${ik}$ ), ldb, b( kp, 1_${ik}$ ), ldb )
                 ! multiply by inv(l(k)), where l(k) is the transformation
                 ! stored in columns k and k+1 of a.
                 if( k<n-1 ) then
                    call stdlib${ii}$_zgeru( n-k-1, nrhs, -cone, ap( kc+2 ), 1_${ik}$, b( k, 1_${ik}$ ),ldb, b( k+2, &
                              1_${ik}$ ), ldb )
                    call stdlib${ii}$_zgeru( n-k-1, nrhs, -cone, ap( kc+n-k+2 ), 1_${ik}$,b( k+1, 1_${ik}$ ), ldb, b( &
                              k+2, 1_${ik}$ ), ldb )
                 end if
                 ! multiply by the inverse of the diagonal block.
                 akm1k = ap( kc+1 )
                 akm1 = ap( kc ) / akm1k
                 ak = ap( kc+n-k+1 ) / akm1k
                 denom = akm1*ak - cone
                 do j = 1, nrhs
                    bkm1 = b( k, j ) / akm1k
                    bk = b( k+1, j ) / akm1k
                    b( k, j ) = ( ak*bkm1-bk ) / denom
                    b( k+1, j ) = ( akm1*bk-bkm1 ) / denom
                 end do
                 kc = kc + 2_${ik}$*( n-k ) + 1_${ik}$
                 k = k + 2_${ik}$
              end if
              go to 60
              80 continue
              ! next solve l**t*x = b, overwriting b with x.
              ! k is the main loop index, decreasing from n to 1 in steps of
              ! 1 or 2, depending on the size of the diagonal blocks.
              k = n
              kc = n*( n+1 ) / 2_${ik}$ + 1_${ik}$
              90 continue
              ! if k < 1, exit from loop.
              if( k<1 )go to 100
              kc = kc - ( n-k+1 )
              if( ipiv( k )>0_${ik}$ ) then
                 ! 1 x 1 diagonal block
                 ! multiply by inv(l**t(k)), where l(k) is the transformation
                 ! stored in column k of a.
                 if( k<n )call stdlib${ii}$_zgemv( 'TRANSPOSE', n-k, nrhs, -cone, b( k+1, 1_${ik}$ ),ldb, ap( &
                           kc+1 ), 1_${ik}$, cone, b( k, 1_${ik}$ ), ldb )
                 ! interchange rows k and ipiv(k).
                 kp = ipiv( k )
                 if( kp/=k )call stdlib${ii}$_zswap( nrhs, b( k, 1_${ik}$ ), ldb, b( kp, 1_${ik}$ ), ldb )
                 k = k - 1_${ik}$
              else
                 ! 2 x 2 diagonal block
                 ! multiply by inv(l**t(k-1)), where l(k-1) is the transformation
                 ! stored in columns k-1 and k of a.
                 if( k<n ) then
                    call stdlib${ii}$_zgemv( 'TRANSPOSE', n-k, nrhs, -cone, b( k+1, 1_${ik}$ ),ldb, ap( kc+1 ),&
                               1_${ik}$, cone, b( k, 1_${ik}$ ), ldb )
                    call stdlib${ii}$_zgemv( 'TRANSPOSE', n-k, nrhs, -cone, b( k+1, 1_${ik}$ ),ldb, ap( kc-( n-&
                              k ) ), 1_${ik}$, cone, b( k-1, 1_${ik}$ ),ldb )
                 end if
                 ! interchange rows k and -ipiv(k).
                 kp = -ipiv( k )
                 if( kp/=k )call stdlib${ii}$_zswap( nrhs, b( k, 1_${ik}$ ), ldb, b( kp, 1_${ik}$ ), ldb )
                 kc = kc - ( n-k+2 )
                 k = k - 2_${ik}$
              end if
              go to 90
              100 continue
           end if
           return
     end subroutine stdlib${ii}$_zsptrs

#:for ck,ct,ci in CMPLX_KINDS_TYPES
#:if not ck in ["sp","dp"]
     pure module subroutine stdlib${ii}$_${ci}$sptrs( uplo, n, nrhs, ap, ipiv, b, ldb, info )
     !! ZSPTRS: solves a system of linear equations A*X = B with a complex
     !! symmetric matrix A stored in packed format using the factorization
     !! A = U*D*U**T or A = L*D*L**T computed by ZSPTRF.
        ! -- 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) :: ldb, n, nrhs
           ! Array Arguments 
           integer(${ik}$), intent(in) :: ipiv(*)
           complex(${ck}$), intent(in) :: ap(*)
           complex(${ck}$), intent(inout) :: b(ldb,*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: upper
           integer(${ik}$) :: j, k, kc, kp
           complex(${ck}$) :: ak, akm1, akm1k, bk, bkm1, denom
           ! Intrinsic Functions 
           ! Executable Statements 
           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( nrhs<0_${ik}$ ) then
              info = -3_${ik}$
           else if( ldb<max( 1_${ik}$, n ) ) then
              info = -7_${ik}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'ZSPTRS', -info )
              return
           end if
           ! quick return if possible
           if( n==0 .or. nrhs==0 )return
           if( upper ) then
              ! solve a*x = b, where a = u*d*u**t.
              ! first solve u*d*x = b, overwriting b with x.
              ! k is the main loop index, decreasing from n to 1 in steps of
              ! 1 or 2, depending on the size of the diagonal blocks.
              k = n
              kc = n*( n+1 ) / 2_${ik}$ + 1_${ik}$
              10 continue
              ! if k < 1, exit from loop.
              if( k<1 )go to 30
              kc = kc - k
              if( ipiv( k )>0_${ik}$ ) then
                 ! 1 x 1 diagonal block
                 ! interchange rows k and ipiv(k).
                 kp = ipiv( k )
                 if( kp/=k )call stdlib${ii}$_${ci}$swap( nrhs, b( k, 1_${ik}$ ), ldb, b( kp, 1_${ik}$ ), ldb )
                 ! multiply by inv(u(k)), where u(k) is the transformation
                 ! stored in column k of a.
                 call stdlib${ii}$_${ci}$geru( k-1, nrhs, -cone, ap( kc ), 1_${ik}$, b( k, 1_${ik}$ ), ldb,b( 1_${ik}$, 1_${ik}$ ), ldb )
                           
                 ! multiply by the inverse of the diagonal block.
                 call stdlib${ii}$_${ci}$scal( nrhs, cone / ap( kc+k-1 ), b( k, 1_${ik}$ ), ldb )
                 k = k - 1_${ik}$
              else
                 ! 2 x 2 diagonal block
                 ! interchange rows k-1 and -ipiv(k).
                 kp = -ipiv( k )
                 if( kp/=k-1 )call stdlib${ii}$_${ci}$swap( nrhs, b( k-1, 1_${ik}$ ), ldb, b( kp, 1_${ik}$ ), ldb )
                 ! multiply by inv(u(k)), where u(k) is the transformation
                 ! stored in columns k-1 and k of a.
                 call stdlib${ii}$_${ci}$geru( k-2, nrhs, -cone, ap( kc ), 1_${ik}$, b( k, 1_${ik}$ ), ldb,b( 1_${ik}$, 1_${ik}$ ), ldb )
                           
                 call stdlib${ii}$_${ci}$geru( k-2, nrhs, -cone, ap( kc-( k-1 ) ), 1_${ik}$,b( k-1, 1_${ik}$ ), ldb, b( 1_${ik}$, &
                           1_${ik}$ ), ldb )
                 ! multiply by the inverse of the diagonal block.
                 akm1k = ap( kc+k-2 )
                 akm1 = ap( kc-1 ) / akm1k
                 ak = ap( kc+k-1 ) / akm1k
                 denom = akm1*ak - cone
                 do j = 1, nrhs
                    bkm1 = b( k-1, j ) / akm1k
                    bk = b( k, j ) / akm1k
                    b( k-1, j ) = ( ak*bkm1-bk ) / denom
                    b( k, j ) = ( akm1*bk-bkm1 ) / denom
                 end do
                 kc = kc - k + 1_${ik}$
                 k = k - 2_${ik}$
              end if
              go to 10
              30 continue
              ! next solve u**t*x = b, overwriting b with x.
              ! k is the main loop index, increasing from 1 to n in steps of
              ! 1 or 2, depending on the size of the diagonal blocks.
              k = 1_${ik}$
              kc = 1_${ik}$
              40 continue
              ! if k > n, exit from loop.
              if( k>n )go to 50
              if( ipiv( k )>0_${ik}$ ) then
                 ! 1 x 1 diagonal block
                 ! multiply by inv(u**t(k)), where u(k) is the transformation
                 ! stored in column k of a.
                 call stdlib${ii}$_${ci}$gemv( 'TRANSPOSE', k-1, nrhs, -cone, b, ldb, ap( kc ),1_${ik}$, cone, b( k,&
                            1_${ik}$ ), ldb )
                 ! interchange rows k and ipiv(k).
                 kp = ipiv( k )
                 if( kp/=k )call stdlib${ii}$_${ci}$swap( nrhs, b( k, 1_${ik}$ ), ldb, b( kp, 1_${ik}$ ), ldb )
                 kc = kc + k
                 k = k + 1_${ik}$
              else
                 ! 2 x 2 diagonal block
                 ! multiply by inv(u**t(k+1)), where u(k+1) is the transformation
                 ! stored in columns k and k+1 of a.
                 call stdlib${ii}$_${ci}$gemv( 'TRANSPOSE', k-1, nrhs, -cone, b, ldb, ap( kc ),1_${ik}$, cone, b( k,&
                            1_${ik}$ ), ldb )
                 call stdlib${ii}$_${ci}$gemv( 'TRANSPOSE', k-1, nrhs, -cone, b, ldb,ap( kc+k ), 1_${ik}$, cone, b( &
                           k+1, 1_${ik}$ ), ldb )
                 ! interchange rows k and -ipiv(k).
                 kp = -ipiv( k )
                 if( kp/=k )call stdlib${ii}$_${ci}$swap( nrhs, b( k, 1_${ik}$ ), ldb, b( kp, 1_${ik}$ ), ldb )
                 kc = kc + 2_${ik}$*k + 1_${ik}$
                 k = k + 2_${ik}$
              end if
              go to 40
              50 continue
           else
              ! solve a*x = b, where a = l*d*l**t.
              ! first solve l*d*x = b, overwriting b with x.
              ! k is the main loop index, increasing from 1 to n in steps of
              ! 1 or 2, depending on the size of the diagonal blocks.
              k = 1_${ik}$
              kc = 1_${ik}$
              60 continue
              ! if k > n, exit from loop.
              if( k>n )go to 80
              if( ipiv( k )>0_${ik}$ ) then
                 ! 1 x 1 diagonal block
                 ! interchange rows k and ipiv(k).
                 kp = ipiv( k )
                 if( kp/=k )call stdlib${ii}$_${ci}$swap( nrhs, b( k, 1_${ik}$ ), ldb, b( kp, 1_${ik}$ ), ldb )
                 ! multiply by inv(l(k)), where l(k) is the transformation
                 ! stored in column k of a.
                 if( k<n )call stdlib${ii}$_${ci}$geru( n-k, nrhs, -cone, ap( kc+1 ), 1_${ik}$, b( k, 1_${ik}$ ),ldb, b( k+&
                           1_${ik}$, 1_${ik}$ ), ldb )
                 ! multiply by the inverse of the diagonal block.
                 call stdlib${ii}$_${ci}$scal( nrhs, cone / ap( kc ), b( k, 1_${ik}$ ), ldb )
                 kc = kc + n - k + 1_${ik}$
                 k = k + 1_${ik}$
              else
                 ! 2 x 2 diagonal block
                 ! interchange rows k+1 and -ipiv(k).
                 kp = -ipiv( k )
                 if( kp/=k+1 )call stdlib${ii}$_${ci}$swap( nrhs, b( k+1, 1_${ik}$ ), ldb, b( kp, 1_${ik}$ ), ldb )
                 ! multiply by inv(l(k)), where l(k) is the transformation
                 ! stored in columns k and k+1 of a.
                 if( k<n-1 ) then
                    call stdlib${ii}$_${ci}$geru( n-k-1, nrhs, -cone, ap( kc+2 ), 1_${ik}$, b( k, 1_${ik}$ ),ldb, b( k+2, &
                              1_${ik}$ ), ldb )
                    call stdlib${ii}$_${ci}$geru( n-k-1, nrhs, -cone, ap( kc+n-k+2 ), 1_${ik}$,b( k+1, 1_${ik}$ ), ldb, b( &
                              k+2, 1_${ik}$ ), ldb )
                 end if
                 ! multiply by the inverse of the diagonal block.
                 akm1k = ap( kc+1 )
                 akm1 = ap( kc ) / akm1k
                 ak = ap( kc+n-k+1 ) / akm1k
                 denom = akm1*ak - cone
                 do j = 1, nrhs
                    bkm1 = b( k, j ) / akm1k
                    bk = b( k+1, j ) / akm1k
                    b( k, j ) = ( ak*bkm1-bk ) / denom
                    b( k+1, j ) = ( akm1*bk-bkm1 ) / denom
                 end do
                 kc = kc + 2_${ik}$*( n-k ) + 1_${ik}$
                 k = k + 2_${ik}$
              end if
              go to 60
              80 continue
              ! next solve l**t*x = b, overwriting b with x.
              ! k is the main loop index, decreasing from n to 1 in steps of
              ! 1 or 2, depending on the size of the diagonal blocks.
              k = n
              kc = n*( n+1 ) / 2_${ik}$ + 1_${ik}$
              90 continue
              ! if k < 1, exit from loop.
              if( k<1 )go to 100
              kc = kc - ( n-k+1 )
              if( ipiv( k )>0_${ik}$ ) then
                 ! 1 x 1 diagonal block
                 ! multiply by inv(l**t(k)), where l(k) is the transformation
                 ! stored in column k of a.
                 if( k<n )call stdlib${ii}$_${ci}$gemv( 'TRANSPOSE', n-k, nrhs, -cone, b( k+1, 1_${ik}$ ),ldb, ap( &
                           kc+1 ), 1_${ik}$, cone, b( k, 1_${ik}$ ), ldb )
                 ! interchange rows k and ipiv(k).
                 kp = ipiv( k )
                 if( kp/=k )call stdlib${ii}$_${ci}$swap( nrhs, b( k, 1_${ik}$ ), ldb, b( kp, 1_${ik}$ ), ldb )
                 k = k - 1_${ik}$
              else
                 ! 2 x 2 diagonal block
                 ! multiply by inv(l**t(k-1)), where l(k-1) is the transformation
                 ! stored in columns k-1 and k of a.
                 if( k<n ) then
                    call stdlib${ii}$_${ci}$gemv( 'TRANSPOSE', n-k, nrhs, -cone, b( k+1, 1_${ik}$ ),ldb, ap( kc+1 ),&
                               1_${ik}$, cone, b( k, 1_${ik}$ ), ldb )
                    call stdlib${ii}$_${ci}$gemv( 'TRANSPOSE', n-k, nrhs, -cone, b( k+1, 1_${ik}$ ),ldb, ap( kc-( n-&
                              k ) ), 1_${ik}$, cone, b( k-1, 1_${ik}$ ),ldb )
                 end if
                 ! interchange rows k and -ipiv(k).
                 kp = -ipiv( k )
                 if( kp/=k )call stdlib${ii}$_${ci}$swap( nrhs, b( k, 1_${ik}$ ), ldb, b( kp, 1_${ik}$ ), ldb )
                 kc = kc - ( n-k+2 )
                 k = k - 2_${ik}$
              end if
              go to 90
              100 continue
           end if
           return
     end subroutine stdlib${ii}$_${ci}$sptrs

#:endif
#:endfor



     pure module subroutine stdlib${ii}$_ssptri( uplo, n, ap, ipiv, work, info )
     !! SSPTRI computes the inverse of a real symmetric indefinite matrix
     !! A in packed storage using the factorization A = U*D*U**T or
     !! A = L*D*L**T computed by SSPTRF.
        ! -- 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) :: n
           ! Array Arguments 
           integer(${ik}$), intent(in) :: ipiv(*)
           real(sp), intent(inout) :: ap(*)
           real(sp), intent(out) :: work(*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: upper
           integer(${ik}$) :: j, k, kc, kcnext, kp, kpc, kstep, kx, npp
           real(sp) :: ak, akkp1, akp1, d, t, temp
           ! 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}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'SSPTRI', -info )
              return
           end if
           ! quick return if possible
           if( n==0 )return
           ! check that the diagonal matrix d is nonsingular.
           if( upper ) then
              ! upper triangular storage: examine d from bottom to top
              kp = n*( n+1 ) / 2_${ik}$
              do info = n, 1, -1
                 if( ipiv( info )>0 .and. ap( kp )==zero )return
                 kp = kp - info
              end do
           else
              ! lower triangular storage: examine d from top to bottom.
              kp = 1_${ik}$
              do info = 1, n
                 if( ipiv( info )>0 .and. ap( kp )==zero )return
                 kp = kp + n - info + 1_${ik}$
              end do
           end if
           info = 0_${ik}$
           if( upper ) then
              ! compute inv(a) from the factorization a = u*d*u**t.
              ! k is the main loop index, increasing from 1 to n in steps of
              ! 1 or 2, depending on the size of the diagonal blocks.
              k = 1_${ik}$
              kc = 1_${ik}$
              30 continue
              ! if k > n, exit from loop.
              if( k>n )go to 50
              kcnext = kc + k
              if( ipiv( k )>0_${ik}$ ) then
                 ! 1 x 1 diagonal block
                 ! invert the diagonal block.
                 ap( kc+k-1 ) = one / ap( kc+k-1 )
                 ! compute column k of the inverse.
                 if( k>1_${ik}$ ) then
                    call stdlib${ii}$_scopy( k-1, ap( kc ), 1_${ik}$, work, 1_${ik}$ )
                    call stdlib${ii}$_sspmv( uplo, k-1, -one, ap, work, 1_${ik}$, zero, ap( kc ),1_${ik}$ )
                    ap( kc+k-1 ) = ap( kc+k-1 ) -stdlib${ii}$_sdot( k-1, work, 1_${ik}$, ap( kc ), 1_${ik}$ )
                 end if
                 kstep = 1_${ik}$
              else
                 ! 2 x 2 diagonal block
                 ! invert the diagonal block.
                 t = abs( ap( kcnext+k-1 ) )
                 ak = ap( kc+k-1 ) / t
                 akp1 = ap( kcnext+k ) / t
                 akkp1 = ap( kcnext+k-1 ) / t
                 d = t*( ak*akp1-one )
                 ap( kc+k-1 ) = akp1 / d
                 ap( kcnext+k ) = ak / d
                 ap( kcnext+k-1 ) = -akkp1 / d
                 ! compute columns k and k+1 of the inverse.
                 if( k>1_${ik}$ ) then
                    call stdlib${ii}$_scopy( k-1, ap( kc ), 1_${ik}$, work, 1_${ik}$ )
                    call stdlib${ii}$_sspmv( uplo, k-1, -one, ap, work, 1_${ik}$, zero, ap( kc ),1_${ik}$ )
                    ap( kc+k-1 ) = ap( kc+k-1 ) -stdlib${ii}$_sdot( k-1, work, 1_${ik}$, ap( kc ), 1_${ik}$ )
                    ap( kcnext+k-1 ) = ap( kcnext+k-1 ) -stdlib${ii}$_sdot( k-1, ap( kc ), 1_${ik}$, ap( &
                              kcnext ),1_${ik}$ )
                    call stdlib${ii}$_scopy( k-1, ap( kcnext ), 1_${ik}$, work, 1_${ik}$ )
                    call stdlib${ii}$_sspmv( uplo, k-1, -one, ap, work, 1_${ik}$, zero,ap( kcnext ), 1_${ik}$ )
                              
                    ap( kcnext+k ) = ap( kcnext+k ) -stdlib${ii}$_sdot( k-1, work, 1_${ik}$, ap( kcnext ), 1_${ik}$ )
                              
                 end if
                 kstep = 2_${ik}$
                 kcnext = kcnext + k + 1_${ik}$
              end if
              kp = abs( ipiv( k ) )
              if( kp/=k ) then
                 ! interchange rows and columns k and kp in the leading
                 ! submatrix a(1:k+1,1:k+1)
                 kpc = ( kp-1 )*kp / 2_${ik}$ + 1_${ik}$
                 call stdlib${ii}$_sswap( kp-1, ap( kc ), 1_${ik}$, ap( kpc ), 1_${ik}$ )
                 kx = kpc + kp - 1_${ik}$
                 do j = kp + 1, k - 1
                    kx = kx + j - 1_${ik}$
                    temp = ap( kc+j-1 )
                    ap( kc+j-1 ) = ap( kx )
                    ap( kx ) = temp
                 end do
                 temp = ap( kc+k-1 )
                 ap( kc+k-1 ) = ap( kpc+kp-1 )
                 ap( kpc+kp-1 ) = temp
                 if( kstep==2_${ik}$ ) then
                    temp = ap( kc+k+k-1 )
                    ap( kc+k+k-1 ) = ap( kc+k+kp-1 )
                    ap( kc+k+kp-1 ) = temp
                 end if
              end if
              k = k + kstep
              kc = kcnext
              go to 30
              50 continue
           else
              ! compute inv(a) from the factorization a = l*d*l**t.
              ! k is the main loop index, increasing from 1 to n in steps of
              ! 1 or 2, depending on the size of the diagonal blocks.
              npp = n*( n+1 ) / 2_${ik}$
              k = n
              kc = npp
              60 continue
              ! if k < 1, exit from loop.
              if( k<1 )go to 80
              kcnext = kc - ( n-k+2 )
              if( ipiv( k )>0_${ik}$ ) then
                 ! 1 x 1 diagonal block
                 ! invert the diagonal block.
                 ap( kc ) = one / ap( kc )
                 ! compute column k of the inverse.
                 if( k<n ) then
                    call stdlib${ii}$_scopy( n-k, ap( kc+1 ), 1_${ik}$, work, 1_${ik}$ )
                    call stdlib${ii}$_sspmv( uplo, n-k, -one, ap( kc+n-k+1 ), work, 1_${ik}$,zero, ap( kc+1 ), &
                              1_${ik}$ )
                    ap( kc ) = ap( kc ) - stdlib${ii}$_sdot( n-k, work, 1_${ik}$, ap( kc+1 ), 1_${ik}$ )
                 end if
                 kstep = 1_${ik}$
              else
                 ! 2 x 2 diagonal block
                 ! invert the diagonal block.
                 t = abs( ap( kcnext+1 ) )
                 ak = ap( kcnext ) / t
                 akp1 = ap( kc ) / t
                 akkp1 = ap( kcnext+1 ) / t
                 d = t*( ak*akp1-one )
                 ap( kcnext ) = akp1 / d
                 ap( kc ) = ak / d
                 ap( kcnext+1 ) = -akkp1 / d
                 ! compute columns k-1 and k of the inverse.
                 if( k<n ) then
                    call stdlib${ii}$_scopy( n-k, ap( kc+1 ), 1_${ik}$, work, 1_${ik}$ )
                    call stdlib${ii}$_sspmv( uplo, n-k, -one, ap( kc+( n-k+1 ) ), work, 1_${ik}$,zero, ap( kc+&
                              1_${ik}$ ), 1_${ik}$ )
                    ap( kc ) = ap( kc ) - stdlib${ii}$_sdot( n-k, work, 1_${ik}$, ap( kc+1 ), 1_${ik}$ )
                    ap( kcnext+1 ) = ap( kcnext+1 ) -stdlib${ii}$_sdot( n-k, ap( kc+1 ), 1_${ik}$,ap( kcnext+2 &
                              ), 1_${ik}$ )
                    call stdlib${ii}$_scopy( n-k, ap( kcnext+2 ), 1_${ik}$, work, 1_${ik}$ )
                    call stdlib${ii}$_sspmv( uplo, n-k, -one, ap( kc+( n-k+1 ) ), work, 1_${ik}$,zero, ap( &
                              kcnext+2 ), 1_${ik}$ )
                    ap( kcnext ) = ap( kcnext ) -stdlib${ii}$_sdot( n-k, work, 1_${ik}$, ap( kcnext+2 ), 1_${ik}$ )
                              
                 end if
                 kstep = 2_${ik}$
                 kcnext = kcnext - ( n-k+3 )
              end if
              kp = abs( ipiv( k ) )
              if( kp/=k ) then
                 ! interchange rows and columns k and kp in the trailing
                 ! submatrix a(k-1:n,k-1:n)
                 kpc = npp - ( n-kp+1 )*( n-kp+2 ) / 2_${ik}$ + 1_${ik}$
                 if( kp<n )call stdlib${ii}$_sswap( n-kp, ap( kc+kp-k+1 ), 1_${ik}$, ap( kpc+1 ), 1_${ik}$ )
                 kx = kc + kp - k
                 do j = k + 1, kp - 1
                    kx = kx + n - j + 1_${ik}$
                    temp = ap( kc+j-k )
                    ap( kc+j-k ) = ap( kx )
                    ap( kx ) = temp
                 end do
                 temp = ap( kc )
                 ap( kc ) = ap( kpc )
                 ap( kpc ) = temp
                 if( kstep==2_${ik}$ ) then
                    temp = ap( kc-n+k-1 )
                    ap( kc-n+k-1 ) = ap( kc-n+kp-1 )
                    ap( kc-n+kp-1 ) = temp
                 end if
              end if
              k = k - kstep
              kc = kcnext
              go to 60
              80 continue
           end if
           return
     end subroutine stdlib${ii}$_ssptri

     pure module subroutine stdlib${ii}$_dsptri( uplo, n, ap, ipiv, work, info )
     !! DSPTRI computes the inverse of a real symmetric indefinite matrix
     !! A in packed storage using the factorization A = U*D*U**T or
     !! A = L*D*L**T computed by DSPTRF.
        ! -- 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) :: n
           ! Array Arguments 
           integer(${ik}$), intent(in) :: ipiv(*)
           real(dp), intent(inout) :: ap(*)
           real(dp), intent(out) :: work(*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: upper
           integer(${ik}$) :: j, k, kc, kcnext, kp, kpc, kstep, kx, npp
           real(dp) :: ak, akkp1, akp1, d, t, temp
           ! 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}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DSPTRI', -info )
              return
           end if
           ! quick return if possible
           if( n==0 )return
           ! check that the diagonal matrix d is nonsingular.
           if( upper ) then
              ! upper triangular storage: examine d from bottom to top
              kp = n*( n+1 ) / 2_${ik}$
              do info = n, 1, -1
                 if( ipiv( info )>0 .and. ap( kp )==zero )return
                 kp = kp - info
              end do
           else
              ! lower triangular storage: examine d from top to bottom.
              kp = 1_${ik}$
              do info = 1, n
                 if( ipiv( info )>0 .and. ap( kp )==zero )return
                 kp = kp + n - info + 1_${ik}$
              end do
           end if
           info = 0_${ik}$
           if( upper ) then
              ! compute inv(a) from the factorization a = u*d*u**t.
              ! k is the main loop index, increasing from 1 to n in steps of
              ! 1 or 2, depending on the size of the diagonal blocks.
              k = 1_${ik}$
              kc = 1_${ik}$
              30 continue
              ! if k > n, exit from loop.
              if( k>n )go to 50
              kcnext = kc + k
              if( ipiv( k )>0_${ik}$ ) then
                 ! 1 x 1 diagonal block
                 ! invert the diagonal block.
                 ap( kc+k-1 ) = one / ap( kc+k-1 )
                 ! compute column k of the inverse.
                 if( k>1_${ik}$ ) then
                    call stdlib${ii}$_dcopy( k-1, ap( kc ), 1_${ik}$, work, 1_${ik}$ )
                    call stdlib${ii}$_dspmv( uplo, k-1, -one, ap, work, 1_${ik}$, zero, ap( kc ),1_${ik}$ )
                    ap( kc+k-1 ) = ap( kc+k-1 ) -stdlib${ii}$_ddot( k-1, work, 1_${ik}$, ap( kc ), 1_${ik}$ )
                 end if
                 kstep = 1_${ik}$
              else
                 ! 2 x 2 diagonal block
                 ! invert the diagonal block.
                 t = abs( ap( kcnext+k-1 ) )
                 ak = ap( kc+k-1 ) / t
                 akp1 = ap( kcnext+k ) / t
                 akkp1 = ap( kcnext+k-1 ) / t
                 d = t*( ak*akp1-one )
                 ap( kc+k-1 ) = akp1 / d
                 ap( kcnext+k ) = ak / d
                 ap( kcnext+k-1 ) = -akkp1 / d
                 ! compute columns k and k+1 of the inverse.
                 if( k>1_${ik}$ ) then
                    call stdlib${ii}$_dcopy( k-1, ap( kc ), 1_${ik}$, work, 1_${ik}$ )
                    call stdlib${ii}$_dspmv( uplo, k-1, -one, ap, work, 1_${ik}$, zero, ap( kc ),1_${ik}$ )
                    ap( kc+k-1 ) = ap( kc+k-1 ) -stdlib${ii}$_ddot( k-1, work, 1_${ik}$, ap( kc ), 1_${ik}$ )
                    ap( kcnext+k-1 ) = ap( kcnext+k-1 ) -stdlib${ii}$_ddot( k-1, ap( kc ), 1_${ik}$, ap( &
                              kcnext ),1_${ik}$ )
                    call stdlib${ii}$_dcopy( k-1, ap( kcnext ), 1_${ik}$, work, 1_${ik}$ )
                    call stdlib${ii}$_dspmv( uplo, k-1, -one, ap, work, 1_${ik}$, zero,ap( kcnext ), 1_${ik}$ )
                              
                    ap( kcnext+k ) = ap( kcnext+k ) -stdlib${ii}$_ddot( k-1, work, 1_${ik}$, ap( kcnext ), 1_${ik}$ )
                              
                 end if
                 kstep = 2_${ik}$
                 kcnext = kcnext + k + 1_${ik}$
              end if
              kp = abs( ipiv( k ) )
              if( kp/=k ) then
                 ! interchange rows and columns k and kp in the leading
                 ! submatrix a(1:k+1,1:k+1)
                 kpc = ( kp-1 )*kp / 2_${ik}$ + 1_${ik}$
                 call stdlib${ii}$_dswap( kp-1, ap( kc ), 1_${ik}$, ap( kpc ), 1_${ik}$ )
                 kx = kpc + kp - 1_${ik}$
                 do j = kp + 1, k - 1
                    kx = kx + j - 1_${ik}$
                    temp = ap( kc+j-1 )
                    ap( kc+j-1 ) = ap( kx )
                    ap( kx ) = temp
                 end do
                 temp = ap( kc+k-1 )
                 ap( kc+k-1 ) = ap( kpc+kp-1 )
                 ap( kpc+kp-1 ) = temp
                 if( kstep==2_${ik}$ ) then
                    temp = ap( kc+k+k-1 )
                    ap( kc+k+k-1 ) = ap( kc+k+kp-1 )
                    ap( kc+k+kp-1 ) = temp
                 end if
              end if
              k = k + kstep
              kc = kcnext
              go to 30
              50 continue
           else
              ! compute inv(a) from the factorization a = l*d*l**t.
              ! k is the main loop index, increasing from 1 to n in steps of
              ! 1 or 2, depending on the size of the diagonal blocks.
              npp = n*( n+1 ) / 2_${ik}$
              k = n
              kc = npp
              60 continue
              ! if k < 1, exit from loop.
              if( k<1 )go to 80
              kcnext = kc - ( n-k+2 )
              if( ipiv( k )>0_${ik}$ ) then
                 ! 1 x 1 diagonal block
                 ! invert the diagonal block.
                 ap( kc ) = one / ap( kc )
                 ! compute column k of the inverse.
                 if( k<n ) then
                    call stdlib${ii}$_dcopy( n-k, ap( kc+1 ), 1_${ik}$, work, 1_${ik}$ )
                    call stdlib${ii}$_dspmv( uplo, n-k, -one, ap( kc+n-k+1 ), work, 1_${ik}$,zero, ap( kc+1 ), &
                              1_${ik}$ )
                    ap( kc ) = ap( kc ) - stdlib${ii}$_ddot( n-k, work, 1_${ik}$, ap( kc+1 ), 1_${ik}$ )
                 end if
                 kstep = 1_${ik}$
              else
                 ! 2 x 2 diagonal block
                 ! invert the diagonal block.
                 t = abs( ap( kcnext+1 ) )
                 ak = ap( kcnext ) / t
                 akp1 = ap( kc ) / t
                 akkp1 = ap( kcnext+1 ) / t
                 d = t*( ak*akp1-one )
                 ap( kcnext ) = akp1 / d
                 ap( kc ) = ak / d
                 ap( kcnext+1 ) = -akkp1 / d
                 ! compute columns k-1 and k of the inverse.
                 if( k<n ) then
                    call stdlib${ii}$_dcopy( n-k, ap( kc+1 ), 1_${ik}$, work, 1_${ik}$ )
                    call stdlib${ii}$_dspmv( uplo, n-k, -one, ap( kc+( n-k+1 ) ), work, 1_${ik}$,zero, ap( kc+&
                              1_${ik}$ ), 1_${ik}$ )
                    ap( kc ) = ap( kc ) - stdlib${ii}$_ddot( n-k, work, 1_${ik}$, ap( kc+1 ), 1_${ik}$ )
                    ap( kcnext+1 ) = ap( kcnext+1 ) -stdlib${ii}$_ddot( n-k, ap( kc+1 ), 1_${ik}$,ap( kcnext+2 &
                              ), 1_${ik}$ )
                    call stdlib${ii}$_dcopy( n-k, ap( kcnext+2 ), 1_${ik}$, work, 1_${ik}$ )
                    call stdlib${ii}$_dspmv( uplo, n-k, -one, ap( kc+( n-k+1 ) ), work, 1_${ik}$,zero, ap( &
                              kcnext+2 ), 1_${ik}$ )
                    ap( kcnext ) = ap( kcnext ) -stdlib${ii}$_ddot( n-k, work, 1_${ik}$, ap( kcnext+2 ), 1_${ik}$ )
                              
                 end if
                 kstep = 2_${ik}$
                 kcnext = kcnext - ( n-k+3 )
              end if
              kp = abs( ipiv( k ) )
              if( kp/=k ) then
                 ! interchange rows and columns k and kp in the trailing
                 ! submatrix a(k-1:n,k-1:n)
                 kpc = npp - ( n-kp+1 )*( n-kp+2 ) / 2_${ik}$ + 1_${ik}$
                 if( kp<n )call stdlib${ii}$_dswap( n-kp, ap( kc+kp-k+1 ), 1_${ik}$, ap( kpc+1 ), 1_${ik}$ )
                 kx = kc + kp - k
                 do j = k + 1, kp - 1
                    kx = kx + n - j + 1_${ik}$
                    temp = ap( kc+j-k )
                    ap( kc+j-k ) = ap( kx )
                    ap( kx ) = temp
                 end do
                 temp = ap( kc )
                 ap( kc ) = ap( kpc )
                 ap( kpc ) = temp
                 if( kstep==2_${ik}$ ) then
                    temp = ap( kc-n+k-1 )
                    ap( kc-n+k-1 ) = ap( kc-n+kp-1 )
                    ap( kc-n+kp-1 ) = temp
                 end if
              end if
              k = k - kstep
              kc = kcnext
              go to 60
              80 continue
           end if
           return
     end subroutine stdlib${ii}$_dsptri

#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
     pure module subroutine stdlib${ii}$_${ri}$sptri( uplo, n, ap, ipiv, work, info )
     !! DSPTRI: computes the inverse of a real symmetric indefinite matrix
     !! A in packed storage using the factorization A = U*D*U**T or
     !! A = L*D*L**T computed by DSPTRF.
        ! -- 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) :: n
           ! Array Arguments 
           integer(${ik}$), intent(in) :: ipiv(*)
           real(${rk}$), intent(inout) :: ap(*)
           real(${rk}$), intent(out) :: work(*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: upper
           integer(${ik}$) :: j, k, kc, kcnext, kp, kpc, kstep, kx, npp
           real(${rk}$) :: ak, akkp1, akp1, d, t, temp
           ! 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}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'DSPTRI', -info )
              return
           end if
           ! quick return if possible
           if( n==0 )return
           ! check that the diagonal matrix d is nonsingular.
           if( upper ) then
              ! upper triangular storage: examine d from bottom to top
              kp = n*( n+1 ) / 2_${ik}$
              do info = n, 1, -1
                 if( ipiv( info )>0 .and. ap( kp )==zero )return
                 kp = kp - info
              end do
           else
              ! lower triangular storage: examine d from top to bottom.
              kp = 1_${ik}$
              do info = 1, n
                 if( ipiv( info )>0 .and. ap( kp )==zero )return
                 kp = kp + n - info + 1_${ik}$
              end do
           end if
           info = 0_${ik}$
           if( upper ) then
              ! compute inv(a) from the factorization a = u*d*u**t.
              ! k is the main loop index, increasing from 1 to n in steps of
              ! 1 or 2, depending on the size of the diagonal blocks.
              k = 1_${ik}$
              kc = 1_${ik}$
              30 continue
              ! if k > n, exit from loop.
              if( k>n )go to 50
              kcnext = kc + k
              if( ipiv( k )>0_${ik}$ ) then
                 ! 1 x 1 diagonal block
                 ! invert the diagonal block.
                 ap( kc+k-1 ) = one / ap( kc+k-1 )
                 ! compute column k of the inverse.
                 if( k>1_${ik}$ ) then
                    call stdlib${ii}$_${ri}$copy( k-1, ap( kc ), 1_${ik}$, work, 1_${ik}$ )
                    call stdlib${ii}$_${ri}$spmv( uplo, k-1, -one, ap, work, 1_${ik}$, zero, ap( kc ),1_${ik}$ )
                    ap( kc+k-1 ) = ap( kc+k-1 ) -stdlib${ii}$_${ri}$dot( k-1, work, 1_${ik}$, ap( kc ), 1_${ik}$ )
                 end if
                 kstep = 1_${ik}$
              else
                 ! 2 x 2 diagonal block
                 ! invert the diagonal block.
                 t = abs( ap( kcnext+k-1 ) )
                 ak = ap( kc+k-1 ) / t
                 akp1 = ap( kcnext+k ) / t
                 akkp1 = ap( kcnext+k-1 ) / t
                 d = t*( ak*akp1-one )
                 ap( kc+k-1 ) = akp1 / d
                 ap( kcnext+k ) = ak / d
                 ap( kcnext+k-1 ) = -akkp1 / d
                 ! compute columns k and k+1 of the inverse.
                 if( k>1_${ik}$ ) then
                    call stdlib${ii}$_${ri}$copy( k-1, ap( kc ), 1_${ik}$, work, 1_${ik}$ )
                    call stdlib${ii}$_${ri}$spmv( uplo, k-1, -one, ap, work, 1_${ik}$, zero, ap( kc ),1_${ik}$ )
                    ap( kc+k-1 ) = ap( kc+k-1 ) -stdlib${ii}$_${ri}$dot( k-1, work, 1_${ik}$, ap( kc ), 1_${ik}$ )
                    ap( kcnext+k-1 ) = ap( kcnext+k-1 ) -stdlib${ii}$_${ri}$dot( k-1, ap( kc ), 1_${ik}$, ap( &
                              kcnext ),1_${ik}$ )
                    call stdlib${ii}$_${ri}$copy( k-1, ap( kcnext ), 1_${ik}$, work, 1_${ik}$ )
                    call stdlib${ii}$_${ri}$spmv( uplo, k-1, -one, ap, work, 1_${ik}$, zero,ap( kcnext ), 1_${ik}$ )
                              
                    ap( kcnext+k ) = ap( kcnext+k ) -stdlib${ii}$_${ri}$dot( k-1, work, 1_${ik}$, ap( kcnext ), 1_${ik}$ )
                              
                 end if
                 kstep = 2_${ik}$
                 kcnext = kcnext + k + 1_${ik}$
              end if
              kp = abs( ipiv( k ) )
              if( kp/=k ) then
                 ! interchange rows and columns k and kp in the leading
                 ! submatrix a(1:k+1,1:k+1)
                 kpc = ( kp-1 )*kp / 2_${ik}$ + 1_${ik}$
                 call stdlib${ii}$_${ri}$swap( kp-1, ap( kc ), 1_${ik}$, ap( kpc ), 1_${ik}$ )
                 kx = kpc + kp - 1_${ik}$
                 do j = kp + 1, k - 1
                    kx = kx + j - 1_${ik}$
                    temp = ap( kc+j-1 )
                    ap( kc+j-1 ) = ap( kx )
                    ap( kx ) = temp
                 end do
                 temp = ap( kc+k-1 )
                 ap( kc+k-1 ) = ap( kpc+kp-1 )
                 ap( kpc+kp-1 ) = temp
                 if( kstep==2_${ik}$ ) then
                    temp = ap( kc+k+k-1 )
                    ap( kc+k+k-1 ) = ap( kc+k+kp-1 )
                    ap( kc+k+kp-1 ) = temp
                 end if
              end if
              k = k + kstep
              kc = kcnext
              go to 30
              50 continue
           else
              ! compute inv(a) from the factorization a = l*d*l**t.
              ! k is the main loop index, increasing from 1 to n in steps of
              ! 1 or 2, depending on the size of the diagonal blocks.
              npp = n*( n+1 ) / 2_${ik}$
              k = n
              kc = npp
              60 continue
              ! if k < 1, exit from loop.
              if( k<1 )go to 80
              kcnext = kc - ( n-k+2 )
              if( ipiv( k )>0_${ik}$ ) then
                 ! 1 x 1 diagonal block
                 ! invert the diagonal block.
                 ap( kc ) = one / ap( kc )
                 ! compute column k of the inverse.
                 if( k<n ) then
                    call stdlib${ii}$_${ri}$copy( n-k, ap( kc+1 ), 1_${ik}$, work, 1_${ik}$ )
                    call stdlib${ii}$_${ri}$spmv( uplo, n-k, -one, ap( kc+n-k+1 ), work, 1_${ik}$,zero, ap( kc+1 ), &
                              1_${ik}$ )
                    ap( kc ) = ap( kc ) - stdlib${ii}$_${ri}$dot( n-k, work, 1_${ik}$, ap( kc+1 ), 1_${ik}$ )
                 end if
                 kstep = 1_${ik}$
              else
                 ! 2 x 2 diagonal block
                 ! invert the diagonal block.
                 t = abs( ap( kcnext+1 ) )
                 ak = ap( kcnext ) / t
                 akp1 = ap( kc ) / t
                 akkp1 = ap( kcnext+1 ) / t
                 d = t*( ak*akp1-one )
                 ap( kcnext ) = akp1 / d
                 ap( kc ) = ak / d
                 ap( kcnext+1 ) = -akkp1 / d
                 ! compute columns k-1 and k of the inverse.
                 if( k<n ) then
                    call stdlib${ii}$_${ri}$copy( n-k, ap( kc+1 ), 1_${ik}$, work, 1_${ik}$ )
                    call stdlib${ii}$_${ri}$spmv( uplo, n-k, -one, ap( kc+( n-k+1 ) ), work, 1_${ik}$,zero, ap( kc+&
                              1_${ik}$ ), 1_${ik}$ )
                    ap( kc ) = ap( kc ) - stdlib${ii}$_${ri}$dot( n-k, work, 1_${ik}$, ap( kc+1 ), 1_${ik}$ )
                    ap( kcnext+1 ) = ap( kcnext+1 ) -stdlib${ii}$_${ri}$dot( n-k, ap( kc+1 ), 1_${ik}$,ap( kcnext+2 &
                              ), 1_${ik}$ )
                    call stdlib${ii}$_${ri}$copy( n-k, ap( kcnext+2 ), 1_${ik}$, work, 1_${ik}$ )
                    call stdlib${ii}$_${ri}$spmv( uplo, n-k, -one, ap( kc+( n-k+1 ) ), work, 1_${ik}$,zero, ap( &
                              kcnext+2 ), 1_${ik}$ )
                    ap( kcnext ) = ap( kcnext ) -stdlib${ii}$_${ri}$dot( n-k, work, 1_${ik}$, ap( kcnext+2 ), 1_${ik}$ )
                              
                 end if
                 kstep = 2_${ik}$
                 kcnext = kcnext - ( n-k+3 )
              end if
              kp = abs( ipiv( k ) )
              if( kp/=k ) then
                 ! interchange rows and columns k and kp in the trailing
                 ! submatrix a(k-1:n,k-1:n)
                 kpc = npp - ( n-kp+1 )*( n-kp+2 ) / 2_${ik}$ + 1_${ik}$
                 if( kp<n )call stdlib${ii}$_${ri}$swap( n-kp, ap( kc+kp-k+1 ), 1_${ik}$, ap( kpc+1 ), 1_${ik}$ )
                 kx = kc + kp - k
                 do j = k + 1, kp - 1
                    kx = kx + n - j + 1_${ik}$
                    temp = ap( kc+j-k )
                    ap( kc+j-k ) = ap( kx )
                    ap( kx ) = temp
                 end do
                 temp = ap( kc )
                 ap( kc ) = ap( kpc )
                 ap( kpc ) = temp
                 if( kstep==2_${ik}$ ) then
                    temp = ap( kc-n+k-1 )
                    ap( kc-n+k-1 ) = ap( kc-n+kp-1 )
                    ap( kc-n+kp-1 ) = temp
                 end if
              end if
              k = k - kstep
              kc = kcnext
              go to 60
              80 continue
           end if
           return
     end subroutine stdlib${ii}$_${ri}$sptri

#:endif
#:endfor

     pure module subroutine stdlib${ii}$_csptri( uplo, n, ap, ipiv, work, info )
     !! CSPTRI computes the inverse of a complex symmetric indefinite matrix
     !! A in packed storage using the factorization A = U*D*U**T or
     !! A = L*D*L**T computed by CSPTRF.
        ! -- 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) :: n
           ! Array Arguments 
           integer(${ik}$), intent(in) :: ipiv(*)
           complex(sp), intent(inout) :: ap(*)
           complex(sp), intent(out) :: work(*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: upper
           integer(${ik}$) :: j, k, kc, kcnext, kp, kpc, kstep, kx, npp
           complex(sp) :: ak, akkp1, akp1, d, t, temp
           ! 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}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'CSPTRI', -info )
              return
           end if
           ! quick return if possible
           if( n==0 )return
           ! check that the diagonal matrix d is nonsingular.
           if( upper ) then
              ! upper triangular storage: examine d from bottom to top
              kp = n*( n+1 ) / 2_${ik}$
              do info = n, 1, -1
                 if( ipiv( info )>0 .and. ap( kp )==czero )return
                 kp = kp - info
              end do
           else
              ! lower triangular storage: examine d from top to bottom.
              kp = 1_${ik}$
              do info = 1, n
                 if( ipiv( info )>0 .and. ap( kp )==czero )return
                 kp = kp + n - info + 1_${ik}$
              end do
           end if
           info = 0_${ik}$
           if( upper ) then
              ! compute inv(a) from the factorization a = u*d*u**t.
              ! k is the main loop index, increasing from 1 to n in steps of
              ! 1 or 2, depending on the size of the diagonal blocks.
              k = 1_${ik}$
              kc = 1_${ik}$
              30 continue
              ! if k > n, exit from loop.
              if( k>n )go to 50
              kcnext = kc + k
              if( ipiv( k )>0_${ik}$ ) then
                 ! 1 x 1 diagonal block
                 ! invert the diagonal block.
                 ap( kc+k-1 ) = cone / ap( kc+k-1 )
                 ! compute column k of the inverse.
                 if( k>1_${ik}$ ) then
                    call stdlib${ii}$_ccopy( k-1, ap( kc ), 1_${ik}$, work, 1_${ik}$ )
                    call stdlib${ii}$_cspmv( uplo, k-1, -cone, ap, work, 1_${ik}$, czero, ap( kc ),1_${ik}$ )
                    ap( kc+k-1 ) = ap( kc+k-1 ) -stdlib${ii}$_cdotu( k-1, work, 1_${ik}$, ap( kc ), 1_${ik}$ )
                 end if
                 kstep = 1_${ik}$
              else
                 ! 2 x 2 diagonal block
                 ! invert the diagonal block.
                 t = ap( kcnext+k-1 )
                 ak = ap( kc+k-1 ) / t
                 akp1 = ap( kcnext+k ) / t
                 akkp1 = ap( kcnext+k-1 ) / t
                 d = t*( ak*akp1-cone )
                 ap( kc+k-1 ) = akp1 / d
                 ap( kcnext+k ) = ak / d
                 ap( kcnext+k-1 ) = -akkp1 / d
                 ! compute columns k and k+1 of the inverse.
                 if( k>1_${ik}$ ) then
                    call stdlib${ii}$_ccopy( k-1, ap( kc ), 1_${ik}$, work, 1_${ik}$ )
                    call stdlib${ii}$_cspmv( uplo, k-1, -cone, ap, work, 1_${ik}$, czero, ap( kc ),1_${ik}$ )
                    ap( kc+k-1 ) = ap( kc+k-1 ) -stdlib${ii}$_cdotu( k-1, work, 1_${ik}$, ap( kc ), 1_${ik}$ )
                    ap( kcnext+k-1 ) = ap( kcnext+k-1 ) -stdlib${ii}$_cdotu( k-1, ap( kc ), 1_${ik}$, ap( &
                              kcnext ),1_${ik}$ )
                    call stdlib${ii}$_ccopy( k-1, ap( kcnext ), 1_${ik}$, work, 1_${ik}$ )
                    call stdlib${ii}$_cspmv( uplo, k-1, -cone, ap, work, 1_${ik}$, czero,ap( kcnext ), 1_${ik}$ )
                              
                    ap( kcnext+k ) = ap( kcnext+k ) -stdlib${ii}$_cdotu( k-1, work, 1_${ik}$, ap( kcnext ), 1_${ik}$ )
                              
                 end if
                 kstep = 2_${ik}$
                 kcnext = kcnext + k + 1_${ik}$
              end if
              kp = abs( ipiv( k ) )
              if( kp/=k ) then
                 ! interchange rows and columns k and kp in the leading
                 ! submatrix a(1:k+1,1:k+1)
                 kpc = ( kp-1 )*kp / 2_${ik}$ + 1_${ik}$
                 call stdlib${ii}$_cswap( kp-1, ap( kc ), 1_${ik}$, ap( kpc ), 1_${ik}$ )
                 kx = kpc + kp - 1_${ik}$
                 do j = kp + 1, k - 1
                    kx = kx + j - 1_${ik}$
                    temp = ap( kc+j-1 )
                    ap( kc+j-1 ) = ap( kx )
                    ap( kx ) = temp
                 end do
                 temp = ap( kc+k-1 )
                 ap( kc+k-1 ) = ap( kpc+kp-1 )
                 ap( kpc+kp-1 ) = temp
                 if( kstep==2_${ik}$ ) then
                    temp = ap( kc+k+k-1 )
                    ap( kc+k+k-1 ) = ap( kc+k+kp-1 )
                    ap( kc+k+kp-1 ) = temp
                 end if
              end if
              k = k + kstep
              kc = kcnext
              go to 30
              50 continue
           else
              ! compute inv(a) from the factorization a = l*d*l**t.
              ! k is the main loop index, increasing from 1 to n in steps of
              ! 1 or 2, depending on the size of the diagonal blocks.
              npp = n*( n+1 ) / 2_${ik}$
              k = n
              kc = npp
              60 continue
              ! if k < 1, exit from loop.
              if( k<1 )go to 80
              kcnext = kc - ( n-k+2 )
              if( ipiv( k )>0_${ik}$ ) then
                 ! 1 x 1 diagonal block
                 ! invert the diagonal block.
                 ap( kc ) = cone / ap( kc )
                 ! compute column k of the inverse.
                 if( k<n ) then
                    call stdlib${ii}$_ccopy( n-k, ap( kc+1 ), 1_${ik}$, work, 1_${ik}$ )
                    call stdlib${ii}$_cspmv( uplo, n-k, -cone, ap( kc+n-k+1 ), work, 1_${ik}$,czero, ap( kc+1 )&
                              , 1_${ik}$ )
                    ap( kc ) = ap( kc ) - stdlib${ii}$_cdotu( n-k, work, 1_${ik}$, ap( kc+1 ),1_${ik}$ )
                 end if
                 kstep = 1_${ik}$
              else
                 ! 2 x 2 diagonal block
                 ! invert the diagonal block.
                 t = ap( kcnext+1 )
                 ak = ap( kcnext ) / t
                 akp1 = ap( kc ) / t
                 akkp1 = ap( kcnext+1 ) / t
                 d = t*( ak*akp1-cone )
                 ap( kcnext ) = akp1 / d
                 ap( kc ) = ak / d
                 ap( kcnext+1 ) = -akkp1 / d
                 ! compute columns k-1 and k of the inverse.
                 if( k<n ) then
                    call stdlib${ii}$_ccopy( n-k, ap( kc+1 ), 1_${ik}$, work, 1_${ik}$ )
                    call stdlib${ii}$_cspmv( uplo, n-k, -cone, ap( kc+( n-k+1 ) ), work, 1_${ik}$,czero, ap( &
                              kc+1 ), 1_${ik}$ )
                    ap( kc ) = ap( kc ) - stdlib${ii}$_cdotu( n-k, work, 1_${ik}$, ap( kc+1 ),1_${ik}$ )
                    ap( kcnext+1 ) = ap( kcnext+1 ) -stdlib${ii}$_cdotu( n-k, ap( kc+1 ), 1_${ik}$,ap( kcnext+&
                              2_${ik}$ ), 1_${ik}$ )
                    call stdlib${ii}$_ccopy( n-k, ap( kcnext+2 ), 1_${ik}$, work, 1_${ik}$ )
                    call stdlib${ii}$_cspmv( uplo, n-k, -cone, ap( kc+( n-k+1 ) ), work, 1_${ik}$,czero, ap( &
                              kcnext+2 ), 1_${ik}$ )
                    ap( kcnext ) = ap( kcnext ) -stdlib${ii}$_cdotu( n-k, work, 1_${ik}$, ap( kcnext+2 ), 1_${ik}$ )
                              
                 end if
                 kstep = 2_${ik}$
                 kcnext = kcnext - ( n-k+3 )
              end if
              kp = abs( ipiv( k ) )
              if( kp/=k ) then
                 ! interchange rows and columns k and kp in the trailing
                 ! submatrix a(k-1:n,k-1:n)
                 kpc = npp - ( n-kp+1 )*( n-kp+2 ) / 2_${ik}$ + 1_${ik}$
                 if( kp<n )call stdlib${ii}$_cswap( n-kp, ap( kc+kp-k+1 ), 1_${ik}$, ap( kpc+1 ), 1_${ik}$ )
                 kx = kc + kp - k
                 do j = k + 1, kp - 1
                    kx = kx + n - j + 1_${ik}$
                    temp = ap( kc+j-k )
                    ap( kc+j-k ) = ap( kx )
                    ap( kx ) = temp
                 end do
                 temp = ap( kc )
                 ap( kc ) = ap( kpc )
                 ap( kpc ) = temp
                 if( kstep==2_${ik}$ ) then
                    temp = ap( kc-n+k-1 )
                    ap( kc-n+k-1 ) = ap( kc-n+kp-1 )
                    ap( kc-n+kp-1 ) = temp
                 end if
              end if
              k = k - kstep
              kc = kcnext
              go to 60
              80 continue
           end if
           return
     end subroutine stdlib${ii}$_csptri

     pure module subroutine stdlib${ii}$_zsptri( uplo, n, ap, ipiv, work, info )
     !! ZSPTRI computes the inverse of a complex symmetric indefinite matrix
     !! A in packed storage using the factorization A = U*D*U**T or
     !! A = L*D*L**T computed by ZSPTRF.
        ! -- 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) :: n
           ! Array Arguments 
           integer(${ik}$), intent(in) :: ipiv(*)
           complex(dp), intent(inout) :: ap(*)
           complex(dp), intent(out) :: work(*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: upper
           integer(${ik}$) :: j, k, kc, kcnext, kp, kpc, kstep, kx, npp
           complex(dp) :: ak, akkp1, akp1, d, t, temp
           ! 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}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'ZSPTRI', -info )
              return
           end if
           ! quick return if possible
           if( n==0 )return
           ! check that the diagonal matrix d is nonsingular.
           if( upper ) then
              ! upper triangular storage: examine d from bottom to top
              kp = n*( n+1 ) / 2_${ik}$
              do info = n, 1, -1
                 if( ipiv( info )>0 .and. ap( kp )==czero )return
                 kp = kp - info
              end do
           else
              ! lower triangular storage: examine d from top to bottom.
              kp = 1_${ik}$
              do info = 1, n
                 if( ipiv( info )>0 .and. ap( kp )==czero )return
                 kp = kp + n - info + 1_${ik}$
              end do
           end if
           info = 0_${ik}$
           if( upper ) then
              ! compute inv(a) from the factorization a = u*d*u**t.
              ! k is the main loop index, increasing from 1 to n in steps of
              ! 1 or 2, depending on the size of the diagonal blocks.
              k = 1_${ik}$
              kc = 1_${ik}$
              30 continue
              ! if k > n, exit from loop.
              if( k>n )go to 50
              kcnext = kc + k
              if( ipiv( k )>0_${ik}$ ) then
                 ! 1 x 1 diagonal block
                 ! invert the diagonal block.
                 ap( kc+k-1 ) = cone / ap( kc+k-1 )
                 ! compute column k of the inverse.
                 if( k>1_${ik}$ ) then
                    call stdlib${ii}$_zcopy( k-1, ap( kc ), 1_${ik}$, work, 1_${ik}$ )
                    call stdlib${ii}$_zspmv( uplo, k-1, -cone, ap, work, 1_${ik}$, czero, ap( kc ),1_${ik}$ )
                    ap( kc+k-1 ) = ap( kc+k-1 ) -stdlib${ii}$_zdotu( k-1, work, 1_${ik}$, ap( kc ), 1_${ik}$ )
                 end if
                 kstep = 1_${ik}$
              else
                 ! 2 x 2 diagonal block
                 ! invert the diagonal block.
                 t = ap( kcnext+k-1 )
                 ak = ap( kc+k-1 ) / t
                 akp1 = ap( kcnext+k ) / t
                 akkp1 = ap( kcnext+k-1 ) / t
                 d = t*( ak*akp1-cone )
                 ap( kc+k-1 ) = akp1 / d
                 ap( kcnext+k ) = ak / d
                 ap( kcnext+k-1 ) = -akkp1 / d
                 ! compute columns k and k+1 of the inverse.
                 if( k>1_${ik}$ ) then
                    call stdlib${ii}$_zcopy( k-1, ap( kc ), 1_${ik}$, work, 1_${ik}$ )
                    call stdlib${ii}$_zspmv( uplo, k-1, -cone, ap, work, 1_${ik}$, czero, ap( kc ),1_${ik}$ )
                    ap( kc+k-1 ) = ap( kc+k-1 ) -stdlib${ii}$_zdotu( k-1, work, 1_${ik}$, ap( kc ), 1_${ik}$ )
                    ap( kcnext+k-1 ) = ap( kcnext+k-1 ) -stdlib${ii}$_zdotu( k-1, ap( kc ), 1_${ik}$, ap( &
                              kcnext ),1_${ik}$ )
                    call stdlib${ii}$_zcopy( k-1, ap( kcnext ), 1_${ik}$, work, 1_${ik}$ )
                    call stdlib${ii}$_zspmv( uplo, k-1, -cone, ap, work, 1_${ik}$, czero,ap( kcnext ), 1_${ik}$ )
                              
                    ap( kcnext+k ) = ap( kcnext+k ) -stdlib${ii}$_zdotu( k-1, work, 1_${ik}$, ap( kcnext ), 1_${ik}$ )
                              
                 end if
                 kstep = 2_${ik}$
                 kcnext = kcnext + k + 1_${ik}$
              end if
              kp = abs( ipiv( k ) )
              if( kp/=k ) then
                 ! interchange rows and columns k and kp in the leading
                 ! submatrix a(1:k+1,1:k+1)
                 kpc = ( kp-1 )*kp / 2_${ik}$ + 1_${ik}$
                 call stdlib${ii}$_zswap( kp-1, ap( kc ), 1_${ik}$, ap( kpc ), 1_${ik}$ )
                 kx = kpc + kp - 1_${ik}$
                 do j = kp + 1, k - 1
                    kx = kx + j - 1_${ik}$
                    temp = ap( kc+j-1 )
                    ap( kc+j-1 ) = ap( kx )
                    ap( kx ) = temp
                 end do
                 temp = ap( kc+k-1 )
                 ap( kc+k-1 ) = ap( kpc+kp-1 )
                 ap( kpc+kp-1 ) = temp
                 if( kstep==2_${ik}$ ) then
                    temp = ap( kc+k+k-1 )
                    ap( kc+k+k-1 ) = ap( kc+k+kp-1 )
                    ap( kc+k+kp-1 ) = temp
                 end if
              end if
              k = k + kstep
              kc = kcnext
              go to 30
              50 continue
           else
              ! compute inv(a) from the factorization a = l*d*l**t.
              ! k is the main loop index, increasing from 1 to n in steps of
              ! 1 or 2, depending on the size of the diagonal blocks.
              npp = n*( n+1 ) / 2_${ik}$
              k = n
              kc = npp
              60 continue
              ! if k < 1, exit from loop.
              if( k<1 )go to 80
              kcnext = kc - ( n-k+2 )
              if( ipiv( k )>0_${ik}$ ) then
                 ! 1 x 1 diagonal block
                 ! invert the diagonal block.
                 ap( kc ) = cone / ap( kc )
                 ! compute column k of the inverse.
                 if( k<n ) then
                    call stdlib${ii}$_zcopy( n-k, ap( kc+1 ), 1_${ik}$, work, 1_${ik}$ )
                    call stdlib${ii}$_zspmv( uplo, n-k, -cone, ap( kc+n-k+1 ), work, 1_${ik}$,czero, ap( kc+1 )&
                              , 1_${ik}$ )
                    ap( kc ) = ap( kc ) - stdlib${ii}$_zdotu( n-k, work, 1_${ik}$, ap( kc+1 ),1_${ik}$ )
                 end if
                 kstep = 1_${ik}$
              else
                 ! 2 x 2 diagonal block
                 ! invert the diagonal block.
                 t = ap( kcnext+1 )
                 ak = ap( kcnext ) / t
                 akp1 = ap( kc ) / t
                 akkp1 = ap( kcnext+1 ) / t
                 d = t*( ak*akp1-cone )
                 ap( kcnext ) = akp1 / d
                 ap( kc ) = ak / d
                 ap( kcnext+1 ) = -akkp1 / d
                 ! compute columns k-1 and k of the inverse.
                 if( k<n ) then
                    call stdlib${ii}$_zcopy( n-k, ap( kc+1 ), 1_${ik}$, work, 1_${ik}$ )
                    call stdlib${ii}$_zspmv( uplo, n-k, -cone, ap( kc+( n-k+1 ) ), work, 1_${ik}$,czero, ap( &
                              kc+1 ), 1_${ik}$ )
                    ap( kc ) = ap( kc ) - stdlib${ii}$_zdotu( n-k, work, 1_${ik}$, ap( kc+1 ),1_${ik}$ )
                    ap( kcnext+1 ) = ap( kcnext+1 ) -stdlib${ii}$_zdotu( n-k, ap( kc+1 ), 1_${ik}$,ap( kcnext+&
                              2_${ik}$ ), 1_${ik}$ )
                    call stdlib${ii}$_zcopy( n-k, ap( kcnext+2 ), 1_${ik}$, work, 1_${ik}$ )
                    call stdlib${ii}$_zspmv( uplo, n-k, -cone, ap( kc+( n-k+1 ) ), work, 1_${ik}$,czero, ap( &
                              kcnext+2 ), 1_${ik}$ )
                    ap( kcnext ) = ap( kcnext ) -stdlib${ii}$_zdotu( n-k, work, 1_${ik}$, ap( kcnext+2 ), 1_${ik}$ )
                              
                 end if
                 kstep = 2_${ik}$
                 kcnext = kcnext - ( n-k+3 )
              end if
              kp = abs( ipiv( k ) )
              if( kp/=k ) then
                 ! interchange rows and columns k and kp in the trailing
                 ! submatrix a(k-1:n,k-1:n)
                 kpc = npp - ( n-kp+1 )*( n-kp+2 ) / 2_${ik}$ + 1_${ik}$
                 if( kp<n )call stdlib${ii}$_zswap( n-kp, ap( kc+kp-k+1 ), 1_${ik}$, ap( kpc+1 ), 1_${ik}$ )
                 kx = kc + kp - k
                 do j = k + 1, kp - 1
                    kx = kx + n - j + 1_${ik}$
                    temp = ap( kc+j-k )
                    ap( kc+j-k ) = ap( kx )
                    ap( kx ) = temp
                 end do
                 temp = ap( kc )
                 ap( kc ) = ap( kpc )
                 ap( kpc ) = temp
                 if( kstep==2_${ik}$ ) then
                    temp = ap( kc-n+k-1 )
                    ap( kc-n+k-1 ) = ap( kc-n+kp-1 )
                    ap( kc-n+kp-1 ) = temp
                 end if
              end if
              k = k - kstep
              kc = kcnext
              go to 60
              80 continue
           end if
           return
     end subroutine stdlib${ii}$_zsptri

#:for ck,ct,ci in CMPLX_KINDS_TYPES
#:if not ck in ["sp","dp"]
     pure module subroutine stdlib${ii}$_${ci}$sptri( uplo, n, ap, ipiv, work, info )
     !! ZSPTRI: computes the inverse of a complex symmetric indefinite matrix
     !! A in packed storage using the factorization A = U*D*U**T or
     !! A = L*D*L**T computed by ZSPTRF.
        ! -- 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) :: n
           ! Array Arguments 
           integer(${ik}$), intent(in) :: ipiv(*)
           complex(${ck}$), intent(inout) :: ap(*)
           complex(${ck}$), intent(out) :: work(*)
        ! =====================================================================
           
           ! Local Scalars 
           logical(lk) :: upper
           integer(${ik}$) :: j, k, kc, kcnext, kp, kpc, kstep, kx, npp
           complex(${ck}$) :: ak, akkp1, akp1, d, t, temp
           ! 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}$
           end if
           if( info/=0_${ik}$ ) then
              call stdlib${ii}$_xerbla( 'ZSPTRI', -info )
              return
           end if
           ! quick return if possible
           if( n==0 )return
           ! check that the diagonal matrix d is nonsingular.
           if( upper ) then
              ! upper triangular storage: examine d from bottom to top
              kp = n*( n+1 ) / 2_${ik}$
              do info = n, 1, -1
                 if( ipiv( info )>0 .and. ap( kp )==czero )return
                 kp = kp - info
              end do
           else
              ! lower triangular storage: examine d from top to bottom.
              kp = 1_${ik}$
              do info = 1, n
                 if( ipiv( info )>0 .and. ap( kp )==czero )return
                 kp = kp + n - info + 1_${ik}$
              end do
           end if
           info = 0_${ik}$
           if( upper ) then
              ! compute inv(a) from the factorization a = u*d*u**t.
              ! k is the main loop index, increasing from 1 to n in steps of
              ! 1 or 2, depending on the size of the diagonal blocks.
              k = 1_${ik}$
              kc = 1_${ik}$
              30 continue
              ! if k > n, exit from loop.
              if( k>n )go to 50
              kcnext = kc + k
              if( ipiv( k )>0_${ik}$ ) then
                 ! 1 x 1 diagonal block
                 ! invert the diagonal block.
                 ap( kc+k-1 ) = cone / ap( kc+k-1 )
                 ! compute column k of the inverse.
                 if( k>1_${ik}$ ) then
                    call stdlib${ii}$_${ci}$copy( k-1, ap( kc ), 1_${ik}$, work, 1_${ik}$ )
                    call stdlib${ii}$_${ci}$spmv( uplo, k-1, -cone, ap, work, 1_${ik}$, czero, ap( kc ),1_${ik}$ )
                    ap( kc+k-1 ) = ap( kc+k-1 ) -stdlib${ii}$_${ci}$dotu( k-1, work, 1_${ik}$, ap( kc ), 1_${ik}$ )
                 end if
                 kstep = 1_${ik}$
              else
                 ! 2 x 2 diagonal block
                 ! invert the diagonal block.
                 t = ap( kcnext+k-1 )
                 ak = ap( kc+k-1 ) / t
                 akp1 = ap( kcnext+k ) / t
                 akkp1 = ap( kcnext+k-1 ) / t
                 d = t*( ak*akp1-cone )
                 ap( kc+k-1 ) = akp1 / d
                 ap( kcnext+k ) = ak / d
                 ap( kcnext+k-1 ) = -akkp1 / d
                 ! compute columns k and k+1 of the inverse.
                 if( k>1_${ik}$ ) then
                    call stdlib${ii}$_${ci}$copy( k-1, ap( kc ), 1_${ik}$, work, 1_${ik}$ )
                    call stdlib${ii}$_${ci}$spmv( uplo, k-1, -cone, ap, work, 1_${ik}$, czero, ap( kc ),1_${ik}$ )
                    ap( kc+k-1 ) = ap( kc+k-1 ) -stdlib${ii}$_${ci}$dotu( k-1, work, 1_${ik}$, ap( kc ), 1_${ik}$ )
                    ap( kcnext+k-1 ) = ap( kcnext+k-1 ) -stdlib${ii}$_${ci}$dotu( k-1, ap( kc ), 1_${ik}$, ap( &
                              kcnext ),1_${ik}$ )