#: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}$ )