#:include "common.fypp"
submodule(stdlib_lapack_solve) stdlib_lapack_solve_chol
implicit none
contains
#:for ik,it,ii in LINALG_INT_KINDS_TYPES
pure module subroutine stdlib${ii}$_sposv( uplo, n, nrhs, a, lda, b, ldb, info )
!! SPOSV computes the solution to a real system of linear equations
!! A * X = B,
!! where A is an N-by-N symmetric positive definite matrix and X and B
!! are N-by-NRHS matrices.
!! The Cholesky decomposition is used to factor A as
!! A = U**T* U, if UPLO = 'U', or
!! A = L * L**T, if UPLO = 'L',
!! where U is an upper triangular matrix and L is a lower triangular
!! matrix. The factored form of A is then used to solve the system of
!! equations A * X = B.
! -- lapack driver routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: uplo
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: lda, ldb, n, nrhs
! Array Arguments
real(sp), intent(inout) :: a(lda,*), b(ldb,*)
! =====================================================================
! Intrinsic Functions
! Executable Statements
! test the input parameters.
info = 0_${ik}$
if( .not.stdlib_lsame( uplo, 'U' ) .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( lda<max( 1_${ik}$, n ) ) then
info = -5_${ik}$
else if( ldb<max( 1_${ik}$, n ) ) then
info = -7_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'SPOSV ', -info )
return
end if
! compute the cholesky factorization a = u**t*u or a = l*l**t.
call stdlib${ii}$_spotrf( uplo, n, a, lda, info )
if( info==0_${ik}$ ) then
! solve the system a*x = b, overwriting b with x.
call stdlib${ii}$_spotrs( uplo, n, nrhs, a, lda, b, ldb, info )
end if
return
end subroutine stdlib${ii}$_sposv
pure module subroutine stdlib${ii}$_dposv( uplo, n, nrhs, a, lda, b, ldb, info )
!! DPOSV computes the solution to a real system of linear equations
!! A * X = B,
!! where A is an N-by-N symmetric positive definite matrix and X and B
!! are N-by-NRHS matrices.
!! The Cholesky decomposition is used to factor A as
!! A = U**T* U, if UPLO = 'U', or
!! A = L * L**T, if UPLO = 'L',
!! where U is an upper triangular matrix and L is a lower triangular
!! matrix. The factored form of A is then used to solve the system of
!! equations A * X = B.
! -- lapack driver routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: uplo
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: lda, ldb, n, nrhs
! Array Arguments
real(dp), intent(inout) :: a(lda,*), b(ldb,*)
! =====================================================================
! Intrinsic Functions
! Executable Statements
! test the input parameters.
info = 0_${ik}$
if( .not.stdlib_lsame( uplo, 'U' ) .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( lda<max( 1_${ik}$, n ) ) then
info = -5_${ik}$
else if( ldb<max( 1_${ik}$, n ) ) then
info = -7_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'DPOSV ', -info )
return
end if
! compute the cholesky factorization a = u**t*u or a = l*l**t.
call stdlib${ii}$_dpotrf( uplo, n, a, lda, info )
if( info==0_${ik}$ ) then
! solve the system a*x = b, overwriting b with x.
call stdlib${ii}$_dpotrs( uplo, n, nrhs, a, lda, b, ldb, info )
end if
return
end subroutine stdlib${ii}$_dposv
#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
pure module subroutine stdlib${ii}$_${ri}$posv( uplo, n, nrhs, a, lda, b, ldb, info )
!! DPOSV: computes the solution to a real system of linear equations
!! A * X = B,
!! where A is an N-by-N symmetric positive definite matrix and X and B
!! are N-by-NRHS matrices.
!! The Cholesky decomposition is used to factor A as
!! A = U**T* U, if UPLO = 'U', or
!! A = L * L**T, if UPLO = 'L',
!! where U is an upper triangular matrix and L is a lower triangular
!! matrix. The factored form of A is then used to solve the system of
!! equations A * X = B.
! -- lapack driver routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_${rk}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: uplo
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: lda, ldb, n, nrhs
! Array Arguments
real(${rk}$), intent(inout) :: a(lda,*), b(ldb,*)
! =====================================================================
! Intrinsic Functions
! Executable Statements
! test the input parameters.
info = 0_${ik}$
if( .not.stdlib_lsame( uplo, 'U' ) .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( lda<max( 1_${ik}$, n ) ) then
info = -5_${ik}$
else if( ldb<max( 1_${ik}$, n ) ) then
info = -7_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'DPOSV ', -info )
return
end if
! compute the cholesky factorization a = u**t*u or a = l*l**t.
call stdlib${ii}$_${ri}$potrf( uplo, n, a, lda, info )
if( info==0_${ik}$ ) then
! solve the system a*x = b, overwriting b with x.
call stdlib${ii}$_${ri}$potrs( uplo, n, nrhs, a, lda, b, ldb, info )
end if
return
end subroutine stdlib${ii}$_${ri}$posv
#:endif
#:endfor
pure module subroutine stdlib${ii}$_cposv( uplo, n, nrhs, a, lda, b, ldb, info )
!! CPOSV computes the solution to a complex system of linear equations
!! A * X = B,
!! where A is an N-by-N Hermitian positive definite matrix and X and B
!! are N-by-NRHS matrices.
!! The Cholesky decomposition is used to factor A as
!! A = U**H* U, if UPLO = 'U', or
!! A = L * L**H, if UPLO = 'L',
!! where U is an upper triangular matrix and L is a lower triangular
!! matrix. The factored form of A is then used to solve the system of
!! equations A * X = B.
! -- lapack driver routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: uplo
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: lda, ldb, n, nrhs
! Array Arguments
complex(sp), intent(inout) :: a(lda,*), b(ldb,*)
! =====================================================================
! Intrinsic Functions
! Executable Statements
! test the input parameters.
info = 0_${ik}$
if( .not.stdlib_lsame( uplo, 'U' ) .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( lda<max( 1_${ik}$, n ) ) then
info = -5_${ik}$
else if( ldb<max( 1_${ik}$, n ) ) then
info = -7_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'CPOSV ', -info )
return
end if
! compute the cholesky factorization a = u**h*u or a = l*l**h.
call stdlib${ii}$_cpotrf( uplo, n, a, lda, info )
if( info==0_${ik}$ ) then
! solve the system a*x = b, overwriting b with x.
call stdlib${ii}$_cpotrs( uplo, n, nrhs, a, lda, b, ldb, info )
end if
return
end subroutine stdlib${ii}$_cposv
pure module subroutine stdlib${ii}$_zposv( uplo, n, nrhs, a, lda, b, ldb, info )
!! ZPOSV computes the solution to a complex system of linear equations
!! A * X = B,
!! where A is an N-by-N Hermitian positive definite matrix and X and B
!! are N-by-NRHS matrices.
!! The Cholesky decomposition is used to factor A as
!! A = U**H* U, if UPLO = 'U', or
!! A = L * L**H, if UPLO = 'L',
!! where U is an upper triangular matrix and L is a lower triangular
!! matrix. The factored form of A is then used to solve the system of
!! equations A * X = B.
! -- lapack driver routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: uplo
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: lda, ldb, n, nrhs
! Array Arguments
complex(dp), intent(inout) :: a(lda,*), b(ldb,*)
! =====================================================================
! Intrinsic Functions
! Executable Statements
! test the input parameters.
info = 0_${ik}$
if( .not.stdlib_lsame( uplo, 'U' ) .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( lda<max( 1_${ik}$, n ) ) then
info = -5_${ik}$
else if( ldb<max( 1_${ik}$, n ) ) then
info = -7_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'ZPOSV ', -info )
return
end if
! compute the cholesky factorization a = u**h *u or a = l*l**h.
call stdlib${ii}$_zpotrf( uplo, n, a, lda, info )
if( info==0_${ik}$ ) then
! solve the system a*x = b, overwriting b with x.
call stdlib${ii}$_zpotrs( uplo, n, nrhs, a, lda, b, ldb, info )
end if
return
end subroutine stdlib${ii}$_zposv
#:for ck,ct,ci in CMPLX_KINDS_TYPES
#:if not ck in ["sp","dp"]
pure module subroutine stdlib${ii}$_${ci}$posv( uplo, n, nrhs, a, lda, b, ldb, info )
!! ZPOSV: computes the solution to a complex system of linear equations
!! A * X = B,
!! where A is an N-by-N Hermitian positive definite matrix and X and B
!! are N-by-NRHS matrices.
!! The Cholesky decomposition is used to factor A as
!! A = U**H* U, if UPLO = 'U', or
!! A = L * L**H, if UPLO = 'L',
!! where U is an upper triangular matrix and L is a lower triangular
!! matrix. The factored form of A is then used to solve the system of
!! equations A * X = B.
! -- lapack driver routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_${ck}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: uplo
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: lda, ldb, n, nrhs
! Array Arguments
complex(${ck}$), intent(inout) :: a(lda,*), b(ldb,*)
! =====================================================================
! Intrinsic Functions
! Executable Statements
! test the input parameters.
info = 0_${ik}$
if( .not.stdlib_lsame( uplo, 'U' ) .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( lda<max( 1_${ik}$, n ) ) then
info = -5_${ik}$
else if( ldb<max( 1_${ik}$, n ) ) then
info = -7_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'ZPOSV ', -info )
return
end if
! compute the cholesky factorization a = u**h *u or a = l*l**h.
call stdlib${ii}$_${ci}$potrf( uplo, n, a, lda, info )
if( info==0_${ik}$ ) then
! solve the system a*x = b, overwriting b with x.
call stdlib${ii}$_${ci}$potrs( uplo, n, nrhs, a, lda, b, ldb, info )
end if
return
end subroutine stdlib${ii}$_${ci}$posv
#:endif
#:endfor
module subroutine stdlib${ii}$_sposvx( fact, uplo, n, nrhs, a, lda, af, ldaf, equed,s, b, ldb, x, ldx, &
!! SPOSVX uses the Cholesky factorization A = U**T*U or A = L*L**T to
!! compute the solution to a real system of linear equations
!! A * X = B,
!! where A is an N-by-N symmetric positive definite matrix and X and B
!! are N-by-NRHS matrices.
!! Error bounds on the solution and a condition estimate are also
!! provided.
rcond, ferr, berr, work,iwork, info )
! -- lapack driver 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(inout) :: equed
character, intent(in) :: fact, uplo
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: lda, ldaf, ldb, ldx, n, nrhs
real(sp), intent(out) :: rcond
! Array Arguments
integer(${ik}$), intent(out) :: iwork(*)
real(sp), intent(inout) :: a(lda,*), af(ldaf,*), b(ldb,*), s(*)
real(sp), intent(out) :: berr(*), ferr(*), work(*), x(ldx,*)
! =====================================================================
! Local Scalars
logical(lk) :: equil, nofact, rcequ
integer(${ik}$) :: i, infequ, j
real(sp) :: amax, anorm, bignum, scond, smax, smin, smlnum
! Intrinsic Functions
! Executable Statements
info = 0_${ik}$
nofact = stdlib_lsame( fact, 'N' )
equil = stdlib_lsame( fact, 'E' )
if( nofact .or. equil ) then
equed = 'N'
rcequ = .false.
else
rcequ = stdlib_lsame( equed, 'Y' )
smlnum = stdlib${ii}$_slamch( 'SAFE MINIMUM' )
bignum = one / smlnum
end if
! test the input parameters.
if( .not.nofact .and. .not.equil .and. .not.stdlib_lsame( fact, 'F' ) )then
info = -1_${ik}$
else if( .not.stdlib_lsame( uplo, 'U' ) .and. .not.stdlib_lsame( uplo, 'L' ) )&
then
info = -2_${ik}$
else if( n<0_${ik}$ ) then
info = -3_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -4_${ik}$
else if( lda<max( 1_${ik}$, n ) ) then
info = -6_${ik}$
else if( ldaf<max( 1_${ik}$, n ) ) then
info = -8_${ik}$
else if( stdlib_lsame( fact, 'F' ) .and. .not.( rcequ .or. stdlib_lsame( equed, 'N' ) )&
) then
info = -9_${ik}$
else
if( rcequ ) then
smin = bignum
smax = zero
do j = 1, n
smin = min( smin, s( j ) )
smax = max( smax, s( j ) )
end do
if( smin<=zero ) then
info = -10_${ik}$
else if( n>0_${ik}$ ) then
scond = max( smin, smlnum ) / min( smax, bignum )
else
scond = one
end if
end if
if( info==0_${ik}$ ) then
if( ldb<max( 1_${ik}$, n ) ) then
info = -12_${ik}$
else if( ldx<max( 1_${ik}$, n ) ) then
info = -14_${ik}$
end if
end if
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'SPOSVX', -info )
return
end if
if( equil ) then
! compute row and column scalings to equilibrate the matrix a.
call stdlib${ii}$_spoequ( n, a, lda, s, scond, amax, infequ )
if( infequ==0_${ik}$ ) then
! equilibrate the matrix.
call stdlib${ii}$_slaqsy( uplo, n, a, lda, s, scond, amax, equed )
rcequ = stdlib_lsame( equed, 'Y' )
end if
end if
! scale the right hand side.
if( rcequ ) then
do j = 1, nrhs
do i = 1, n
b( i, j ) = s( i )*b( i, j )
end do
end do
end if
if( nofact .or. equil ) then
! compute the cholesky factorization a = u**t *u or a = l*l**t.
call stdlib${ii}$_slacpy( uplo, n, n, a, lda, af, ldaf )
call stdlib${ii}$_spotrf( uplo, n, af, ldaf, info )
! return if info is non-zero.
if( info>0_${ik}$ )then
rcond = zero
return
end if
end if
! compute the norm of the matrix a.
anorm = stdlib${ii}$_slansy( '1', uplo, n, a, lda, work )
! compute the reciprocal of the condition number of a.
call stdlib${ii}$_spocon( uplo, n, af, ldaf, anorm, rcond, work, iwork, info )
! compute the solution matrix x.
call stdlib${ii}$_slacpy( 'FULL', n, nrhs, b, ldb, x, ldx )
call stdlib${ii}$_spotrs( uplo, n, nrhs, af, ldaf, x, ldx, info )
! use iterative refinement to improve the computed solution and
! compute error bounds and backward error estimates for it.
call stdlib${ii}$_sporfs( uplo, n, nrhs, a, lda, af, ldaf, b, ldb, x, ldx,ferr, berr, work, &
iwork, info )
! transform the solution matrix x to a solution of the original
! system.
if( rcequ ) then
do j = 1, nrhs
do i = 1, n
x( i, j ) = s( i )*x( i, j )
end do
end do
do j = 1, nrhs
ferr( j ) = ferr( j ) / scond
end do
end if
! set info = n+1 if the matrix is singular to working precision.
if( rcond<stdlib${ii}$_slamch( 'EPSILON' ) )info = n + 1_${ik}$
return
end subroutine stdlib${ii}$_sposvx
module subroutine stdlib${ii}$_dposvx( fact, uplo, n, nrhs, a, lda, af, ldaf, equed,s, b, ldb, x, ldx, &
!! DPOSVX uses the Cholesky factorization A = U**T*U or A = L*L**T to
!! compute the solution to a real system of linear equations
!! A * X = B,
!! where A is an N-by-N symmetric positive definite matrix and X and B
!! are N-by-NRHS matrices.
!! Error bounds on the solution and a condition estimate are also
!! provided.
rcond, ferr, berr, work,iwork, info )
! -- lapack driver 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(inout) :: equed
character, intent(in) :: fact, uplo
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: lda, ldaf, ldb, ldx, n, nrhs
real(dp), intent(out) :: rcond
! Array Arguments
integer(${ik}$), intent(out) :: iwork(*)
real(dp), intent(inout) :: a(lda,*), af(ldaf,*), b(ldb,*), s(*)
real(dp), intent(out) :: berr(*), ferr(*), work(*), x(ldx,*)
! =====================================================================
! Local Scalars
logical(lk) :: equil, nofact, rcequ
integer(${ik}$) :: i, infequ, j
real(dp) :: amax, anorm, bignum, scond, smax, smin, smlnum
! Intrinsic Functions
! Executable Statements
info = 0_${ik}$
nofact = stdlib_lsame( fact, 'N' )
equil = stdlib_lsame( fact, 'E' )
if( nofact .or. equil ) then
equed = 'N'
rcequ = .false.
else
rcequ = stdlib_lsame( equed, 'Y' )
smlnum = stdlib${ii}$_dlamch( 'SAFE MINIMUM' )
bignum = one / smlnum
end if
! test the input parameters.
if( .not.nofact .and. .not.equil .and. .not.stdlib_lsame( fact, 'F' ) )then
info = -1_${ik}$
else if( .not.stdlib_lsame( uplo, 'U' ) .and. .not.stdlib_lsame( uplo, 'L' ) )&
then
info = -2_${ik}$
else if( n<0_${ik}$ ) then
info = -3_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -4_${ik}$
else if( lda<max( 1_${ik}$, n ) ) then
info = -6_${ik}$
else if( ldaf<max( 1_${ik}$, n ) ) then
info = -8_${ik}$
else if( stdlib_lsame( fact, 'F' ) .and. .not.( rcequ .or. stdlib_lsame( equed, 'N' ) )&
) then
info = -9_${ik}$
else
if( rcequ ) then
smin = bignum
smax = zero
do j = 1, n
smin = min( smin, s( j ) )
smax = max( smax, s( j ) )
end do
if( smin<=zero ) then
info = -10_${ik}$
else if( n>0_${ik}$ ) then
scond = max( smin, smlnum ) / min( smax, bignum )
else
scond = one
end if
end if
if( info==0_${ik}$ ) then
if( ldb<max( 1_${ik}$, n ) ) then
info = -12_${ik}$
else if( ldx<max( 1_${ik}$, n ) ) then
info = -14_${ik}$
end if
end if
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'DPOSVX', -info )
return
end if
if( equil ) then
! compute row and column scalings to equilibrate the matrix a.
call stdlib${ii}$_dpoequ( n, a, lda, s, scond, amax, infequ )
if( infequ==0_${ik}$ ) then
! equilibrate the matrix.
call stdlib${ii}$_dlaqsy( uplo, n, a, lda, s, scond, amax, equed )
rcequ = stdlib_lsame( equed, 'Y' )
end if
end if
! scale the right hand side.
if( rcequ ) then
do j = 1, nrhs
do i = 1, n
b( i, j ) = s( i )*b( i, j )
end do
end do
end if
if( nofact .or. equil ) then
! compute the cholesky factorization a = u**t *u or a = l*l**t.
call stdlib${ii}$_dlacpy( uplo, n, n, a, lda, af, ldaf )
call stdlib${ii}$_dpotrf( uplo, n, af, ldaf, info )
! return if info is non-zero.
if( info>0_${ik}$ )then
rcond = zero
return
end if
end if
! compute the norm of the matrix a.
anorm = stdlib${ii}$_dlansy( '1', uplo, n, a, lda, work )
! compute the reciprocal of the condition number of a.
call stdlib${ii}$_dpocon( uplo, n, af, ldaf, anorm, rcond, work, iwork, info )
! compute the solution matrix x.
call stdlib${ii}$_dlacpy( 'FULL', n, nrhs, b, ldb, x, ldx )
call stdlib${ii}$_dpotrs( uplo, n, nrhs, af, ldaf, x, ldx, info )
! use iterative refinement to improve the computed solution and
! compute error bounds and backward error estimates for it.
call stdlib${ii}$_dporfs( uplo, n, nrhs, a, lda, af, ldaf, b, ldb, x, ldx,ferr, berr, work, &
iwork, info )
! transform the solution matrix x to a solution of the original
! system.
if( rcequ ) then
do j = 1, nrhs
do i = 1, n
x( i, j ) = s( i )*x( i, j )
end do
end do
do j = 1, nrhs
ferr( j ) = ferr( j ) / scond
end do
end if
! set info = n+1 if the matrix is singular to working precision.
if( rcond<stdlib${ii}$_dlamch( 'EPSILON' ) )info = n + 1_${ik}$
return
end subroutine stdlib${ii}$_dposvx
#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
module subroutine stdlib${ii}$_${ri}$posvx( fact, uplo, n, nrhs, a, lda, af, ldaf, equed,s, b, ldb, x, ldx, &
!! DPOSVX: uses the Cholesky factorization A = U**T*U or A = L*L**T to
!! compute the solution to a real system of linear equations
!! A * X = B,
!! where A is an N-by-N symmetric positive definite matrix and X and B
!! are N-by-NRHS matrices.
!! Error bounds on the solution and a condition estimate are also
!! provided.
rcond, ferr, berr, work,iwork, info )
! -- lapack driver 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(inout) :: equed
character, intent(in) :: fact, uplo
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: lda, ldaf, ldb, ldx, n, nrhs
real(${rk}$), intent(out) :: rcond
! Array Arguments
integer(${ik}$), intent(out) :: iwork(*)
real(${rk}$), intent(inout) :: a(lda,*), af(ldaf,*), b(ldb,*), s(*)
real(${rk}$), intent(out) :: berr(*), ferr(*), work(*), x(ldx,*)
! =====================================================================
! Local Scalars
logical(lk) :: equil, nofact, rcequ
integer(${ik}$) :: i, infequ, j
real(${rk}$) :: amax, anorm, bignum, scond, smax, smin, smlnum
! Intrinsic Functions
! Executable Statements
info = 0_${ik}$
nofact = stdlib_lsame( fact, 'N' )
equil = stdlib_lsame( fact, 'E' )
if( nofact .or. equil ) then
equed = 'N'
rcequ = .false.
else
rcequ = stdlib_lsame( equed, 'Y' )
smlnum = stdlib${ii}$_${ri}$lamch( 'SAFE MINIMUM' )
bignum = one / smlnum
end if
! test the input parameters.
if( .not.nofact .and. .not.equil .and. .not.stdlib_lsame( fact, 'F' ) )then
info = -1_${ik}$
else if( .not.stdlib_lsame( uplo, 'U' ) .and. .not.stdlib_lsame( uplo, 'L' ) )&
then
info = -2_${ik}$
else if( n<0_${ik}$ ) then
info = -3_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -4_${ik}$
else if( lda<max( 1_${ik}$, n ) ) then
info = -6_${ik}$
else if( ldaf<max( 1_${ik}$, n ) ) then
info = -8_${ik}$
else if( stdlib_lsame( fact, 'F' ) .and. .not.( rcequ .or. stdlib_lsame( equed, 'N' ) )&
) then
info = -9_${ik}$
else
if( rcequ ) then
smin = bignum
smax = zero
do j = 1, n
smin = min( smin, s( j ) )
smax = max( smax, s( j ) )
end do
if( smin<=zero ) then
info = -10_${ik}$
else if( n>0_${ik}$ ) then
scond = max( smin, smlnum ) / min( smax, bignum )
else
scond = one
end if
end if
if( info==0_${ik}$ ) then
if( ldb<max( 1_${ik}$, n ) ) then
info = -12_${ik}$
else if( ldx<max( 1_${ik}$, n ) ) then
info = -14_${ik}$
end if
end if
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'DPOSVX', -info )
return
end if
if( equil ) then
! compute row and column scalings to equilibrate the matrix a.
call stdlib${ii}$_${ri}$poequ( n, a, lda, s, scond, amax, infequ )
if( infequ==0_${ik}$ ) then
! equilibrate the matrix.
call stdlib${ii}$_${ri}$laqsy( uplo, n, a, lda, s, scond, amax, equed )
rcequ = stdlib_lsame( equed, 'Y' )
end if
end if
! scale the right hand side.
if( rcequ ) then
do j = 1, nrhs
do i = 1, n
b( i, j ) = s( i )*b( i, j )
end do
end do
end if
if( nofact .or. equil ) then
! compute the cholesky factorization a = u**t *u or a = l*l**t.
call stdlib${ii}$_${ri}$lacpy( uplo, n, n, a, lda, af, ldaf )
call stdlib${ii}$_${ri}$potrf( uplo, n, af, ldaf, info )
! return if info is non-zero.
if( info>0_${ik}$ )then
rcond = zero
return
end if
end if
! compute the norm of the matrix a.
anorm = stdlib${ii}$_${ri}$lansy( '1', uplo, n, a, lda, work )
! compute the reciprocal of the condition number of a.
call stdlib${ii}$_${ri}$pocon( uplo, n, af, ldaf, anorm, rcond, work, iwork, info )
! compute the solution matrix x.
call stdlib${ii}$_${ri}$lacpy( 'FULL', n, nrhs, b, ldb, x, ldx )
call stdlib${ii}$_${ri}$potrs( uplo, n, nrhs, af, ldaf, x, ldx, info )
! use iterative refinement to improve the computed solution and
! compute error bounds and backward error estimates for it.
call stdlib${ii}$_${ri}$porfs( uplo, n, nrhs, a, lda, af, ldaf, b, ldb, x, ldx,ferr, berr, work, &
iwork, info )
! transform the solution matrix x to a solution of the original
! system.
if( rcequ ) then
do j = 1, nrhs
do i = 1, n
x( i, j ) = s( i )*x( i, j )
end do
end do
do j = 1, nrhs
ferr( j ) = ferr( j ) / scond
end do
end if
! set info = n+1 if the matrix is singular to working precision.
if( rcond<stdlib${ii}$_${ri}$lamch( 'EPSILON' ) )info = n + 1_${ik}$
return
end subroutine stdlib${ii}$_${ri}$posvx
#:endif
#:endfor
module subroutine stdlib${ii}$_cposvx( fact, uplo, n, nrhs, a, lda, af, ldaf, equed,s, b, ldb, x, ldx, &
!! CPOSVX uses the Cholesky factorization A = U**H*U or A = L*L**H to
!! compute the solution to a complex system of linear equations
!! A * X = B,
!! where A is an N-by-N Hermitian positive definite matrix and X and B
!! are N-by-NRHS matrices.
!! Error bounds on the solution and a condition estimate are also
!! provided.
rcond, ferr, berr, work,rwork, info )
! -- lapack driver 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(inout) :: equed
character, intent(in) :: fact, uplo
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: lda, ldaf, ldb, ldx, n, nrhs
real(sp), intent(out) :: rcond
! Array Arguments
real(sp), intent(out) :: berr(*), ferr(*), rwork(*)
real(sp), intent(inout) :: s(*)
complex(sp), intent(inout) :: a(lda,*), af(ldaf,*), b(ldb,*)
complex(sp), intent(out) :: work(*), x(ldx,*)
! =====================================================================
! Local Scalars
logical(lk) :: equil, nofact, rcequ
integer(${ik}$) :: i, infequ, j
real(sp) :: amax, anorm, bignum, scond, smax, smin, smlnum
! Intrinsic Functions
! Executable Statements
info = 0_${ik}$
nofact = stdlib_lsame( fact, 'N' )
equil = stdlib_lsame( fact, 'E' )
if( nofact .or. equil ) then
equed = 'N'
rcequ = .false.
else
rcequ = stdlib_lsame( equed, 'Y' )
smlnum = stdlib${ii}$_slamch( 'SAFE MINIMUM' )
bignum = one / smlnum
end if
! test the input parameters.
if( .not.nofact .and. .not.equil .and. .not.stdlib_lsame( fact, 'F' ) )then
info = -1_${ik}$
else if( .not.stdlib_lsame( uplo, 'U' ) .and. .not.stdlib_lsame( uplo, 'L' ) )&
then
info = -2_${ik}$
else if( n<0_${ik}$ ) then
info = -3_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -4_${ik}$
else if( lda<max( 1_${ik}$, n ) ) then
info = -6_${ik}$
else if( ldaf<max( 1_${ik}$, n ) ) then
info = -8_${ik}$
else if( stdlib_lsame( fact, 'F' ) .and. .not.( rcequ .or. stdlib_lsame( equed, 'N' ) )&
) then
info = -9_${ik}$
else
if( rcequ ) then
smin = bignum
smax = zero
do j = 1, n
smin = min( smin, s( j ) )
smax = max( smax, s( j ) )
end do
if( smin<=zero ) then
info = -10_${ik}$
else if( n>0_${ik}$ ) then
scond = max( smin, smlnum ) / min( smax, bignum )
else
scond = one
end if
end if
if( info==0_${ik}$ ) then
if( ldb<max( 1_${ik}$, n ) ) then
info = -12_${ik}$
else if( ldx<max( 1_${ik}$, n ) ) then
info = -14_${ik}$
end if
end if
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'CPOSVX', -info )
return
end if
if( equil ) then
! compute row and column scalings to equilibrate the matrix a.
call stdlib${ii}$_cpoequ( n, a, lda, s, scond, amax, infequ )
if( infequ==0_${ik}$ ) then
! equilibrate the matrix.
call stdlib${ii}$_claqhe( uplo, n, a, lda, s, scond, amax, equed )
rcequ = stdlib_lsame( equed, 'Y' )
end if
end if
! scale the right hand side.
if( rcequ ) then
do j = 1, nrhs
do i = 1, n
b( i, j ) = s( i )*b( i, j )
end do
end do
end if
if( nofact .or. equil ) then
! compute the cholesky factorization a = u**h *u or a = l*l**h.
call stdlib${ii}$_clacpy( uplo, n, n, a, lda, af, ldaf )
call stdlib${ii}$_cpotrf( uplo, n, af, ldaf, info )
! return if info is non-zero.
if( info>0_${ik}$ )then
rcond = zero
return
end if
end if
! compute the norm of the matrix a.
anorm = stdlib${ii}$_clanhe( '1', uplo, n, a, lda, rwork )
! compute the reciprocal of the condition number of a.
call stdlib${ii}$_cpocon( uplo, n, af, ldaf, anorm, rcond, work, rwork, info )
! compute the solution matrix x.
call stdlib${ii}$_clacpy( 'FULL', n, nrhs, b, ldb, x, ldx )
call stdlib${ii}$_cpotrs( uplo, n, nrhs, af, ldaf, x, ldx, info )
! use iterative refinement to improve the computed solution and
! compute error bounds and backward error estimates for it.
call stdlib${ii}$_cporfs( uplo, n, nrhs, a, lda, af, ldaf, b, ldb, x, ldx,ferr, berr, work, &
rwork, info )
! transform the solution matrix x to a solution of the original
! system.
if( rcequ ) then
do j = 1, nrhs
do i = 1, n
x( i, j ) = s( i )*x( i, j )
end do
end do
do j = 1, nrhs
ferr( j ) = ferr( j ) / scond
end do
end if
! set info = n+1 if the matrix is singular to working precision.
if( rcond<stdlib${ii}$_slamch( 'EPSILON' ) )info = n + 1_${ik}$
return
end subroutine stdlib${ii}$_cposvx
module subroutine stdlib${ii}$_zposvx( fact, uplo, n, nrhs, a, lda, af, ldaf, equed,s, b, ldb, x, ldx, &
!! ZPOSVX uses the Cholesky factorization A = U**H*U or A = L*L**H to
!! compute the solution to a complex system of linear equations
!! A * X = B,
!! where A is an N-by-N Hermitian positive definite matrix and X and B
!! are N-by-NRHS matrices.
!! Error bounds on the solution and a condition estimate are also
!! provided.
rcond, ferr, berr, work,rwork, info )
! -- lapack driver 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(inout) :: equed
character, intent(in) :: fact, uplo
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: lda, ldaf, ldb, ldx, n, nrhs
real(dp), intent(out) :: rcond
! Array Arguments
real(dp), intent(out) :: berr(*), ferr(*), rwork(*)
real(dp), intent(inout) :: s(*)
complex(dp), intent(inout) :: a(lda,*), af(ldaf,*), b(ldb,*)
complex(dp), intent(out) :: work(*), x(ldx,*)
! =====================================================================
! Local Scalars
logical(lk) :: equil, nofact, rcequ
integer(${ik}$) :: i, infequ, j
real(dp) :: amax, anorm, bignum, scond, smax, smin, smlnum
! Intrinsic Functions
! Executable Statements
info = 0_${ik}$
nofact = stdlib_lsame( fact, 'N' )
equil = stdlib_lsame( fact, 'E' )
if( nofact .or. equil ) then
equed = 'N'
rcequ = .false.
else
rcequ = stdlib_lsame( equed, 'Y' )
smlnum = stdlib${ii}$_dlamch( 'SAFE MINIMUM' )
bignum = one / smlnum
end if
! test the input parameters.
if( .not.nofact .and. .not.equil .and. .not.stdlib_lsame( fact, 'F' ) )then
info = -1_${ik}$
else if( .not.stdlib_lsame( uplo, 'U' ) .and. .not.stdlib_lsame( uplo, 'L' ) )&
then
info = -2_${ik}$
else if( n<0_${ik}$ ) then
info = -3_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -4_${ik}$
else if( lda<max( 1_${ik}$, n ) ) then
info = -6_${ik}$
else if( ldaf<max( 1_${ik}$, n ) ) then
info = -8_${ik}$
else if( stdlib_lsame( fact, 'F' ) .and. .not.( rcequ .or. stdlib_lsame( equed, 'N' ) )&
) then
info = -9_${ik}$
else
if( rcequ ) then
smin = bignum
smax = zero
do j = 1, n
smin = min( smin, s( j ) )
smax = max( smax, s( j ) )
end do
if( smin<=zero ) then
info = -10_${ik}$
else if( n>0_${ik}$ ) then
scond = max( smin, smlnum ) / min( smax, bignum )
else
scond = one
end if
end if
if( info==0_${ik}$ ) then
if( ldb<max( 1_${ik}$, n ) ) then
info = -12_${ik}$
else if( ldx<max( 1_${ik}$, n ) ) then
info = -14_${ik}$
end if
end if
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'ZPOSVX', -info )
return
end if
if( equil ) then
! compute row and column scalings to equilibrate the matrix a.
call stdlib${ii}$_zpoequ( n, a, lda, s, scond, amax, infequ )
if( infequ==0_${ik}$ ) then
! equilibrate the matrix.
call stdlib${ii}$_zlaqhe( uplo, n, a, lda, s, scond, amax, equed )
rcequ = stdlib_lsame( equed, 'Y' )
end if
end if
! scale the right hand side.
if( rcequ ) then
do j = 1, nrhs
do i = 1, n
b( i, j ) = s( i )*b( i, j )
end do
end do
end if
if( nofact .or. equil ) then
! compute the cholesky factorization a = u**h *u or a = l*l**h.
call stdlib${ii}$_zlacpy( uplo, n, n, a, lda, af, ldaf )
call stdlib${ii}$_zpotrf( uplo, n, af, ldaf, info )
! return if info is non-zero.
if( info>0_${ik}$ )then
rcond = zero
return
end if
end if
! compute the norm of the matrix a.
anorm = stdlib${ii}$_zlanhe( '1', uplo, n, a, lda, rwork )
! compute the reciprocal of the condition number of a.
call stdlib${ii}$_zpocon( uplo, n, af, ldaf, anorm, rcond, work, rwork, info )
! compute the solution matrix x.
call stdlib${ii}$_zlacpy( 'FULL', n, nrhs, b, ldb, x, ldx )
call stdlib${ii}$_zpotrs( uplo, n, nrhs, af, ldaf, x, ldx, info )
! use iterative refinement to improve the computed solution and
! compute error bounds and backward error estimates for it.
call stdlib${ii}$_zporfs( uplo, n, nrhs, a, lda, af, ldaf, b, ldb, x, ldx,ferr, berr, work, &
rwork, info )
! transform the solution matrix x to a solution of the original
! system.
if( rcequ ) then
do j = 1, nrhs
do i = 1, n
x( i, j ) = s( i )*x( i, j )
end do
end do
do j = 1, nrhs
ferr( j ) = ferr( j ) / scond
end do
end if
! set info = n+1 if the matrix is singular to working precision.
if( rcond<stdlib${ii}$_dlamch( 'EPSILON' ) )info = n + 1_${ik}$
return
end subroutine stdlib${ii}$_zposvx
#:for ck,ct,ci in CMPLX_KINDS_TYPES
#:if not ck in ["sp","dp"]
module subroutine stdlib${ii}$_${ci}$posvx( fact, uplo, n, nrhs, a, lda, af, ldaf, equed,s, b, ldb, x, ldx, &
!! ZPOSVX: uses the Cholesky factorization A = U**H*U or A = L*L**H to
!! compute the solution to a complex system of linear equations
!! A * X = B,
!! where A is an N-by-N Hermitian positive definite matrix and X and B
!! are N-by-NRHS matrices.
!! Error bounds on the solution and a condition estimate are also
!! provided.
rcond, ferr, berr, work,rwork, info )
! -- lapack driver 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(inout) :: equed
character, intent(in) :: fact, uplo
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: lda, ldaf, ldb, ldx, n, nrhs
real(${ck}$), intent(out) :: rcond
! Array Arguments
real(${ck}$), intent(out) :: berr(*), ferr(*), rwork(*)
real(${ck}$), intent(inout) :: s(*)
complex(${ck}$), intent(inout) :: a(lda,*), af(ldaf,*), b(ldb,*)
complex(${ck}$), intent(out) :: work(*), x(ldx,*)
! =====================================================================
! Local Scalars
logical(lk) :: equil, nofact, rcequ
integer(${ik}$) :: i, infequ, j
real(${ck}$) :: amax, anorm, bignum, scond, smax, smin, smlnum
! Intrinsic Functions
! Executable Statements
info = 0_${ik}$
nofact = stdlib_lsame( fact, 'N' )
equil = stdlib_lsame( fact, 'E' )
if( nofact .or. equil ) then
equed = 'N'
rcequ = .false.
else
rcequ = stdlib_lsame( equed, 'Y' )
smlnum = stdlib${ii}$_${c2ri(ci)}$lamch( 'SAFE MINIMUM' )
bignum = one / smlnum
end if
! test the input parameters.
if( .not.nofact .and. .not.equil .and. .not.stdlib_lsame( fact, 'F' ) )then
info = -1_${ik}$
else if( .not.stdlib_lsame( uplo, 'U' ) .and. .not.stdlib_lsame( uplo, 'L' ) )&
then
info = -2_${ik}$
else if( n<0_${ik}$ ) then
info = -3_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -4_${ik}$
else if( lda<max( 1_${ik}$, n ) ) then
info = -6_${ik}$
else if( ldaf<max( 1_${ik}$, n ) ) then
info = -8_${ik}$
else if( stdlib_lsame( fact, 'F' ) .and. .not.( rcequ .or. stdlib_lsame( equed, 'N' ) )&
) then
info = -9_${ik}$
else
if( rcequ ) then
smin = bignum
smax = zero
do j = 1, n
smin = min( smin, s( j ) )
smax = max( smax, s( j ) )
end do
if( smin<=zero ) then
info = -10_${ik}$
else if( n>0_${ik}$ ) then
scond = max( smin, smlnum ) / min( smax, bignum )
else
scond = one
end if
end if
if( info==0_${ik}$ ) then
if( ldb<max( 1_${ik}$, n ) ) then
info = -12_${ik}$
else if( ldx<max( 1_${ik}$, n ) ) then
info = -14_${ik}$
end if
end if
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'ZPOSVX', -info )
return
end if
if( equil ) then
! compute row and column scalings to equilibrate the matrix a.
call stdlib${ii}$_${ci}$poequ( n, a, lda, s, scond, amax, infequ )
if( infequ==0_${ik}$ ) then
! equilibrate the matrix.
call stdlib${ii}$_${ci}$laqhe( uplo, n, a, lda, s, scond, amax, equed )
rcequ = stdlib_lsame( equed, 'Y' )
end if
end if
! scale the right hand side.
if( rcequ ) then
do j = 1, nrhs
do i = 1, n
b( i, j ) = s( i )*b( i, j )
end do
end do
end if
if( nofact .or. equil ) then
! compute the cholesky factorization a = u**h *u or a = l*l**h.
call stdlib${ii}$_${ci}$lacpy( uplo, n, n, a, lda, af, ldaf )
call stdlib${ii}$_${ci}$potrf( uplo, n, af, ldaf, info )
! return if info is non-zero.
if( info>0_${ik}$ )then
rcond = zero
return
end if
end if
! compute the norm of the matrix a.
anorm = stdlib${ii}$_${ci}$lanhe( '1', uplo, n, a, lda, rwork )
! compute the reciprocal of the condition number of a.
call stdlib${ii}$_${ci}$pocon( uplo, n, af, ldaf, anorm, rcond, work, rwork, info )
! compute the solution matrix x.
call stdlib${ii}$_${ci}$lacpy( 'FULL', n, nrhs, b, ldb, x, ldx )
call stdlib${ii}$_${ci}$potrs( uplo, n, nrhs, af, ldaf, x, ldx, info )
! use iterative refinement to improve the computed solution and
! compute error bounds and backward error estimates for it.
call stdlib${ii}$_${ci}$porfs( uplo, n, nrhs, a, lda, af, ldaf, b, ldb, x, ldx,ferr, berr, work, &
rwork, info )
! transform the solution matrix x to a solution of the original
! system.
if( rcequ ) then
do j = 1, nrhs
do i = 1, n
x( i, j ) = s( i )*x( i, j )
end do
end do
do j = 1, nrhs
ferr( j ) = ferr( j ) / scond
end do
end if
! set info = n+1 if the matrix is singular to working precision.
if( rcond<stdlib${ii}$_${c2ri(ci)}$lamch( 'EPSILON' ) )info = n + 1_${ik}$
return
end subroutine stdlib${ii}$_${ci}$posvx
#:endif
#:endfor
pure module subroutine stdlib${ii}$_sppsv( uplo, n, nrhs, ap, b, ldb, info )
!! SPPSV computes the solution to a real system of linear equations
!! A * X = B,
!! where A is an N-by-N symmetric positive definite matrix stored in
!! packed format and X and B are N-by-NRHS matrices.
!! The Cholesky decomposition is used to factor A as
!! A = U**T* U, if UPLO = 'U', or
!! A = L * L**T, if UPLO = 'L',
!! where U is an upper triangular matrix and L is a lower triangular
!! matrix. The factored form of A is then used to solve the system of
!! equations A * X = B.
! -- lapack driver 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
real(sp), intent(inout) :: ap(*), b(ldb,*)
! =====================================================================
! Intrinsic Functions
! Executable Statements
! test the input parameters.
info = 0_${ik}$
if( .not.stdlib_lsame( uplo, 'U' ) .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 = -6_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'SPPSV ', -info )
return
end if
! compute the cholesky factorization a = u**t*u or a = l*l**t.
call stdlib${ii}$_spptrf( uplo, n, ap, info )
if( info==0_${ik}$ ) then
! solve the system a*x = b, overwriting b with x.
call stdlib${ii}$_spptrs( uplo, n, nrhs, ap, b, ldb, info )
end if
return
end subroutine stdlib${ii}$_sppsv
pure module subroutine stdlib${ii}$_dppsv( uplo, n, nrhs, ap, b, ldb, info )
!! DPPSV computes the solution to a real system of linear equations
!! A * X = B,
!! where A is an N-by-N symmetric positive definite matrix stored in
!! packed format and X and B are N-by-NRHS matrices.
!! The Cholesky decomposition is used to factor A as
!! A = U**T* U, if UPLO = 'U', or
!! A = L * L**T, if UPLO = 'L',
!! where U is an upper triangular matrix and L is a lower triangular
!! matrix. The factored form of A is then used to solve the system of
!! equations A * X = B.
! -- lapack driver 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
real(dp), intent(inout) :: ap(*), b(ldb,*)
! =====================================================================
! Intrinsic Functions
! Executable Statements
! test the input parameters.
info = 0_${ik}$
if( .not.stdlib_lsame( uplo, 'U' ) .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 = -6_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'DPPSV ', -info )
return
end if
! compute the cholesky factorization a = u**t*u or a = l*l**t.
call stdlib${ii}$_dpptrf( uplo, n, ap, info )
if( info==0_${ik}$ ) then
! solve the system a*x = b, overwriting b with x.
call stdlib${ii}$_dpptrs( uplo, n, nrhs, ap, b, ldb, info )
end if
return
end subroutine stdlib${ii}$_dppsv
#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
pure module subroutine stdlib${ii}$_${ri}$ppsv( uplo, n, nrhs, ap, b, ldb, info )
!! DPPSV: computes the solution to a real system of linear equations
!! A * X = B,
!! where A is an N-by-N symmetric positive definite matrix stored in
!! packed format and X and B are N-by-NRHS matrices.
!! The Cholesky decomposition is used to factor A as
!! A = U**T* U, if UPLO = 'U', or
!! A = L * L**T, if UPLO = 'L',
!! where U is an upper triangular matrix and L is a lower triangular
!! matrix. The factored form of A is then used to solve the system of
!! equations A * X = B.
! -- lapack driver 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
real(${rk}$), intent(inout) :: ap(*), b(ldb,*)
! =====================================================================
! Intrinsic Functions
! Executable Statements
! test the input parameters.
info = 0_${ik}$
if( .not.stdlib_lsame( uplo, 'U' ) .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 = -6_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'DPPSV ', -info )
return
end if
! compute the cholesky factorization a = u**t*u or a = l*l**t.
call stdlib${ii}$_${ri}$pptrf( uplo, n, ap, info )
if( info==0_${ik}$ ) then
! solve the system a*x = b, overwriting b with x.
call stdlib${ii}$_${ri}$pptrs( uplo, n, nrhs, ap, b, ldb, info )
end if
return
end subroutine stdlib${ii}$_${ri}$ppsv
#:endif
#:endfor
pure module subroutine stdlib${ii}$_cppsv( uplo, n, nrhs, ap, b, ldb, info )
!! CPPSV computes the solution to a complex system of linear equations
!! A * X = B,
!! where A is an N-by-N Hermitian positive definite matrix stored in
!! packed format and X and B are N-by-NRHS matrices.
!! The Cholesky decomposition is used to factor A as
!! A = U**H * U, if UPLO = 'U', or
!! A = L * L**H, if UPLO = 'L',
!! where U is an upper triangular matrix and L is a lower triangular
!! matrix. The factored form of A is then used to solve the system of
!! equations A * X = B.
! -- lapack driver 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
complex(sp), intent(inout) :: ap(*), b(ldb,*)
! =====================================================================
! Intrinsic Functions
! Executable Statements
! test the input parameters.
info = 0_${ik}$
if( .not.stdlib_lsame( uplo, 'U' ) .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 = -6_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'CPPSV ', -info )
return
end if
! compute the cholesky factorization a = u**h *u or a = l*l**h.
call stdlib${ii}$_cpptrf( uplo, n, ap, info )
if( info==0_${ik}$ ) then
! solve the system a*x = b, overwriting b with x.
call stdlib${ii}$_cpptrs( uplo, n, nrhs, ap, b, ldb, info )
end if
return
end subroutine stdlib${ii}$_cppsv
pure module subroutine stdlib${ii}$_zppsv( uplo, n, nrhs, ap, b, ldb, info )
!! ZPPSV computes the solution to a complex system of linear equations
!! A * X = B,
!! where A is an N-by-N Hermitian positive definite matrix stored in
!! packed format and X and B are N-by-NRHS matrices.
!! The Cholesky decomposition is used to factor A as
!! A = U**H * U, if UPLO = 'U', or
!! A = L * L**H, if UPLO = 'L',
!! where U is an upper triangular matrix and L is a lower triangular
!! matrix. The factored form of A is then used to solve the system of
!! equations A * X = B.
! -- lapack driver 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
complex(dp), intent(inout) :: ap(*), b(ldb,*)
! =====================================================================
! Intrinsic Functions
! Executable Statements
! test the input parameters.
info = 0_${ik}$
if( .not.stdlib_lsame( uplo, 'U' ) .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 = -6_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'ZPPSV ', -info )
return
end if
! compute the cholesky factorization a = u**h *u or a = l*l**h.
call stdlib${ii}$_zpptrf( uplo, n, ap, info )
if( info==0_${ik}$ ) then
! solve the system a*x = b, overwriting b with x.
call stdlib${ii}$_zpptrs( uplo, n, nrhs, ap, b, ldb, info )
end if
return
end subroutine stdlib${ii}$_zppsv
#:for ck,ct,ci in CMPLX_KINDS_TYPES
#:if not ck in ["sp","dp"]
pure module subroutine stdlib${ii}$_${ci}$ppsv( uplo, n, nrhs, ap, b, ldb, info )
!! ZPPSV: computes the solution to a complex system of linear equations
!! A * X = B,
!! where A is an N-by-N Hermitian positive definite matrix stored in
!! packed format and X and B are N-by-NRHS matrices.
!! The Cholesky decomposition is used to factor A as
!! A = U**H * U, if UPLO = 'U', or
!! A = L * L**H, if UPLO = 'L',
!! where U is an upper triangular matrix and L is a lower triangular
!! matrix. The factored form of A is then used to solve the system of
!! equations A * X = B.
! -- lapack driver 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
complex(${ck}$), intent(inout) :: ap(*), b(ldb,*)
! =====================================================================
! Intrinsic Functions
! Executable Statements
! test the input parameters.
info = 0_${ik}$
if( .not.stdlib_lsame( uplo, 'U' ) .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 = -6_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'ZPPSV ', -info )
return
end if
! compute the cholesky factorization a = u**h *u or a = l*l**h.
call stdlib${ii}$_${ci}$pptrf( uplo, n, ap, info )
if( info==0_${ik}$ ) then
! solve the system a*x = b, overwriting b with x.
call stdlib${ii}$_${ci}$pptrs( uplo, n, nrhs, ap, b, ldb, info )
end if
return
end subroutine stdlib${ii}$_${ci}$ppsv
#:endif
#:endfor
module subroutine stdlib${ii}$_sppsvx( fact, uplo, n, nrhs, ap, afp, equed, s, b, ldb,x, ldx, rcond, ferr,&
!! SPPSVX uses the Cholesky factorization A = U**T*U or A = L*L**T to
!! compute the solution to a real system of linear equations
!! A * X = B,
!! where A is an N-by-N symmetric positive definite matrix stored in
!! packed format and X and B are N-by-NRHS matrices.
!! Error bounds on the solution and a condition estimate are also
!! provided.
berr, work, iwork, info )
! -- lapack driver 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(inout) :: equed
character, intent(in) :: fact, uplo
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: ldb, ldx, n, nrhs
real(sp), intent(out) :: rcond
! Array Arguments
integer(${ik}$), intent(out) :: iwork(*)
real(sp), intent(inout) :: afp(*), ap(*), b(ldb,*), s(*)
real(sp), intent(out) :: berr(*), ferr(*), work(*), x(ldx,*)
! =====================================================================
! Local Scalars
logical(lk) :: equil, nofact, rcequ
integer(${ik}$) :: i, infequ, j
real(sp) :: amax, anorm, bignum, scond, smax, smin, smlnum
! Intrinsic Functions
! Executable Statements
info = 0_${ik}$
nofact = stdlib_lsame( fact, 'N' )
equil = stdlib_lsame( fact, 'E' )
if( nofact .or. equil ) then
equed = 'N'
rcequ = .false.
else
rcequ = stdlib_lsame( equed, 'Y' )
smlnum = stdlib${ii}$_slamch( 'SAFE MINIMUM' )
bignum = one / smlnum
end if
! test the input parameters.
if( .not.nofact .and. .not.equil .and. .not.stdlib_lsame( fact, 'F' ) )then
info = -1_${ik}$
else if( .not.stdlib_lsame( uplo, 'U' ) .and. .not.stdlib_lsame( uplo, 'L' ) )&
then
info = -2_${ik}$
else if( n<0_${ik}$ ) then
info = -3_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -4_${ik}$
else if( stdlib_lsame( fact, 'F' ) .and. .not.( rcequ .or. stdlib_lsame( equed, 'N' ) )&
) then
info = -7_${ik}$
else
if( rcequ ) then
smin = bignum
smax = zero
do j = 1, n
smin = min( smin, s( j ) )
smax = max( smax, s( j ) )
end do
if( smin<=zero ) then
info = -8_${ik}$
else if( n>0_${ik}$ ) then
scond = max( smin, smlnum ) / min( smax, bignum )
else
scond = one
end if
end if
if( info==0_${ik}$ ) then
if( ldb<max( 1_${ik}$, n ) ) then
info = -10_${ik}$
else if( ldx<max( 1_${ik}$, n ) ) then
info = -12_${ik}$
end if
end if
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'SPPSVX', -info )
return
end if
if( equil ) then
! compute row and column scalings to equilibrate the matrix a.
call stdlib${ii}$_sppequ( uplo, n, ap, s, scond, amax, infequ )
if( infequ==0_${ik}$ ) then
! equilibrate the matrix.
call stdlib${ii}$_slaqsp( uplo, n, ap, s, scond, amax, equed )
rcequ = stdlib_lsame( equed, 'Y' )
end if
end if
! scale the right-hand side.
if( rcequ ) then
do j = 1, nrhs
do i = 1, n
b( i, j ) = s( i )*b( i, j )
end do
end do
end if
if( nofact .or. equil ) then
! compute the cholesky factorization a = u**t * u or a = l * l**t.
call stdlib${ii}$_scopy( n*( n+1 ) / 2_${ik}$, ap, 1_${ik}$, afp, 1_${ik}$ )
call stdlib${ii}$_spptrf( uplo, n, afp, info )
! return if info is non-zero.
if( info>0_${ik}$ )then
rcond = zero
return
end if
end if
! compute the norm of the matrix a.
anorm = stdlib${ii}$_slansp( 'I', uplo, n, ap, work )
! compute the reciprocal of the condition number of a.
call stdlib${ii}$_sppcon( uplo, n, afp, anorm, rcond, work, iwork, info )
! compute the solution matrix x.
call stdlib${ii}$_slacpy( 'FULL', n, nrhs, b, ldb, x, ldx )
call stdlib${ii}$_spptrs( uplo, n, nrhs, afp, x, ldx, info )
! use iterative refinement to improve the computed solution and
! compute error bounds and backward error estimates for it.
call stdlib${ii}$_spprfs( uplo, n, nrhs, ap, afp, b, ldb, x, ldx, ferr, berr,work, iwork, &
info )
! transform the solution matrix x to a solution of the original
! system.
if( rcequ ) then
do j = 1, nrhs
do i = 1, n
x( i, j ) = s( i )*x( i, j )
end do
end do
do j = 1, nrhs
ferr( j ) = ferr( j ) / scond
end do
end if
! set info = n+1 if the matrix is singular to working precision.
if( rcond<stdlib${ii}$_slamch( 'EPSILON' ) )info = n + 1_${ik}$
return
end subroutine stdlib${ii}$_sppsvx
module subroutine stdlib${ii}$_dppsvx( fact, uplo, n, nrhs, ap, afp, equed, s, b, ldb,x, ldx, rcond, ferr,&
!! DPPSVX uses the Cholesky factorization A = U**T*U or A = L*L**T to
!! compute the solution to a real system of linear equations
!! A * X = B,
!! where A is an N-by-N symmetric positive definite matrix stored in
!! packed format and X and B are N-by-NRHS matrices.
!! Error bounds on the solution and a condition estimate are also
!! provided.
berr, work, iwork, info )
! -- lapack driver 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(inout) :: equed
character, intent(in) :: fact, uplo
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: ldb, ldx, n, nrhs
real(dp), intent(out) :: rcond
! Array Arguments
integer(${ik}$), intent(out) :: iwork(*)
real(dp), intent(inout) :: afp(*), ap(*), b(ldb,*), s(*)
real(dp), intent(out) :: berr(*), ferr(*), work(*), x(ldx,*)
! =====================================================================
! Local Scalars
logical(lk) :: equil, nofact, rcequ
integer(${ik}$) :: i, infequ, j
real(dp) :: amax, anorm, bignum, scond, smax, smin, smlnum
! Intrinsic Functions
! Executable Statements
info = 0_${ik}$
nofact = stdlib_lsame( fact, 'N' )
equil = stdlib_lsame( fact, 'E' )
if( nofact .or. equil ) then
equed = 'N'
rcequ = .false.
else
rcequ = stdlib_lsame( equed, 'Y' )
smlnum = stdlib${ii}$_dlamch( 'SAFE MINIMUM' )
bignum = one / smlnum
end if
! test the input parameters.
if( .not.nofact .and. .not.equil .and. .not.stdlib_lsame( fact, 'F' ) )then
info = -1_${ik}$
else if( .not.stdlib_lsame( uplo, 'U' ) .and. .not.stdlib_lsame( uplo, 'L' ) )&
then
info = -2_${ik}$
else if( n<0_${ik}$ ) then
info = -3_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -4_${ik}$
else if( stdlib_lsame( fact, 'F' ) .and. .not.( rcequ .or. stdlib_lsame( equed, 'N' ) )&
) then
info = -7_${ik}$
else
if( rcequ ) then
smin = bignum
smax = zero
do j = 1, n
smin = min( smin, s( j ) )
smax = max( smax, s( j ) )
end do
if( smin<=zero ) then
info = -8_${ik}$
else if( n>0_${ik}$ ) then
scond = max( smin, smlnum ) / min( smax, bignum )
else
scond = one
end if
end if
if( info==0_${ik}$ ) then
if( ldb<max( 1_${ik}$, n ) ) then
info = -10_${ik}$
else if( ldx<max( 1_${ik}$, n ) ) then
info = -12_${ik}$
end if
end if
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'DPPSVX', -info )
return
end if
if( equil ) then
! compute row and column scalings to equilibrate the matrix a.
call stdlib${ii}$_dppequ( uplo, n, ap, s, scond, amax, infequ )
if( infequ==0_${ik}$ ) then
! equilibrate the matrix.
call stdlib${ii}$_dlaqsp( uplo, n, ap, s, scond, amax, equed )
rcequ = stdlib_lsame( equed, 'Y' )
end if
end if
! scale the right-hand side.
if( rcequ ) then
do j = 1, nrhs
do i = 1, n
b( i, j ) = s( i )*b( i, j )
end do
end do
end if
if( nofact .or. equil ) then
! compute the cholesky factorization a = u**t * u or a = l * l**t.
call stdlib${ii}$_dcopy( n*( n+1 ) / 2_${ik}$, ap, 1_${ik}$, afp, 1_${ik}$ )
call stdlib${ii}$_dpptrf( uplo, n, afp, info )
! return if info is non-zero.
if( info>0_${ik}$ )then
rcond = zero
return
end if
end if
! compute the norm of the matrix a.
anorm = stdlib${ii}$_dlansp( 'I', uplo, n, ap, work )
! compute the reciprocal of the condition number of a.
call stdlib${ii}$_dppcon( uplo, n, afp, anorm, rcond, work, iwork, info )
! compute the solution matrix x.
call stdlib${ii}$_dlacpy( 'FULL', n, nrhs, b, ldb, x, ldx )
call stdlib${ii}$_dpptrs( uplo, n, nrhs, afp, x, ldx, info )
! use iterative refinement to improve the computed solution and
! compute error bounds and backward error estimates for it.
call stdlib${ii}$_dpprfs( uplo, n, nrhs, ap, afp, b, ldb, x, ldx, ferr, berr,work, iwork, &
info )
! transform the solution matrix x to a solution of the original
! system.
if( rcequ ) then
do j = 1, nrhs
do i = 1, n
x( i, j ) = s( i )*x( i, j )
end do
end do
do j = 1, nrhs
ferr( j ) = ferr( j ) / scond
end do
end if
! set info = n+1 if the matrix is singular to working precision.
if( rcond<stdlib${ii}$_dlamch( 'EPSILON' ) )info = n + 1_${ik}$
return
end subroutine stdlib${ii}$_dppsvx
#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
module subroutine stdlib${ii}$_${ri}$ppsvx( fact, uplo, n, nrhs, ap, afp, equed, s, b, ldb,x, ldx, rcond, ferr,&
!! DPPSVX: uses the Cholesky factorization A = U**T*U or A = L*L**T to
!! compute the solution to a real system of linear equations
!! A * X = B,
!! where A is an N-by-N symmetric positive definite matrix stored in
!! packed format and X and B are N-by-NRHS matrices.
!! Error bounds on the solution and a condition estimate are also
!! provided.
berr, work, iwork, info )
! -- lapack driver 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(inout) :: equed
character, intent(in) :: fact, uplo
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: ldb, ldx, n, nrhs
real(${rk}$), intent(out) :: rcond
! Array Arguments
integer(${ik}$), intent(out) :: iwork(*)
real(${rk}$), intent(inout) :: afp(*), ap(*), b(ldb,*), s(*)
real(${rk}$), intent(out) :: berr(*), ferr(*), work(*), x(ldx,*)
! =====================================================================
! Local Scalars
logical(lk) :: equil, nofact, rcequ
integer(${ik}$) :: i, infequ, j
real(${rk}$) :: amax, anorm, bignum, scond, smax, smin, smlnum
! Intrinsic Functions
! Executable Statements
info = 0_${ik}$
nofact = stdlib_lsame( fact, 'N' )
equil = stdlib_lsame( fact, 'E' )
if( nofact .or. equil ) then
equed = 'N'
rcequ = .false.
else
rcequ = stdlib_lsame( equed, 'Y' )
smlnum = stdlib${ii}$_${ri}$lamch( 'SAFE MINIMUM' )
bignum = one / smlnum
end if
! test the input parameters.
if( .not.nofact .and. .not.equil .and. .not.stdlib_lsame( fact, 'F' ) )then
info = -1_${ik}$
else if( .not.stdlib_lsame( uplo, 'U' ) .and. .not.stdlib_lsame( uplo, 'L' ) )&
then
info = -2_${ik}$
else if( n<0_${ik}$ ) then
info = -3_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -4_${ik}$
else if( stdlib_lsame( fact, 'F' ) .and. .not.( rcequ .or. stdlib_lsame( equed, 'N' ) )&
) then
info = -7_${ik}$
else
if( rcequ ) then
smin = bignum
smax = zero
do j = 1, n
smin = min( smin, s( j ) )
smax = max( smax, s( j ) )
end do
if( smin<=zero ) then
info = -8_${ik}$
else if( n>0_${ik}$ ) then
scond = max( smin, smlnum ) / min( smax, bignum )
else
scond = one
end if
end if
if( info==0_${ik}$ ) then
if( ldb<max( 1_${ik}$, n ) ) then
info = -10_${ik}$
else if( ldx<max( 1_${ik}$, n ) ) then
info = -12_${ik}$
end if
end if
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'DPPSVX', -info )
return
end if
if( equil ) then
! compute row and column scalings to equilibrate the matrix a.
call stdlib${ii}$_${ri}$ppequ( uplo, n, ap, s, scond, amax, infequ )
if( infequ==0_${ik}$ ) then
! equilibrate the matrix.
call stdlib${ii}$_${ri}$laqsp( uplo, n, ap, s, scond, amax, equed )
rcequ = stdlib_lsame( equed, 'Y' )
end if
end if
! scale the right-hand side.
if( rcequ ) then
do j = 1, nrhs
do i = 1, n
b( i, j ) = s( i )*b( i, j )
end do
end do
end if
if( nofact .or. equil ) then
! compute the cholesky factorization a = u**t * u or a = l * l**t.
call stdlib${ii}$_${ri}$copy( n*( n+1 ) / 2_${ik}$, ap, 1_${ik}$, afp, 1_${ik}$ )
call stdlib${ii}$_${ri}$pptrf( uplo, n, afp, info )
! return if info is non-zero.
if( info>0_${ik}$ )then
rcond = zero
return
end if
end if
! compute the norm of the matrix a.
anorm = stdlib${ii}$_${ri}$lansp( 'I', uplo, n, ap, work )
! compute the reciprocal of the condition number of a.
call stdlib${ii}$_${ri}$ppcon( uplo, n, afp, anorm, rcond, work, iwork, info )
! compute the solution matrix x.
call stdlib${ii}$_${ri}$lacpy( 'FULL', n, nrhs, b, ldb, x, ldx )
call stdlib${ii}$_${ri}$pptrs( uplo, n, nrhs, afp, x, ldx, info )
! use iterative refinement to improve the computed solution and
! compute error bounds and backward error estimates for it.
call stdlib${ii}$_${ri}$pprfs( uplo, n, nrhs, ap, afp, b, ldb, x, ldx, ferr, berr,work, iwork, &
info )
! transform the solution matrix x to a solution of the original
! system.
if( rcequ ) then
do j = 1, nrhs
do i = 1, n
x( i, j ) = s( i )*x( i, j )
end do
end do
do j = 1, nrhs
ferr( j ) = ferr( j ) / scond
end do
end if
! set info = n+1 if the matrix is singular to working precision.
if( rcond<stdlib${ii}$_${ri}$lamch( 'EPSILON' ) )info = n + 1_${ik}$
return
end subroutine stdlib${ii}$_${ri}$ppsvx
#:endif
#:endfor
module subroutine stdlib${ii}$_cppsvx( fact, uplo, n, nrhs, ap, afp, equed, s, b, ldb,x, ldx, rcond, ferr,&
!! CPPSVX uses the Cholesky factorization A = U**H*U or A = L*L**H to
!! compute the solution to a complex system of linear equations
!! A * X = B,
!! where A is an N-by-N Hermitian positive definite matrix stored in
!! packed format and X and B are N-by-NRHS matrices.
!! Error bounds on the solution and a condition estimate are also
!! provided.
berr, work, rwork, info )
! -- lapack driver 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(inout) :: equed
character, intent(in) :: fact, uplo
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: ldb, ldx, n, nrhs
real(sp), intent(out) :: rcond
! Array Arguments
real(sp), intent(out) :: berr(*), ferr(*), rwork(*)
real(sp), intent(inout) :: s(*)
complex(sp), intent(inout) :: afp(*), ap(*), b(ldb,*)
complex(sp), intent(out) :: work(*), x(ldx,*)
! =====================================================================
! Local Scalars
logical(lk) :: equil, nofact, rcequ
integer(${ik}$) :: i, infequ, j
real(sp) :: amax, anorm, bignum, scond, smax, smin, smlnum
! Intrinsic Functions
! Executable Statements
info = 0_${ik}$
nofact = stdlib_lsame( fact, 'N' )
equil = stdlib_lsame( fact, 'E' )
if( nofact .or. equil ) then
equed = 'N'
rcequ = .false.
else
rcequ = stdlib_lsame( equed, 'Y' )
smlnum = stdlib${ii}$_slamch( 'SAFE MINIMUM' )
bignum = one / smlnum
end if
! test the input parameters.
if( .not.nofact .and. .not.equil .and. .not.stdlib_lsame( fact, 'F' ) )then
info = -1_${ik}$
else if( .not.stdlib_lsame( uplo, 'U' ) .and. .not.stdlib_lsame( uplo, 'L' ) )&
then
info = -2_${ik}$
else if( n<0_${ik}$ ) then
info = -3_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -4_${ik}$
else if( stdlib_lsame( fact, 'F' ) .and. .not.( rcequ .or. stdlib_lsame( equed, 'N' ) )&
) then
info = -7_${ik}$
else
if( rcequ ) then
smin = bignum
smax = zero
do j = 1, n
smin = min( smin, s( j ) )
smax = max( smax, s( j ) )
end do
if( smin<=zero ) then
info = -8_${ik}$
else if( n>0_${ik}$ ) then
scond = max( smin, smlnum ) / min( smax, bignum )
else
scond = one
end if
end if
if( info==0_${ik}$ ) then
if( ldb<max( 1_${ik}$, n ) ) then
info = -10_${ik}$
else if( ldx<max( 1_${ik}$, n ) ) then
info = -12_${ik}$
end if
end if
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'CPPSVX', -info )
return
end if
if( equil ) then
! compute row and column scalings to equilibrate the matrix a.
call stdlib${ii}$_cppequ( uplo, n, ap, s, scond, amax, infequ )
if( infequ==0_${ik}$ ) then
! equilibrate the matrix.
call stdlib${ii}$_claqhp( uplo, n, ap, s, scond, amax, equed )
rcequ = stdlib_lsame( equed, 'Y' )
end if
end if
! scale the right-hand side.
if( rcequ ) then
do j = 1, nrhs
do i = 1, n
b( i, j ) = s( i )*b( i, j )
end do
end do
end if
if( nofact .or. equil ) then
! compute the cholesky factorization a = u**h * u or a = l * l**h.
call stdlib${ii}$_ccopy( n*( n+1 ) / 2_${ik}$, ap, 1_${ik}$, afp, 1_${ik}$ )
call stdlib${ii}$_cpptrf( uplo, n, afp, info )
! return if info is non-zero.
if( info>0_${ik}$ )then
rcond = zero
return
end if
end if
! compute the norm of the matrix a.
anorm = stdlib${ii}$_clanhp( 'I', uplo, n, ap, rwork )
! compute the reciprocal of the condition number of a.
call stdlib${ii}$_cppcon( uplo, n, afp, anorm, rcond, work, rwork, info )
! compute the solution matrix x.
call stdlib${ii}$_clacpy( 'FULL', n, nrhs, b, ldb, x, ldx )
call stdlib${ii}$_cpptrs( uplo, n, nrhs, afp, x, ldx, info )
! use iterative refinement to improve the computed solution and
! compute error bounds and backward error estimates for it.
call stdlib${ii}$_cpprfs( uplo, n, nrhs, ap, afp, b, ldb, x, ldx, ferr, berr,work, rwork, &
info )
! transform the solution matrix x to a solution of the original
! system.
if( rcequ ) then
do j = 1, nrhs
do i = 1, n
x( i, j ) = s( i )*x( i, j )
end do
end do
do j = 1, nrhs
ferr( j ) = ferr( j ) / scond
end do
end if
! set info = n+1 if the matrix is singular to working precision.
if( rcond<stdlib${ii}$_slamch( 'EPSILON' ) )info = n + 1_${ik}$
return
end subroutine stdlib${ii}$_cppsvx
module subroutine stdlib${ii}$_zppsvx( fact, uplo, n, nrhs, ap, afp, equed, s, b, ldb,x, ldx, rcond, ferr,&
!! ZPPSVX uses the Cholesky factorization A = U**H * U or A = L * L**H to
!! compute the solution to a complex system of linear equations
!! A * X = B,
!! where A is an N-by-N Hermitian positive definite matrix stored in
!! packed format and X and B are N-by-NRHS matrices.
!! Error bounds on the solution and a condition estimate are also
!! provided.
berr, work, rwork, info )
! -- lapack driver 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(inout) :: equed
character, intent(in) :: fact, uplo
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: ldb, ldx, n, nrhs
real(dp), intent(out) :: rcond
! Array Arguments
real(dp), intent(out) :: berr(*), ferr(*), rwork(*)
real(dp), intent(inout) :: s(*)
complex(dp), intent(inout) :: afp(*), ap(*), b(ldb,*)
complex(dp), intent(out) :: work(*), x(ldx,*)
! =====================================================================
! Local Scalars
logical(lk) :: equil, nofact, rcequ
integer(${ik}$) :: i, infequ, j
real(dp) :: amax, anorm, bignum, scond, smax, smin, smlnum
! Intrinsic Functions
! Executable Statements
info = 0_${ik}$
nofact = stdlib_lsame( fact, 'N' )
equil = stdlib_lsame( fact, 'E' )
if( nofact .or. equil ) then
equed = 'N'
rcequ = .false.
else
rcequ = stdlib_lsame( equed, 'Y' )
smlnum = stdlib${ii}$_dlamch( 'SAFE MINIMUM' )
bignum = one / smlnum
end if
! test the input parameters.
if( .not.nofact .and. .not.equil .and. .not.stdlib_lsame( fact, 'F' ) )then
info = -1_${ik}$
else if( .not.stdlib_lsame( uplo, 'U' ) .and. .not.stdlib_lsame( uplo, 'L' ) )&
then
info = -2_${ik}$
else if( n<0_${ik}$ ) then
info = -3_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -4_${ik}$
else if( stdlib_lsame( fact, 'F' ) .and. .not.( rcequ .or. stdlib_lsame( equed, 'N' ) )&
) then
info = -7_${ik}$
else
if( rcequ ) then
smin = bignum
smax = zero
do j = 1, n
smin = min( smin, s( j ) )
smax = max( smax, s( j ) )
end do
if( smin<=zero ) then
info = -8_${ik}$
else if( n>0_${ik}$ ) then
scond = max( smin, smlnum ) / min( smax, bignum )
else
scond = one
end if
end if
if( info==0_${ik}$ ) then
if( ldb<max( 1_${ik}$, n ) ) then
info = -10_${ik}$
else if( ldx<max( 1_${ik}$, n ) ) then
info = -12_${ik}$
end if
end if
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'ZPPSVX', -info )
return
end if
if( equil ) then
! compute row and column scalings to equilibrate the matrix a.
call stdlib${ii}$_zppequ( uplo, n, ap, s, scond, amax, infequ )
if( infequ==0_${ik}$ ) then
! equilibrate the matrix.
call stdlib${ii}$_zlaqhp( uplo, n, ap, s, scond, amax, equed )
rcequ = stdlib_lsame( equed, 'Y' )
end if
end if
! scale the right-hand side.
if( rcequ ) then
do j = 1, nrhs
do i = 1, n
b( i, j ) = s( i )*b( i, j )
end do
end do
end if
if( nofact .or. equil ) then
! compute the cholesky factorization a = u**h * u or a = l * l**h.
call stdlib${ii}$_zcopy( n*( n+1 ) / 2_${ik}$, ap, 1_${ik}$, afp, 1_${ik}$ )
call stdlib${ii}$_zpptrf( uplo, n, afp, info )
! return if info is non-zero.
if( info>0_${ik}$ )then
rcond = zero
return
end if
end if
! compute the norm of the matrix a.
anorm = stdlib${ii}$_zlanhp( 'I', uplo, n, ap, rwork )
! compute the reciprocal of the condition number of a.
call stdlib${ii}$_zppcon( uplo, n, afp, anorm, rcond, work, rwork, info )
! compute the solution matrix x.
call stdlib${ii}$_zlacpy( 'FULL', n, nrhs, b, ldb, x, ldx )
call stdlib${ii}$_zpptrs( uplo, n, nrhs, afp, x, ldx, info )
! use iterative refinement to improve the computed solution and
! compute error bounds and backward error estimates for it.
call stdlib${ii}$_zpprfs( uplo, n, nrhs, ap, afp, b, ldb, x, ldx, ferr, berr,work, rwork, &
info )
! transform the solution matrix x to a solution of the original
! system.
if( rcequ ) then
do j = 1, nrhs
do i = 1, n
x( i, j ) = s( i )*x( i, j )
end do
end do
do j = 1, nrhs
ferr( j ) = ferr( j ) / scond
end do
end if
! set info = n+1 if the matrix is singular to working precision.
if( rcond<stdlib${ii}$_dlamch( 'EPSILON' ) )info = n + 1_${ik}$
return
end subroutine stdlib${ii}$_zppsvx
#:for ck,ct,ci in CMPLX_KINDS_TYPES
#:if not ck in ["sp","dp"]
module subroutine stdlib${ii}$_${ci}$ppsvx( fact, uplo, n, nrhs, ap, afp, equed, s, b, ldb,x, ldx, rcond, ferr,&
!! ZPPSVX: uses the Cholesky factorization A = U**H * U or A = L * L**H to
!! compute the solution to a complex system of linear equations
!! A * X = B,
!! where A is an N-by-N Hermitian positive definite matrix stored in
!! packed format and X and B are N-by-NRHS matrices.
!! Error bounds on the solution and a condition estimate are also
!! provided.
berr, work, rwork, info )
! -- lapack driver 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(inout) :: equed
character, intent(in) :: fact, uplo
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: ldb, ldx, n, nrhs
real(${ck}$), intent(out) :: rcond
! Array Arguments
real(${ck}$), intent(out) :: berr(*), ferr(*), rwork(*)
real(${ck}$), intent(inout) :: s(*)
complex(${ck}$), intent(inout) :: afp(*), ap(*), b(ldb,*)
complex(${ck}$), intent(out) :: work(*), x(ldx,*)
! =====================================================================
! Local Scalars
logical(lk) :: equil, nofact, rcequ
integer(${ik}$) :: i, infequ, j
real(${ck}$) :: amax, anorm, bignum, scond, smax, smin, smlnum
! Intrinsic Functions
! Executable Statements
info = 0_${ik}$
nofact = stdlib_lsame( fact, 'N' )
equil = stdlib_lsame( fact, 'E' )
if( nofact .or. equil ) then
equed = 'N'
rcequ = .false.
else
rcequ = stdlib_lsame( equed, 'Y' )
smlnum = stdlib${ii}$_${c2ri(ci)}$lamch( 'SAFE MINIMUM' )
bignum = one / smlnum
end if
! test the input parameters.
if( .not.nofact .and. .not.equil .and. .not.stdlib_lsame( fact, 'F' ) )then
info = -1_${ik}$
else if( .not.stdlib_lsame( uplo, 'U' ) .and. .not.stdlib_lsame( uplo, 'L' ) )&
then
info = -2_${ik}$
else if( n<0_${ik}$ ) then
info = -3_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -4_${ik}$
else if( stdlib_lsame( fact, 'F' ) .and. .not.( rcequ .or. stdlib_lsame( equed, 'N' ) )&
) then
info = -7_${ik}$
else
if( rcequ ) then
smin = bignum
smax = zero
do j = 1, n
smin = min( smin, s( j ) )
smax = max( smax, s( j ) )
end do
if( smin<=zero ) then
info = -8_${ik}$
else if( n>0_${ik}$ ) then
scond = max( smin, smlnum ) / min( smax, bignum )
else
scond = one
end if
end if
if( info==0_${ik}$ ) then
if( ldb<max( 1_${ik}$, n ) ) then
info = -10_${ik}$
else if( ldx<max( 1_${ik}$, n ) ) then
info = -12_${ik}$
end if
end if
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'ZPPSVX', -info )
return
end if
if( equil ) then
! compute row and column scalings to equilibrate the matrix a.
call stdlib${ii}$_${ci}$ppequ( uplo, n, ap, s, scond, amax, infequ )
if( infequ==0_${ik}$ ) then
! equilibrate the matrix.
call stdlib${ii}$_${ci}$laqhp( uplo, n, ap, s, scond, amax, equed )
rcequ = stdlib_lsame( equed, 'Y' )
end if
end if
! scale the right-hand side.
if( rcequ ) then
do j = 1, nrhs
do i = 1, n
b( i, j ) = s( i )*b( i, j )
end do
end do
end if
if( nofact .or. equil ) then
! compute the cholesky factorization a = u**h * u or a = l * l**h.
call stdlib${ii}$_${ci}$copy( n*( n+1 ) / 2_${ik}$, ap, 1_${ik}$, afp, 1_${ik}$ )
call stdlib${ii}$_${ci}$pptrf( uplo, n, afp, info )
! return if info is non-zero.
if( info>0_${ik}$ )then
rcond = zero
return
end if
end if
! compute the norm of the matrix a.
anorm = stdlib${ii}$_${ci}$lanhp( 'I', uplo, n, ap, rwork )
! compute the reciprocal of the condition number of a.
call stdlib${ii}$_${ci}$ppcon( uplo, n, afp, anorm, rcond, work, rwork, info )
! compute the solution matrix x.
call stdlib${ii}$_${ci}$lacpy( 'FULL', n, nrhs, b, ldb, x, ldx )
call stdlib${ii}$_${ci}$pptrs( uplo, n, nrhs, afp, x, ldx, info )
! use iterative refinement to improve the computed solution and
! compute error bounds and backward error estimates for it.
call stdlib${ii}$_${ci}$pprfs( uplo, n, nrhs, ap, afp, b, ldb, x, ldx, ferr, berr,work, rwork, &
info )
! transform the solution matrix x to a solution of the original
! system.
if( rcequ ) then
do j = 1, nrhs
do i = 1, n
x( i, j ) = s( i )*x( i, j )
end do
end do
do j = 1, nrhs
ferr( j ) = ferr( j ) / scond
end do
end if
! set info = n+1 if the matrix is singular to working precision.
if( rcond<stdlib${ii}$_${c2ri(ci)}$lamch( 'EPSILON' ) )info = n + 1_${ik}$
return
end subroutine stdlib${ii}$_${ci}$ppsvx
#:endif
#:endfor
pure module subroutine stdlib${ii}$_spbsv( uplo, n, kd, nrhs, ab, ldab, b, ldb, info )
!! SPBSV computes the solution to a real system of linear equations
!! A * X = B,
!! where A is an N-by-N symmetric positive definite band matrix and X
!! and B are N-by-NRHS matrices.
!! The Cholesky decomposition is used to factor A as
!! A = U**T * U, if UPLO = 'U', or
!! A = L * L**T, if UPLO = 'L',
!! where U is an upper triangular band matrix, and L is a lower
!! triangular band matrix, with the same number of superdiagonals or
!! subdiagonals as A. The factored form of A is then used to solve the
!! system of equations A * X = B.
! -- lapack driver 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) :: kd, ldab, ldb, n, nrhs
! Array Arguments
real(sp), intent(inout) :: ab(ldab,*), b(ldb,*)
! =====================================================================
! Intrinsic Functions
! Executable Statements
! test the input parameters.
info = 0_${ik}$
if( .not.stdlib_lsame( uplo, 'U' ) .and. .not.stdlib_lsame( uplo, 'L' ) ) then
info = -1_${ik}$
else if( n<0_${ik}$ ) then
info = -2_${ik}$
else if( kd<0_${ik}$ ) then
info = -3_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -4_${ik}$
else if( ldab<kd+1 ) then
info = -6_${ik}$
else if( ldb<max( 1_${ik}$, n ) ) then
info = -8_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'SPBSV ', -info )
return
end if
! compute the cholesky factorization a = u**t*u or a = l*l**t.
call stdlib${ii}$_spbtrf( uplo, n, kd, ab, ldab, info )
if( info==0_${ik}$ ) then
! solve the system a*x = b, overwriting b with x.
call stdlib${ii}$_spbtrs( uplo, n, kd, nrhs, ab, ldab, b, ldb, info )
end if
return
end subroutine stdlib${ii}$_spbsv
pure module subroutine stdlib${ii}$_dpbsv( uplo, n, kd, nrhs, ab, ldab, b, ldb, info )
!! DPBSV computes the solution to a real system of linear equations
!! A * X = B,
!! where A is an N-by-N symmetric positive definite band matrix and X
!! and B are N-by-NRHS matrices.
!! The Cholesky decomposition is used to factor A as
!! A = U**T * U, if UPLO = 'U', or
!! A = L * L**T, if UPLO = 'L',
!! where U is an upper triangular band matrix, and L is a lower
!! triangular band matrix, with the same number of superdiagonals or
!! subdiagonals as A. The factored form of A is then used to solve the
!! system of equations A * X = B.
! -- lapack driver 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) :: kd, ldab, ldb, n, nrhs
! Array Arguments
real(dp), intent(inout) :: ab(ldab,*), b(ldb,*)
! =====================================================================
! Intrinsic Functions
! Executable Statements
! test the input parameters.
info = 0_${ik}$
if( .not.stdlib_lsame( uplo, 'U' ) .and. .not.stdlib_lsame( uplo, 'L' ) ) then
info = -1_${ik}$
else if( n<0_${ik}$ ) then
info = -2_${ik}$
else if( kd<0_${ik}$ ) then
info = -3_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -4_${ik}$
else if( ldab<kd+1 ) then
info = -6_${ik}$
else if( ldb<max( 1_${ik}$, n ) ) then
info = -8_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'DPBSV ', -info )
return
end if
! compute the cholesky factorization a = u**t*u or a = l*l**t.
call stdlib${ii}$_dpbtrf( uplo, n, kd, ab, ldab, info )
if( info==0_${ik}$ ) then
! solve the system a*x = b, overwriting b with x.
call stdlib${ii}$_dpbtrs( uplo, n, kd, nrhs, ab, ldab, b, ldb, info )
end if
return
end subroutine stdlib${ii}$_dpbsv
#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
pure module subroutine stdlib${ii}$_${ri}$pbsv( uplo, n, kd, nrhs, ab, ldab, b, ldb, info )
!! DPBSV: computes the solution to a real system of linear equations
!! A * X = B,
!! where A is an N-by-N symmetric positive definite band matrix and X
!! and B are N-by-NRHS matrices.
!! The Cholesky decomposition is used to factor A as
!! A = U**T * U, if UPLO = 'U', or
!! A = L * L**T, if UPLO = 'L',
!! where U is an upper triangular band matrix, and L is a lower
!! triangular band matrix, with the same number of superdiagonals or
!! subdiagonals as A. The factored form of A is then used to solve the
!! system of equations A * X = B.
! -- lapack driver 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) :: kd, ldab, ldb, n, nrhs
! Array Arguments
real(${rk}$), intent(inout) :: ab(ldab,*), b(ldb,*)
! =====================================================================
! Intrinsic Functions
! Executable Statements
! test the input parameters.
info = 0_${ik}$
if( .not.stdlib_lsame( uplo, 'U' ) .and. .not.stdlib_lsame( uplo, 'L' ) ) then
info = -1_${ik}$
else if( n<0_${ik}$ ) then
info = -2_${ik}$
else if( kd<0_${ik}$ ) then
info = -3_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -4_${ik}$
else if( ldab<kd+1 ) then
info = -6_${ik}$
else if( ldb<max( 1_${ik}$, n ) ) then
info = -8_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'DPBSV ', -info )
return
end if
! compute the cholesky factorization a = u**t*u or a = l*l**t.
call stdlib${ii}$_${ri}$pbtrf( uplo, n, kd, ab, ldab, info )
if( info==0_${ik}$ ) then
! solve the system a*x = b, overwriting b with x.
call stdlib${ii}$_${ri}$pbtrs( uplo, n, kd, nrhs, ab, ldab, b, ldb, info )
end if
return
end subroutine stdlib${ii}$_${ri}$pbsv
#:endif
#:endfor
pure module subroutine stdlib${ii}$_cpbsv( uplo, n, kd, nrhs, ab, ldab, b, ldb, info )
!! CPBSV computes the solution to a complex system of linear equations
!! A * X = B,
!! where A is an N-by-N Hermitian positive definite band matrix and X
!! and B are N-by-NRHS matrices.
!! The Cholesky decomposition is used to factor A as
!! A = U**H * U, if UPLO = 'U', or
!! A = L * L**H, if UPLO = 'L',
!! where U is an upper triangular band matrix, and L is a lower
!! triangular band matrix, with the same number of superdiagonals or
!! subdiagonals as A. The factored form of A is then used to solve the
!! system of equations A * X = B.
! -- lapack driver 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) :: kd, ldab, ldb, n, nrhs
! Array Arguments
complex(sp), intent(inout) :: ab(ldab,*), b(ldb,*)
! =====================================================================
! Intrinsic Functions
! Executable Statements
! test the input parameters.
info = 0_${ik}$
if( .not.stdlib_lsame( uplo, 'U' ) .and. .not.stdlib_lsame( uplo, 'L' ) ) then
info = -1_${ik}$
else if( n<0_${ik}$ ) then
info = -2_${ik}$
else if( kd<0_${ik}$ ) then
info = -3_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -4_${ik}$
else if( ldab<kd+1 ) then
info = -6_${ik}$
else if( ldb<max( 1_${ik}$, n ) ) then
info = -8_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'CPBSV ', -info )
return
end if
! compute the cholesky factorization a = u**h*u or a = l*l**h.
call stdlib${ii}$_cpbtrf( uplo, n, kd, ab, ldab, info )
if( info==0_${ik}$ ) then
! solve the system a*x = b, overwriting b with x.
call stdlib${ii}$_cpbtrs( uplo, n, kd, nrhs, ab, ldab, b, ldb, info )
end if
return
end subroutine stdlib${ii}$_cpbsv
pure module subroutine stdlib${ii}$_zpbsv( uplo, n, kd, nrhs, ab, ldab, b, ldb, info )
!! ZPBSV computes the solution to a complex system of linear equations
!! A * X = B,
!! where A is an N-by-N Hermitian positive definite band matrix and X
!! and B are N-by-NRHS matrices.
!! The Cholesky decomposition is used to factor A as
!! A = U**H * U, if UPLO = 'U', or
!! A = L * L**H, if UPLO = 'L',
!! where U is an upper triangular band matrix, and L is a lower
!! triangular band matrix, with the same number of superdiagonals or
!! subdiagonals as A. The factored form of A is then used to solve the
!! system of equations A * X = B.
! -- lapack driver 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) :: kd, ldab, ldb, n, nrhs
! Array Arguments
complex(dp), intent(inout) :: ab(ldab,*), b(ldb,*)
! =====================================================================
! Intrinsic Functions
! Executable Statements
! test the input parameters.
info = 0_${ik}$
if( .not.stdlib_lsame( uplo, 'U' ) .and. .not.stdlib_lsame( uplo, 'L' ) ) then
info = -1_${ik}$
else if( n<0_${ik}$ ) then
info = -2_${ik}$
else if( kd<0_${ik}$ ) then
info = -3_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -4_${ik}$
else if( ldab<kd+1 ) then
info = -6_${ik}$
else if( ldb<max( 1_${ik}$, n ) ) then
info = -8_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'ZPBSV ', -info )
return
end if
! compute the cholesky factorization a = u**h *u or a = l*l**h.
call stdlib${ii}$_zpbtrf( uplo, n, kd, ab, ldab, info )
if( info==0_${ik}$ ) then
! solve the system a*x = b, overwriting b with x.
call stdlib${ii}$_zpbtrs( uplo, n, kd, nrhs, ab, ldab, b, ldb, info )
end if
return
end subroutine stdlib${ii}$_zpbsv
#:for ck,ct,ci in CMPLX_KINDS_TYPES
#:if not ck in ["sp","dp"]
pure module subroutine stdlib${ii}$_${ci}$pbsv( uplo, n, kd, nrhs, ab, ldab, b, ldb, info )
!! ZPBSV: computes the solution to a complex system of linear equations
!! A * X = B,
!! where A is an N-by-N Hermitian positive definite band matrix and X
!! and B are N-by-NRHS matrices.
!! The Cholesky decomposition is used to factor A as
!! A = U**H * U, if UPLO = 'U', or
!! A = L * L**H, if UPLO = 'L',
!! where U is an upper triangular band matrix, and L is a lower
!! triangular band matrix, with the same number of superdiagonals or
!! subdiagonals as A. The factored form of A is then used to solve the
!! system of equations A * X = B.
! -- lapack driver 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) :: kd, ldab, ldb, n, nrhs
! Array Arguments
complex(${ck}$), intent(inout) :: ab(ldab,*), b(ldb,*)
! =====================================================================
! Intrinsic Functions
! Executable Statements
! test the input parameters.
info = 0_${ik}$
if( .not.stdlib_lsame( uplo, 'U' ) .and. .not.stdlib_lsame( uplo, 'L' ) ) then
info = -1_${ik}$
else if( n<0_${ik}$ ) then
info = -2_${ik}$
else if( kd<0_${ik}$ ) then
info = -3_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -4_${ik}$
else if( ldab<kd+1 ) then
info = -6_${ik}$
else if( ldb<max( 1_${ik}$, n ) ) then
info = -8_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'ZPBSV ', -info )
return
end if
! compute the cholesky factorization a = u**h *u or a = l*l**h.
call stdlib${ii}$_${ci}$pbtrf( uplo, n, kd, ab, ldab, info )
if( info==0_${ik}$ ) then
! solve the system a*x = b, overwriting b with x.
call stdlib${ii}$_${ci}$pbtrs( uplo, n, kd, nrhs, ab, ldab, b, ldb, info )
end if
return
end subroutine stdlib${ii}$_${ci}$pbsv
#:endif
#:endfor
module subroutine stdlib${ii}$_spbsvx( fact, uplo, n, kd, nrhs, ab, ldab, afb, ldafb,equed, s, b, ldb, x, &
!! SPBSVX uses the Cholesky factorization A = U**T*U or A = L*L**T to
!! compute the solution to a real system of linear equations
!! A * X = B,
!! where A is an N-by-N symmetric positive definite band matrix and X
!! and B are N-by-NRHS matrices.
!! Error bounds on the solution and a condition estimate are also
!! provided.
ldx, rcond, ferr, berr,work, iwork, info )
! -- lapack driver 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(inout) :: equed
character, intent(in) :: fact, uplo
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: kd, ldab, ldafb, ldb, ldx, n, nrhs
real(sp), intent(out) :: rcond
! Array Arguments
integer(${ik}$), intent(out) :: iwork(*)
real(sp), intent(inout) :: ab(ldab,*), afb(ldafb,*), b(ldb,*), s(*)
real(sp), intent(out) :: berr(*), ferr(*), work(*), x(ldx,*)
! =====================================================================
! Local Scalars
logical(lk) :: equil, nofact, rcequ, upper
integer(${ik}$) :: i, infequ, j, j1, j2
real(sp) :: amax, anorm, bignum, scond, smax, smin, smlnum
! Intrinsic Functions
! Executable Statements
info = 0_${ik}$
nofact = stdlib_lsame( fact, 'N' )
equil = stdlib_lsame( fact, 'E' )
upper = stdlib_lsame( uplo, 'U' )
if( nofact .or. equil ) then
equed = 'N'
rcequ = .false.
else
rcequ = stdlib_lsame( equed, 'Y' )
smlnum = stdlib${ii}$_slamch( 'SAFE MINIMUM' )
bignum = one / smlnum
end if
! test the input parameters.
if( .not.nofact .and. .not.equil .and. .not.stdlib_lsame( fact, 'F' ) )then
info = -1_${ik}$
else if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then
info = -2_${ik}$
else if( n<0_${ik}$ ) then
info = -3_${ik}$
else if( kd<0_${ik}$ ) then
info = -4_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -5_${ik}$
else if( ldab<kd+1 ) then
info = -7_${ik}$
else if( ldafb<kd+1 ) then
info = -9_${ik}$
else if( stdlib_lsame( fact, 'F' ) .and. .not.( rcequ .or. stdlib_lsame( equed, 'N' ) )&
) then
info = -10_${ik}$
else
if( rcequ ) then
smin = bignum
smax = zero
do j = 1, n
smin = min( smin, s( j ) )
smax = max( smax, s( j ) )
end do
if( smin<=zero ) then
info = -11_${ik}$
else if( n>0_${ik}$ ) then
scond = max( smin, smlnum ) / min( smax, bignum )
else
scond = one
end if
end if
if( info==0_${ik}$ ) then
if( ldb<max( 1_${ik}$, n ) ) then
info = -13_${ik}$
else if( ldx<max( 1_${ik}$, n ) ) then
info = -15_${ik}$
end if
end if
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'SPBSVX', -info )
return
end if
if( equil ) then
! compute row and column scalings to equilibrate the matrix a.
call stdlib${ii}$_spbequ( uplo, n, kd, ab, ldab, s, scond, amax, infequ )
if( infequ==0_${ik}$ ) then
! equilibrate the matrix.
call stdlib${ii}$_slaqsb( uplo, n, kd, ab, ldab, s, scond, amax, equed )
rcequ = stdlib_lsame( equed, 'Y' )
end if
end if
! scale the right-hand side.
if( rcequ ) then
do j = 1, nrhs
do i = 1, n
b( i, j ) = s( i )*b( i, j )
end do
end do
end if
if( nofact .or. equil ) then
! compute the cholesky factorization a = u**t *u or a = l*l**t.
if( upper ) then
do j = 1, n
j1 = max( j-kd, 1_${ik}$ )
call stdlib${ii}$_scopy( j-j1+1, ab( kd+1-j+j1, j ), 1_${ik}$,afb( kd+1-j+j1, j ), 1_${ik}$ )
end do
else
do j = 1, n
j2 = min( j+kd, n )
call stdlib${ii}$_scopy( j2-j+1, ab( 1_${ik}$, j ), 1_${ik}$, afb( 1_${ik}$, j ), 1_${ik}$ )
end do
end if
call stdlib${ii}$_spbtrf( uplo, n, kd, afb, ldafb, info )
! return if info is non-zero.
if( info>0_${ik}$ )then
rcond = zero
return
end if
end if
! compute the norm of the matrix a.
anorm = stdlib${ii}$_slansb( '1', uplo, n, kd, ab, ldab, work )
! compute the reciprocal of the condition number of a.
call stdlib${ii}$_spbcon( uplo, n, kd, afb, ldafb, anorm, rcond, work, iwork,info )
! compute the solution matrix x.
call stdlib${ii}$_slacpy( 'FULL', n, nrhs, b, ldb, x, ldx )
call stdlib${ii}$_spbtrs( uplo, n, kd, nrhs, afb, ldafb, x, ldx, info )
! use iterative refinement to improve the computed solution and
! compute error bounds and backward error estimates for it.
call stdlib${ii}$_spbrfs( uplo, n, kd, nrhs, ab, ldab, afb, ldafb, b, ldb, x,ldx, ferr, berr,&
work, iwork, info )
! transform the solution matrix x to a solution of the original
! system.
if( rcequ ) then
do j = 1, nrhs
do i = 1, n
x( i, j ) = s( i )*x( i, j )
end do
end do
do j = 1, nrhs
ferr( j ) = ferr( j ) / scond
end do
end if
! set info = n+1 if the matrix is singular to working precision.
if( rcond<stdlib${ii}$_slamch( 'EPSILON' ) )info = n + 1_${ik}$
return
end subroutine stdlib${ii}$_spbsvx
module subroutine stdlib${ii}$_dpbsvx( fact, uplo, n, kd, nrhs, ab, ldab, afb, ldafb,equed, s, b, ldb, x, &
!! DPBSVX uses the Cholesky factorization A = U**T*U or A = L*L**T to
!! compute the solution to a real system of linear equations
!! A * X = B,
!! where A is an N-by-N symmetric positive definite band matrix and X
!! and B are N-by-NRHS matrices.
!! Error bounds on the solution and a condition estimate are also
!! provided.
ldx, rcond, ferr, berr,work, iwork, info )
! -- lapack driver 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(inout) :: equed
character, intent(in) :: fact, uplo
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: kd, ldab, ldafb, ldb, ldx, n, nrhs
real(dp), intent(out) :: rcond
! Array Arguments
integer(${ik}$), intent(out) :: iwork(*)
real(dp), intent(inout) :: ab(ldab,*), afb(ldafb,*), b(ldb,*), s(*)
real(dp), intent(out) :: berr(*), ferr(*), work(*), x(ldx,*)
! =====================================================================
! Local Scalars
logical(lk) :: equil, nofact, rcequ, upper
integer(${ik}$) :: i, infequ, j, j1, j2
real(dp) :: amax, anorm, bignum, scond, smax, smin, smlnum
! Intrinsic Functions
! Executable Statements
info = 0_${ik}$
nofact = stdlib_lsame( fact, 'N' )
equil = stdlib_lsame( fact, 'E' )
upper = stdlib_lsame( uplo, 'U' )
if( nofact .or. equil ) then
equed = 'N'
rcequ = .false.
else
rcequ = stdlib_lsame( equed, 'Y' )
smlnum = stdlib${ii}$_dlamch( 'SAFE MINIMUM' )
bignum = one / smlnum
end if
! test the input parameters.
if( .not.nofact .and. .not.equil .and. .not.stdlib_lsame( fact, 'F' ) )then
info = -1_${ik}$
else if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then
info = -2_${ik}$
else if( n<0_${ik}$ ) then
info = -3_${ik}$
else if( kd<0_${ik}$ ) then
info = -4_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -5_${ik}$
else if( ldab<kd+1 ) then
info = -7_${ik}$
else if( ldafb<kd+1 ) then
info = -9_${ik}$
else if( stdlib_lsame( fact, 'F' ) .and. .not.( rcequ .or. stdlib_lsame( equed, 'N' ) )&
) then
info = -10_${ik}$
else
if( rcequ ) then
smin = bignum
smax = zero
do j = 1, n
smin = min( smin, s( j ) )
smax = max( smax, s( j ) )
end do
if( smin<=zero ) then
info = -11_${ik}$
else if( n>0_${ik}$ ) then
scond = max( smin, smlnum ) / min( smax, bignum )
else
scond = one
end if
end if
if( info==0_${ik}$ ) then
if( ldb<max( 1_${ik}$, n ) ) then
info = -13_${ik}$
else if( ldx<max( 1_${ik}$, n ) ) then
info = -15_${ik}$
end if
end if
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'DPBSVX', -info )
return
end if
if( equil ) then
! compute row and column scalings to equilibrate the matrix a.
call stdlib${ii}$_dpbequ( uplo, n, kd, ab, ldab, s, scond, amax, infequ )
if( infequ==0_${ik}$ ) then
! equilibrate the matrix.
call stdlib${ii}$_dlaqsb( uplo, n, kd, ab, ldab, s, scond, amax, equed )
rcequ = stdlib_lsame( equed, 'Y' )
end if
end if
! scale the right-hand side.
if( rcequ ) then
do j = 1, nrhs
do i = 1, n
b( i, j ) = s( i )*b( i, j )
end do
end do
end if
if( nofact .or. equil ) then
! compute the cholesky factorization a = u**t *u or a = l*l**t.
if( upper ) then
do j = 1, n
j1 = max( j-kd, 1_${ik}$ )
call stdlib${ii}$_dcopy( j-j1+1, ab( kd+1-j+j1, j ), 1_${ik}$,afb( kd+1-j+j1, j ), 1_${ik}$ )
end do
else
do j = 1, n
j2 = min( j+kd, n )
call stdlib${ii}$_dcopy( j2-j+1, ab( 1_${ik}$, j ), 1_${ik}$, afb( 1_${ik}$, j ), 1_${ik}$ )
end do
end if
call stdlib${ii}$_dpbtrf( uplo, n, kd, afb, ldafb, info )
! return if info is non-zero.
if( info>0_${ik}$ )then
rcond = zero
return
end if
end if
! compute the norm of the matrix a.
anorm = stdlib${ii}$_dlansb( '1', uplo, n, kd, ab, ldab, work )
! compute the reciprocal of the condition number of a.
call stdlib${ii}$_dpbcon( uplo, n, kd, afb, ldafb, anorm, rcond, work, iwork,info )
! compute the solution matrix x.
call stdlib${ii}$_dlacpy( 'FULL', n, nrhs, b, ldb, x, ldx )
call stdlib${ii}$_dpbtrs( uplo, n, kd, nrhs, afb, ldafb, x, ldx, info )
! use iterative refinement to improve the computed solution and
! compute error bounds and backward error estimates for it.
call stdlib${ii}$_dpbrfs( uplo, n, kd, nrhs, ab, ldab, afb, ldafb, b, ldb, x,ldx, ferr, berr,&
work, iwork, info )
! transform the solution matrix x to a solution of the original
! system.
if( rcequ ) then
do j = 1, nrhs
do i = 1, n
x( i, j ) = s( i )*x( i, j )
end do
end do
do j = 1, nrhs
ferr( j ) = ferr( j ) / scond
end do
end if
! set info = n+1 if the matrix is singular to working precision.
if( rcond<stdlib${ii}$_dlamch( 'EPSILON' ) )info = n + 1_${ik}$
return
end subroutine stdlib${ii}$_dpbsvx
#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
module subroutine stdlib${ii}$_${ri}$pbsvx( fact, uplo, n, kd, nrhs, ab, ldab, afb, ldafb,equed, s, b, ldb, x, &
!! DPBSVX: uses the Cholesky factorization A = U**T*U or A = L*L**T to
!! compute the solution to a real system of linear equations
!! A * X = B,
!! where A is an N-by-N symmetric positive definite band matrix and X
!! and B are N-by-NRHS matrices.
!! Error bounds on the solution and a condition estimate are also
!! provided.
ldx, rcond, ferr, berr,work, iwork, info )
! -- lapack driver 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(inout) :: equed
character, intent(in) :: fact, uplo
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: kd, ldab, ldafb, ldb, ldx, n, nrhs
real(${rk}$), intent(out) :: rcond
! Array Arguments
integer(${ik}$), intent(out) :: iwork(*)
real(${rk}$), intent(inout) :: ab(ldab,*), afb(ldafb,*), b(ldb,*), s(*)
real(${rk}$), intent(out) :: berr(*), ferr(*), work(*), x(ldx,*)
! =====================================================================
! Local Scalars
logical(lk) :: equil, nofact, rcequ, upper
integer(${ik}$) :: i, infequ, j, j1, j2
real(${rk}$) :: amax, anorm, bignum, scond, smax, smin, smlnum
! Intrinsic Functions
! Executable Statements
info = 0_${ik}$
nofact = stdlib_lsame( fact, 'N' )
equil = stdlib_lsame( fact, 'E' )
upper = stdlib_lsame( uplo, 'U' )
if( nofact .or. equil ) then
equed = 'N'
rcequ = .false.
else
rcequ = stdlib_lsame( equed, 'Y' )
smlnum = stdlib${ii}$_${ri}$lamch( 'SAFE MINIMUM' )
bignum = one / smlnum
end if
! test the input parameters.
if( .not.nofact .and. .not.equil .and. .not.stdlib_lsame( fact, 'F' ) )then
info = -1_${ik}$
else if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then
info = -2_${ik}$
else if( n<0_${ik}$ ) then
info = -3_${ik}$
else if( kd<0_${ik}$ ) then
info = -4_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -5_${ik}$
else if( ldab<kd+1 ) then
info = -7_${ik}$
else if( ldafb<kd+1 ) then
info = -9_${ik}$
else if( stdlib_lsame( fact, 'F' ) .and. .not.( rcequ .or. stdlib_lsame( equed, 'N' ) )&
) then
info = -10_${ik}$
else
if( rcequ ) then
smin = bignum
smax = zero
do j = 1, n
smin = min( smin, s( j ) )
smax = max( smax, s( j ) )
end do
if( smin<=zero ) then
info = -11_${ik}$
else if( n>0_${ik}$ ) then
scond = max( smin, smlnum ) / min( smax, bignum )
else
scond = one
end if
end if
if( info==0_${ik}$ ) then
if( ldb<max( 1_${ik}$, n ) ) then
info = -13_${ik}$
else if( ldx<max( 1_${ik}$, n ) ) then
info = -15_${ik}$
end if
end if
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'DPBSVX', -info )
return
end if
if( equil ) then
! compute row and column scalings to equilibrate the matrix a.
call stdlib${ii}$_${ri}$pbequ( uplo, n, kd, ab, ldab, s, scond, amax, infequ )
if( infequ==0_${ik}$ ) then
! equilibrate the matrix.
call stdlib${ii}$_${ri}$laqsb( uplo, n, kd, ab, ldab, s, scond, amax, equed )
rcequ = stdlib_lsame( equed, 'Y' )
end if
end if
! scale the right-hand side.
if( rcequ ) then
do j = 1, nrhs
do i = 1, n
b( i, j ) = s( i )*b( i, j )
end do
end do
end if
if( nofact .or. equil ) then
! compute the cholesky factorization a = u**t *u or a = l*l**t.
if( upper ) then
do j = 1, n
j1 = max( j-kd, 1_${ik}$ )
call stdlib${ii}$_${ri}$copy( j-j1+1, ab( kd+1-j+j1, j ), 1_${ik}$,afb( kd+1-j+j1, j ), 1_${ik}$ )
end do
else
do j = 1, n
j2 = min( j+kd, n )
call stdlib${ii}$_${ri}$copy( j2-j+1, ab( 1_${ik}$, j ), 1_${ik}$, afb( 1_${ik}$, j ), 1_${ik}$ )
end do
end if
call stdlib${ii}$_${ri}$pbtrf( uplo, n, kd, afb, ldafb, info )
! return if info is non-zero.
if( info>0_${ik}$ )then
rcond = zero
return
end if
end if
! compute the norm of the matrix a.
anorm = stdlib${ii}$_${ri}$lansb( '1', uplo, n, kd, ab, ldab, work )
! compute the reciprocal of the condition number of a.
call stdlib${ii}$_${ri}$pbcon( uplo, n, kd, afb, ldafb, anorm, rcond, work, iwork,info )
! compute the solution matrix x.
call stdlib${ii}$_${ri}$lacpy( 'FULL', n, nrhs, b, ldb, x, ldx )
call stdlib${ii}$_${ri}$pbtrs( uplo, n, kd, nrhs, afb, ldafb, x, ldx, info )
! use iterative refinement to improve the computed solution and
! compute error bounds and backward error estimates for it.
call stdlib${ii}$_${ri}$pbrfs( uplo, n, kd, nrhs, ab, ldab, afb, ldafb, b, ldb, x,ldx, ferr, berr,&
work, iwork, info )
! transform the solution matrix x to a solution of the original
! system.
if( rcequ ) then
do j = 1, nrhs
do i = 1, n
x( i, j ) = s( i )*x( i, j )
end do
end do
do j = 1, nrhs
ferr( j ) = ferr( j ) / scond
end do
end if
! set info = n+1 if the matrix is singular to working precision.
if( rcond<stdlib${ii}$_${ri}$lamch( 'EPSILON' ) )info = n + 1_${ik}$
return
end subroutine stdlib${ii}$_${ri}$pbsvx
#:endif
#:endfor
module subroutine stdlib${ii}$_cpbsvx( fact, uplo, n, kd, nrhs, ab, ldab, afb, ldafb,equed, s, b, ldb, x, &
!! CPBSVX uses the Cholesky factorization A = U**H*U or A = L*L**H to
!! compute the solution to a complex system of linear equations
!! A * X = B,
!! where A is an N-by-N Hermitian positive definite band matrix and X
!! and B are N-by-NRHS matrices.
!! Error bounds on the solution and a condition estimate are also
!! provided.
ldx, rcond, ferr, berr,work, rwork, info )
! -- lapack driver 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(inout) :: equed
character, intent(in) :: fact, uplo
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: kd, ldab, ldafb, ldb, ldx, n, nrhs
real(sp), intent(out) :: rcond
! Array Arguments
real(sp), intent(out) :: berr(*), ferr(*), rwork(*)
real(sp), intent(inout) :: s(*)
complex(sp), intent(inout) :: ab(ldab,*), afb(ldafb,*), b(ldb,*)
complex(sp), intent(out) :: work(*), x(ldx,*)
! =====================================================================
! Local Scalars
logical(lk) :: equil, nofact, rcequ, upper
integer(${ik}$) :: i, infequ, j, j1, j2
real(sp) :: amax, anorm, bignum, scond, smax, smin, smlnum
! Intrinsic Functions
! Executable Statements
info = 0_${ik}$
nofact = stdlib_lsame( fact, 'N' )
equil = stdlib_lsame( fact, 'E' )
upper = stdlib_lsame( uplo, 'U' )
if( nofact .or. equil ) then
equed = 'N'
rcequ = .false.
else
rcequ = stdlib_lsame( equed, 'Y' )
smlnum = stdlib${ii}$_slamch( 'SAFE MINIMUM' )
bignum = one / smlnum
end if
! test the input parameters.
if( .not.nofact .and. .not.equil .and. .not.stdlib_lsame( fact, 'F' ) )then
info = -1_${ik}$
else if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then
info = -2_${ik}$
else if( n<0_${ik}$ ) then
info = -3_${ik}$
else if( kd<0_${ik}$ ) then
info = -4_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -5_${ik}$
else if( ldab<kd+1 ) then
info = -7_${ik}$
else if( ldafb<kd+1 ) then
info = -9_${ik}$
else if( stdlib_lsame( fact, 'F' ) .and. .not.( rcequ .or. stdlib_lsame( equed, 'N' ) )&
) then
info = -10_${ik}$
else
if( rcequ ) then
smin = bignum
smax = zero
do j = 1, n
smin = min( smin, s( j ) )
smax = max( smax, s( j ) )
end do
if( smin<=zero ) then
info = -11_${ik}$
else if( n>0_${ik}$ ) then
scond = max( smin, smlnum ) / min( smax, bignum )
else
scond = one
end if
end if
if( info==0_${ik}$ ) then
if( ldb<max( 1_${ik}$, n ) ) then
info = -13_${ik}$
else if( ldx<max( 1_${ik}$, n ) ) then
info = -15_${ik}$
end if
end if
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'CPBSVX', -info )
return
end if
if( equil ) then
! compute row and column scalings to equilibrate the matrix a.
call stdlib${ii}$_cpbequ( uplo, n, kd, ab, ldab, s, scond, amax, infequ )
if( infequ==0_${ik}$ ) then
! equilibrate the matrix.
call stdlib${ii}$_claqhb( uplo, n, kd, ab, ldab, s, scond, amax, equed )
rcequ = stdlib_lsame( equed, 'Y' )
end if
end if
! scale the right-hand side.
if( rcequ ) then
do j = 1, nrhs
do i = 1, n
b( i, j ) = s( i )*b( i, j )
end do
end do
end if
if( nofact .or. equil ) then
! compute the cholesky factorization a = u**h *u or a = l*l**h.
if( upper ) then
do j = 1, n
j1 = max( j-kd, 1_${ik}$ )
call stdlib${ii}$_ccopy( j-j1+1, ab( kd+1-j+j1, j ), 1_${ik}$,afb( kd+1-j+j1, j ), 1_${ik}$ )
end do
else
do j = 1, n
j2 = min( j+kd, n )
call stdlib${ii}$_ccopy( j2-j+1, ab( 1_${ik}$, j ), 1_${ik}$, afb( 1_${ik}$, j ), 1_${ik}$ )
end do
end if
call stdlib${ii}$_cpbtrf( uplo, n, kd, afb, ldafb, info )
! return if info is non-zero.
if( info>0_${ik}$ )then
rcond = zero
return
end if
end if
! compute the norm of the matrix a.
anorm = stdlib${ii}$_clanhb( '1', uplo, n, kd, ab, ldab, rwork )
! compute the reciprocal of the condition number of a.
call stdlib${ii}$_cpbcon( uplo, n, kd, afb, ldafb, anorm, rcond, work, rwork,info )
! compute the solution matrix x.
call stdlib${ii}$_clacpy( 'FULL', n, nrhs, b, ldb, x, ldx )
call stdlib${ii}$_cpbtrs( uplo, n, kd, nrhs, afb, ldafb, x, ldx, info )
! use iterative refinement to improve the computed solution and
! compute error bounds and backward error estimates for it.
call stdlib${ii}$_cpbrfs( uplo, n, kd, nrhs, ab, ldab, afb, ldafb, b, ldb, x,ldx, ferr, berr,&
work, rwork, info )
! transform the solution matrix x to a solution of the original
! system.
if( rcequ ) then
do j = 1, nrhs
do i = 1, n
x( i, j ) = s( i )*x( i, j )
end do
end do
do j = 1, nrhs
ferr( j ) = ferr( j ) / scond
end do
end if
! set info = n+1 if the matrix is singular to working precision.
if( rcond<stdlib${ii}$_slamch( 'EPSILON' ) )info = n + 1_${ik}$
return
end subroutine stdlib${ii}$_cpbsvx
module subroutine stdlib${ii}$_zpbsvx( fact, uplo, n, kd, nrhs, ab, ldab, afb, ldafb,equed, s, b, ldb, x, &
!! ZPBSVX uses the Cholesky factorization A = U**H*U or A = L*L**H to
!! compute the solution to a complex system of linear equations
!! A * X = B,
!! where A is an N-by-N Hermitian positive definite band matrix and X
!! and B are N-by-NRHS matrices.
!! Error bounds on the solution and a condition estimate are also
!! provided.
ldx, rcond, ferr, berr,work, rwork, info )
! -- lapack driver 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(inout) :: equed
character, intent(in) :: fact, uplo
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: kd, ldab, ldafb, ldb, ldx, n, nrhs
real(dp), intent(out) :: rcond
! Array Arguments
real(dp), intent(out) :: berr(*), ferr(*), rwork(*)
real(dp), intent(inout) :: s(*)
complex(dp), intent(inout) :: ab(ldab,*), afb(ldafb,*), b(ldb,*)
complex(dp), intent(out) :: work(*), x(ldx,*)
! =====================================================================
! Local Scalars
logical(lk) :: equil, nofact, rcequ, upper
integer(${ik}$) :: i, infequ, j, j1, j2
real(dp) :: amax, anorm, bignum, scond, smax, smin, smlnum
! Intrinsic Functions
! Executable Statements
info = 0_${ik}$
nofact = stdlib_lsame( fact, 'N' )
equil = stdlib_lsame( fact, 'E' )
upper = stdlib_lsame( uplo, 'U' )
if( nofact .or. equil ) then
equed = 'N'
rcequ = .false.
else
rcequ = stdlib_lsame( equed, 'Y' )
smlnum = stdlib${ii}$_dlamch( 'SAFE MINIMUM' )
bignum = one / smlnum
end if
! test the input parameters.
if( .not.nofact .and. .not.equil .and. .not.stdlib_lsame( fact, 'F' ) )then
info = -1_${ik}$
else if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then
info = -2_${ik}$
else if( n<0_${ik}$ ) then
info = -3_${ik}$
else if( kd<0_${ik}$ ) then
info = -4_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -5_${ik}$
else if( ldab<kd+1 ) then
info = -7_${ik}$
else if( ldafb<kd+1 ) then
info = -9_${ik}$
else if( stdlib_lsame( fact, 'F' ) .and. .not.( rcequ .or. stdlib_lsame( equed, 'N' ) )&
) then
info = -10_${ik}$
else
if( rcequ ) then
smin = bignum
smax = zero
do j = 1, n
smin = min( smin, s( j ) )
smax = max( smax, s( j ) )
end do
if( smin<=zero ) then
info = -11_${ik}$
else if( n>0_${ik}$ ) then
scond = max( smin, smlnum ) / min( smax, bignum )
else
scond = one
end if
end if
if( info==0_${ik}$ ) then
if( ldb<max( 1_${ik}$, n ) ) then
info = -13_${ik}$
else if( ldx<max( 1_${ik}$, n ) ) then
info = -15_${ik}$
end if
end if
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'ZPBSVX', -info )
return
end if
if( equil ) then
! compute row and column scalings to equilibrate the matrix a.
call stdlib${ii}$_zpbequ( uplo, n, kd, ab, ldab, s, scond, amax, infequ )
if( infequ==0_${ik}$ ) then
! equilibrate the matrix.
call stdlib${ii}$_zlaqhb( uplo, n, kd, ab, ldab, s, scond, amax, equed )
rcequ = stdlib_lsame( equed, 'Y' )
end if
end if
! scale the right-hand side.
if( rcequ ) then
do j = 1, nrhs
do i = 1, n
b( i, j ) = s( i )*b( i, j )
end do
end do
end if
if( nofact .or. equil ) then
! compute the cholesky factorization a = u**h *u or a = l*l**h.
if( upper ) then
do j = 1, n
j1 = max( j-kd, 1_${ik}$ )
call stdlib${ii}$_zcopy( j-j1+1, ab( kd+1-j+j1, j ), 1_${ik}$,afb( kd+1-j+j1, j ), 1_${ik}$ )
end do
else
do j = 1, n
j2 = min( j+kd, n )
call stdlib${ii}$_zcopy( j2-j+1, ab( 1_${ik}$, j ), 1_${ik}$, afb( 1_${ik}$, j ), 1_${ik}$ )
end do
end if
call stdlib${ii}$_zpbtrf( uplo, n, kd, afb, ldafb, info )
! return if info is non-zero.
if( info>0_${ik}$ )then
rcond = zero
return
end if
end if
! compute the norm of the matrix a.
anorm = stdlib${ii}$_zlanhb( '1', uplo, n, kd, ab, ldab, rwork )
! compute the reciprocal of the condition number of a.
call stdlib${ii}$_zpbcon( uplo, n, kd, afb, ldafb, anorm, rcond, work, rwork,info )
! compute the solution matrix x.
call stdlib${ii}$_zlacpy( 'FULL', n, nrhs, b, ldb, x, ldx )
call stdlib${ii}$_zpbtrs( uplo, n, kd, nrhs, afb, ldafb, x, ldx, info )
! use iterative refinement to improve the computed solution and
! compute error bounds and backward error estimates for it.
call stdlib${ii}$_zpbrfs( uplo, n, kd, nrhs, ab, ldab, afb, ldafb, b, ldb, x,ldx, ferr, berr,&
work, rwork, info )
! transform the solution matrix x to a solution of the original
! system.
if( rcequ ) then
do j = 1, nrhs
do i = 1, n
x( i, j ) = s( i )*x( i, j )
end do
end do
do j = 1, nrhs
ferr( j ) = ferr( j ) / scond
end do
end if
! set info = n+1 if the matrix is singular to working precision.
if( rcond<stdlib${ii}$_dlamch( 'EPSILON' ) )info = n + 1_${ik}$
return
end subroutine stdlib${ii}$_zpbsvx
#:for ck,ct,ci in CMPLX_KINDS_TYPES
#:if not ck in ["sp","dp"]
module subroutine stdlib${ii}$_${ci}$pbsvx( fact, uplo, n, kd, nrhs, ab, ldab, afb, ldafb,equed, s, b, ldb, x, &
!! ZPBSVX: uses the Cholesky factorization A = U**H*U or A = L*L**H to
!! compute the solution to a complex system of linear equations
!! A * X = B,
!! where A is an N-by-N Hermitian positive definite band matrix and X
!! and B are N-by-NRHS matrices.
!! Error bounds on the solution and a condition estimate are also
!! provided.
ldx, rcond, ferr, berr,work, rwork, info )
! -- lapack driver 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(inout) :: equed
character, intent(in) :: fact, uplo
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: kd, ldab, ldafb, ldb, ldx, n, nrhs
real(${ck}$), intent(out) :: rcond
! Array Arguments
real(${ck}$), intent(out) :: berr(*), ferr(*), rwork(*)
real(${ck}$), intent(inout) :: s(*)
complex(${ck}$), intent(inout) :: ab(ldab,*), afb(ldafb,*), b(ldb,*)
complex(${ck}$), intent(out) :: work(*), x(ldx,*)
! =====================================================================
! Local Scalars
logical(lk) :: equil, nofact, rcequ, upper
integer(${ik}$) :: i, infequ, j, j1, j2
real(${ck}$) :: amax, anorm, bignum, scond, smax, smin, smlnum
! Intrinsic Functions
! Executable Statements
info = 0_${ik}$
nofact = stdlib_lsame( fact, 'N' )
equil = stdlib_lsame( fact, 'E' )
upper = stdlib_lsame( uplo, 'U' )
if( nofact .or. equil ) then
equed = 'N'
rcequ = .false.
else
rcequ = stdlib_lsame( equed, 'Y' )
smlnum = stdlib${ii}$_${c2ri(ci)}$lamch( 'SAFE MINIMUM' )
bignum = one / smlnum
end if
! test the input parameters.
if( .not.nofact .and. .not.equil .and. .not.stdlib_lsame( fact, 'F' ) )then
info = -1_${ik}$
else if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then
info = -2_${ik}$
else if( n<0_${ik}$ ) then
info = -3_${ik}$
else if( kd<0_${ik}$ ) then
info = -4_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -5_${ik}$
else if( ldab<kd+1 ) then
info = -7_${ik}$
else if( ldafb<kd+1 ) then
info = -9_${ik}$
else if( stdlib_lsame( fact, 'F' ) .and. .not.( rcequ .or. stdlib_lsame( equed, 'N' ) )&
) then
info = -10_${ik}$
else
if( rcequ ) then
smin = bignum
smax = zero
do j = 1, n
smin = min( smin, s( j ) )
smax = max( smax, s( j ) )
end do
if( smin<=zero ) then
info = -11_${ik}$
else if( n>0_${ik}$ ) then
scond = max( smin, smlnum ) / min( smax, bignum )
else
scond = one
end if
end if
if( info==0_${ik}$ ) then
if( ldb<max( 1_${ik}$, n ) ) then
info = -13_${ik}$
else if( ldx<max( 1_${ik}$, n ) ) then
info = -15_${ik}$
end if
end if
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'ZPBSVX', -info )
return
end if
if( equil ) then
! compute row and column scalings to equilibrate the matrix a.
call stdlib${ii}$_${ci}$pbequ( uplo, n, kd, ab, ldab, s, scond, amax, infequ )
if( infequ==0_${ik}$ ) then
! equilibrate the matrix.
call stdlib${ii}$_${ci}$laqhb( uplo, n, kd, ab, ldab, s, scond, amax, equed )
rcequ = stdlib_lsame( equed, 'Y' )
end if
end if
! scale the right-hand side.
if( rcequ ) then
do j = 1, nrhs
do i = 1, n
b( i, j ) = s( i )*b( i, j )
end do
end do
end if
if( nofact .or. equil ) then
! compute the cholesky factorization a = u**h *u or a = l*l**h.
if( upper ) then
do j = 1, n
j1 = max( j-kd, 1_${ik}$ )
call stdlib${ii}$_${ci}$copy( j-j1+1, ab( kd+1-j+j1, j ), 1_${ik}$,afb( kd+1-j+j1, j ), 1_${ik}$ )
end do
else
do j = 1, n
j2 = min( j+kd, n )
call stdlib${ii}$_${ci}$copy( j2-j+1, ab( 1_${ik}$, j ), 1_${ik}$, afb( 1_${ik}$, j ), 1_${ik}$ )
end do
end if
call stdlib${ii}$_${ci}$pbtrf( uplo, n, kd, afb, ldafb, info )
! return if info is non-zero.
if( info>0_${ik}$ )then
rcond = zero
return
end if
end if
! compute the norm of the matrix a.
anorm = stdlib${ii}$_${ci}$lanhb( '1', uplo, n, kd, ab, ldab, rwork )
! compute the reciprocal of the condition number of a.
call stdlib${ii}$_${ci}$pbcon( uplo, n, kd, afb, ldafb, anorm, rcond, work, rwork,info )
! compute the solution matrix x.
call stdlib${ii}$_${ci}$lacpy( 'FULL', n, nrhs, b, ldb, x, ldx )
call stdlib${ii}$_${ci}$pbtrs( uplo, n, kd, nrhs, afb, ldafb, x, ldx, info )
! use iterative refinement to improve the computed solution and
! compute error bounds and backward error estimates for it.
call stdlib${ii}$_${ci}$pbrfs( uplo, n, kd, nrhs, ab, ldab, afb, ldafb, b, ldb, x,ldx, ferr, berr,&
work, rwork, info )
! transform the solution matrix x to a solution of the original
! system.
if( rcequ ) then
do j = 1, nrhs
do i = 1, n
x( i, j ) = s( i )*x( i, j )
end do
end do
do j = 1, nrhs
ferr( j ) = ferr( j ) / scond
end do
end if
! set info = n+1 if the matrix is singular to working precision.
if( rcond<stdlib${ii}$_${c2ri(ci)}$lamch( 'EPSILON' ) )info = n + 1_${ik}$
return
end subroutine stdlib${ii}$_${ci}$pbsvx
#:endif
#:endfor
pure module subroutine stdlib${ii}$_sptsv( n, nrhs, d, e, b, ldb, info )
!! SPTSV computes the solution to a real system of linear equations
!! A*X = B, where A is an N-by-N symmetric positive definite tridiagonal
!! matrix, and X and B are N-by-NRHS matrices.
!! A is factored as A = L*D*L**T, and the factored form of A is then
!! used to solve the system of equations.
! -- lapack driver 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
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: ldb, n, nrhs
! Array Arguments
real(sp), intent(inout) :: b(ldb,*), d(*), e(*)
! =====================================================================
! Intrinsic Functions
! Executable Statements
! test the input parameters.
info = 0_${ik}$
if( n<0_${ik}$ ) then
info = -1_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -2_${ik}$
else if( ldb<max( 1_${ik}$, n ) ) then
info = -6_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'SPTSV ', -info )
return
end if
! compute the l*d*l**t (or u**t*d*u) factorization of a.
call stdlib${ii}$_spttrf( n, d, e, info )
if( info==0_${ik}$ ) then
! solve the system a*x = b, overwriting b with x.
call stdlib${ii}$_spttrs( n, nrhs, d, e, b, ldb, info )
end if
return
end subroutine stdlib${ii}$_sptsv
pure module subroutine stdlib${ii}$_dptsv( n, nrhs, d, e, b, ldb, info )
!! DPTSV computes the solution to a real system of linear equations
!! A*X = B, where A is an N-by-N symmetric positive definite tridiagonal
!! matrix, and X and B are N-by-NRHS matrices.
!! A is factored as A = L*D*L**T, and the factored form of A is then
!! used to solve the system of equations.
! -- lapack driver 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
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: ldb, n, nrhs
! Array Arguments
real(dp), intent(inout) :: b(ldb,*), d(*), e(*)
! =====================================================================
! Intrinsic Functions
! Executable Statements
! test the input parameters.
info = 0_${ik}$
if( n<0_${ik}$ ) then
info = -1_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -2_${ik}$
else if( ldb<max( 1_${ik}$, n ) ) then
info = -6_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'DPTSV ', -info )
return
end if
! compute the l*d*l**t (or u**t*d*u) factorization of a.
call stdlib${ii}$_dpttrf( n, d, e, info )
if( info==0_${ik}$ ) then
! solve the system a*x = b, overwriting b with x.
call stdlib${ii}$_dpttrs( n, nrhs, d, e, b, ldb, info )
end if
return
end subroutine stdlib${ii}$_dptsv
#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
pure module subroutine stdlib${ii}$_${ri}$ptsv( n, nrhs, d, e, b, ldb, info )
!! DPTSV: computes the solution to a real system of linear equations
!! A*X = B, where A is an N-by-N symmetric positive definite tridiagonal
!! matrix, and X and B are N-by-NRHS matrices.
!! A is factored as A = L*D*L**T, and the factored form of A is then
!! used to solve the system of equations.
! -- lapack driver 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
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: ldb, n, nrhs
! Array Arguments
real(${rk}$), intent(inout) :: b(ldb,*), d(*), e(*)
! =====================================================================
! Intrinsic Functions
! Executable Statements
! test the input parameters.
info = 0_${ik}$
if( n<0_${ik}$ ) then
info = -1_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -2_${ik}$
else if( ldb<max( 1_${ik}$, n ) ) then
info = -6_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'DPTSV ', -info )
return
end if
! compute the l*d*l**t (or u**t*d*u) factorization of a.
call stdlib${ii}$_${ri}$pttrf( n, d, e, info )
if( info==0_${ik}$ ) then
! solve the system a*x = b, overwriting b with x.
call stdlib${ii}$_${ri}$pttrs( n, nrhs, d, e, b, ldb, info )
end if
return
end subroutine stdlib${ii}$_${ri}$ptsv
#:endif
#:endfor
pure module subroutine stdlib${ii}$_cptsv( n, nrhs, d, e, b, ldb, info )
!! CPTSV computes the solution to a complex system of linear equations
!! A*X = B, where A is an N-by-N Hermitian positive definite tridiagonal
!! matrix, and X and B are N-by-NRHS matrices.
!! A is factored as A = L*D*L**H, and the factored form of A is then
!! used to solve the system of equations.
! -- lapack driver 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
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: ldb, n, nrhs
! Array Arguments
real(sp), intent(inout) :: d(*)
complex(sp), intent(inout) :: b(ldb,*), e(*)
! =====================================================================
! Intrinsic Functions
! Executable Statements
! test the input parameters.
info = 0_${ik}$
if( n<0_${ik}$ ) then
info = -1_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -2_${ik}$
else if( ldb<max( 1_${ik}$, n ) ) then
info = -6_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'CPTSV ', -info )
return
end if
! compute the l*d*l**h (or u**h*d*u) factorization of a.
call stdlib${ii}$_cpttrf( n, d, e, info )
if( info==0_${ik}$ ) then
! solve the system a*x = b, overwriting b with x.
call stdlib${ii}$_cpttrs( 'LOWER', n, nrhs, d, e, b, ldb, info )
end if
return
end subroutine stdlib${ii}$_cptsv
pure module subroutine stdlib${ii}$_zptsv( n, nrhs, d, e, b, ldb, info )
!! ZPTSV computes the solution to a complex system of linear equations
!! A*X = B, where A is an N-by-N Hermitian positive definite tridiagonal
!! matrix, and X and B are N-by-NRHS matrices.
!! A is factored as A = L*D*L**H, and the factored form of A is then
!! used to solve the system of equations.
! -- lapack driver 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
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: ldb, n, nrhs
! Array Arguments
real(dp), intent(inout) :: d(*)
complex(dp), intent(inout) :: b(ldb,*), e(*)
! =====================================================================
! Intrinsic Functions
! Executable Statements
! test the input parameters.
info = 0_${ik}$
if( n<0_${ik}$ ) then
info = -1_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -2_${ik}$
else if( ldb<max( 1_${ik}$, n ) ) then
info = -6_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'ZPTSV ', -info )
return
end if
! compute the l*d*l**h (or u**h*d*u) factorization of a.
call stdlib${ii}$_zpttrf( n, d, e, info )
if( info==0_${ik}$ ) then
! solve the system a*x = b, overwriting b with x.
call stdlib${ii}$_zpttrs( 'LOWER', n, nrhs, d, e, b, ldb, info )
end if
return
end subroutine stdlib${ii}$_zptsv
#:for ck,ct,ci in CMPLX_KINDS_TYPES
#:if not ck in ["sp","dp"]
pure module subroutine stdlib${ii}$_${ci}$ptsv( n, nrhs, d, e, b, ldb, info )
!! ZPTSV: computes the solution to a complex system of linear equations
!! A*X = B, where A is an N-by-N Hermitian positive definite tridiagonal
!! matrix, and X and B are N-by-NRHS matrices.
!! A is factored as A = L*D*L**H, and the factored form of A is then
!! used to solve the system of equations.
! -- lapack driver 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
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: ldb, n, nrhs
! Array Arguments
real(${ck}$), intent(inout) :: d(*)
complex(${ck}$), intent(inout) :: b(ldb,*), e(*)
! =====================================================================
! Intrinsic Functions
! Executable Statements
! test the input parameters.
info = 0_${ik}$
if( n<0_${ik}$ ) then
info = -1_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -2_${ik}$
else if( ldb<max( 1_${ik}$, n ) ) then
info = -6_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'ZPTSV ', -info )
return
end if
! compute the l*d*l**h (or u**h*d*u) factorization of a.
call stdlib${ii}$_${ci}$pttrf( n, d, e, info )
if( info==0_${ik}$ ) then
! solve the system a*x = b, overwriting b with x.
call stdlib${ii}$_${ci}$pttrs( 'LOWER', n, nrhs, d, e, b, ldb, info )
end if
return
end subroutine stdlib${ii}$_${ci}$ptsv
#:endif
#:endfor
pure module subroutine stdlib${ii}$_sptsvx( fact, n, nrhs, d, e, df, ef, b, ldb, x, ldx,rcond, ferr, berr,&
!! SPTSVX uses the factorization A = L*D*L**T to compute the solution
!! to a real system of linear equations A*X = B, where A is an N-by-N
!! symmetric positive definite tridiagonal matrix and X and B are
!! N-by-NRHS matrices.
!! Error bounds on the solution and a condition estimate are also
!! provided.
work, info )
! -- lapack driver 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) :: fact
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: ldb, ldx, n, nrhs
real(sp), intent(out) :: rcond
! Array Arguments
real(sp), intent(in) :: b(ldb,*), d(*), e(*)
real(sp), intent(out) :: berr(*), ferr(*), work(*), x(ldx,*)
real(sp), intent(inout) :: df(*), ef(*)
! =====================================================================
! Local Scalars
logical(lk) :: nofact
real(sp) :: anorm
! Intrinsic Functions
! Executable Statements
! test the input parameters.
info = 0_${ik}$
nofact = stdlib_lsame( fact, 'N' )
if( .not.nofact .and. .not.stdlib_lsame( fact, 'F' ) ) 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 = -9_${ik}$
else if( ldx<max( 1_${ik}$, n ) ) then
info = -11_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'SPTSVX', -info )
return
end if
if( nofact ) then
! compute the l*d*l**t (or u**t*d*u) factorization of a.
call stdlib${ii}$_scopy( n, d, 1_${ik}$, df, 1_${ik}$ )
if( n>1_${ik}$ )call stdlib${ii}$_scopy( n-1, e, 1_${ik}$, ef, 1_${ik}$ )
call stdlib${ii}$_spttrf( n, df, ef, info )
! return if info is non-zero.
if( info>0_${ik}$ )then
rcond = zero
return
end if
end if
! compute the norm of the matrix a.
anorm = stdlib${ii}$_slanst( '1', n, d, e )
! compute the reciprocal of the condition number of a.
call stdlib${ii}$_sptcon( n, df, ef, anorm, rcond, work, info )
! compute the solution vectors x.
call stdlib${ii}$_slacpy( 'FULL', n, nrhs, b, ldb, x, ldx )
call stdlib${ii}$_spttrs( n, nrhs, df, ef, x, ldx, info )
! use iterative refinement to improve the computed solutions and
! compute error bounds and backward error estimates for them.
call stdlib${ii}$_sptrfs( n, nrhs, d, e, df, ef, b, ldb, x, ldx, ferr, berr,work, info )
! set info = n+1 if the matrix is singular to working precision.
if( rcond<stdlib${ii}$_slamch( 'EPSILON' ) )info = n + 1_${ik}$
return
end subroutine stdlib${ii}$_sptsvx
pure module subroutine stdlib${ii}$_dptsvx( fact, n, nrhs, d, e, df, ef, b, ldb, x, ldx,rcond, ferr, berr,&
!! DPTSVX uses the factorization A = L*D*L**T to compute the solution
!! to a real system of linear equations A*X = B, where A is an N-by-N
!! symmetric positive definite tridiagonal matrix and X and B are
!! N-by-NRHS matrices.
!! Error bounds on the solution and a condition estimate are also
!! provided.
work, info )
! -- lapack driver 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) :: fact
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: ldb, ldx, n, nrhs
real(dp), intent(out) :: rcond
! Array Arguments
real(dp), intent(in) :: b(ldb,*), d(*), e(*)
real(dp), intent(out) :: berr(*), ferr(*), work(*), x(ldx,*)
real(dp), intent(inout) :: df(*), ef(*)
! =====================================================================
! Local Scalars
logical(lk) :: nofact
real(dp) :: anorm
! Intrinsic Functions
! Executable Statements
! test the input parameters.
info = 0_${ik}$
nofact = stdlib_lsame( fact, 'N' )
if( .not.nofact .and. .not.stdlib_lsame( fact, 'F' ) ) 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 = -9_${ik}$
else if( ldx<max( 1_${ik}$, n ) ) then
info = -11_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'DPTSVX', -info )
return
end if
if( nofact ) then
! compute the l*d*l**t (or u**t*d*u) factorization of a.
call stdlib${ii}$_dcopy( n, d, 1_${ik}$, df, 1_${ik}$ )
if( n>1_${ik}$ )call stdlib${ii}$_dcopy( n-1, e, 1_${ik}$, ef, 1_${ik}$ )
call stdlib${ii}$_dpttrf( n, df, ef, info )
! return if info is non-zero.
if( info>0_${ik}$ )then
rcond = zero
return
end if
end if
! compute the norm of the matrix a.
anorm = stdlib${ii}$_dlanst( '1', n, d, e )
! compute the reciprocal of the condition number of a.
call stdlib${ii}$_dptcon( n, df, ef, anorm, rcond, work, info )
! compute the solution vectors x.
call stdlib${ii}$_dlacpy( 'FULL', n, nrhs, b, ldb, x, ldx )
call stdlib${ii}$_dpttrs( n, nrhs, df, ef, x, ldx, info )
! use iterative refinement to improve the computed solutions and
! compute error bounds and backward error estimates for them.
call stdlib${ii}$_dptrfs( n, nrhs, d, e, df, ef, b, ldb, x, ldx, ferr, berr,work, info )
! set info = n+1 if the matrix is singular to working precision.
if( rcond<stdlib${ii}$_dlamch( 'EPSILON' ) )info = n + 1_${ik}$
return
end subroutine stdlib${ii}$_dptsvx
#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
pure module subroutine stdlib${ii}$_${ri}$ptsvx( fact, n, nrhs, d, e, df, ef, b, ldb, x, ldx,rcond, ferr, berr,&
!! DPTSVX: uses the factorization A = L*D*L**T to compute the solution
!! to a real system of linear equations A*X = B, where A is an N-by-N
!! symmetric positive definite tridiagonal matrix and X and B are
!! N-by-NRHS matrices.
!! Error bounds on the solution and a condition estimate are also
!! provided.
work, info )
! -- lapack driver 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) :: fact
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: ldb, ldx, n, nrhs
real(${rk}$), intent(out) :: rcond
! Array Arguments
real(${rk}$), intent(in) :: b(ldb,*), d(*), e(*)
real(${rk}$), intent(out) :: berr(*), ferr(*), work(*), x(ldx,*)
real(${rk}$), intent(inout) :: df(*), ef(*)
! =====================================================================
! Local Scalars
logical(lk) :: nofact
real(${rk}$) :: anorm
! Intrinsic Functions
! Executable Statements
! test the input parameters.
info = 0_${ik}$
nofact = stdlib_lsame( fact, 'N' )
if( .not.nofact .and. .not.stdlib_lsame( fact, 'F' ) ) 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 = -9_${ik}$
else if( ldx<max( 1_${ik}$, n ) ) then
info = -11_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'DPTSVX', -info )
return
end if
if( nofact ) then
! compute the l*d*l**t (or u**t*d*u) factorization of a.
call stdlib${ii}$_${ri}$copy( n, d, 1_${ik}$, df, 1_${ik}$ )
if( n>1_${ik}$ )call stdlib${ii}$_${ri}$copy( n-1, e, 1_${ik}$, ef, 1_${ik}$ )
call stdlib${ii}$_${ri}$pttrf( n, df, ef, info )
! return if info is non-zero.
if( info>0_${ik}$ )then
rcond = zero
return
end if
end if
! compute the norm of the matrix a.
anorm = stdlib${ii}$_${ri}$lanst( '1', n, d, e )
! compute the reciprocal of the condition number of a.
call stdlib${ii}$_${ri}$ptcon( n, df, ef, anorm, rcond, work, info )
! compute the solution vectors x.
call stdlib${ii}$_${ri}$lacpy( 'FULL', n, nrhs, b, ldb, x, ldx )
call stdlib${ii}$_${ri}$pttrs( n, nrhs, df, ef, x, ldx, info )
! use iterative refinement to improve the computed solutions and
! compute error bounds and backward error estimates for them.
call stdlib${ii}$_${ri}$ptrfs( n, nrhs, d, e, df, ef, b, ldb, x, ldx, ferr, berr,work, info )
! set info = n+1 if the matrix is singular to working precision.
if( rcond<stdlib${ii}$_${ri}$lamch( 'EPSILON' ) )info = n + 1_${ik}$
return
end subroutine stdlib${ii}$_${ri}$ptsvx
#:endif
#:endfor
pure module subroutine stdlib${ii}$_cptsvx( fact, n, nrhs, d, e, df, ef, b, ldb, x, ldx,rcond, ferr, berr,&
!! CPTSVX uses the factorization A = L*D*L**H to compute the solution
!! to a complex system of linear equations A*X = B, where A is an
!! N-by-N Hermitian positive definite tridiagonal matrix and X and B
!! are N-by-NRHS matrices.
!! Error bounds on the solution and a condition estimate are also
!! provided.
work, rwork, info )
! -- lapack driver 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) :: fact
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: ldb, ldx, n, nrhs
real(sp), intent(out) :: rcond
! Array Arguments
real(sp), intent(out) :: berr(*), ferr(*), rwork(*)
real(sp), intent(in) :: d(*)
real(sp), intent(inout) :: df(*)
complex(sp), intent(in) :: b(ldb,*), e(*)
complex(sp), intent(inout) :: ef(*)
complex(sp), intent(out) :: work(*), x(ldx,*)
! =====================================================================
! Local Scalars
logical(lk) :: nofact
real(sp) :: anorm
! Intrinsic Functions
! Executable Statements
! test the input parameters.
info = 0_${ik}$
nofact = stdlib_lsame( fact, 'N' )
if( .not.nofact .and. .not.stdlib_lsame( fact, 'F' ) ) 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 = -9_${ik}$
else if( ldx<max( 1_${ik}$, n ) ) then
info = -11_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'CPTSVX', -info )
return
end if
if( nofact ) then
! compute the l*d*l**h (or u**h*d*u) factorization of a.
call stdlib${ii}$_scopy( n, d, 1_${ik}$, df, 1_${ik}$ )
if( n>1_${ik}$ )call stdlib${ii}$_ccopy( n-1, e, 1_${ik}$, ef, 1_${ik}$ )
call stdlib${ii}$_cpttrf( n, df, ef, info )
! return if info is non-zero.
if( info>0_${ik}$ )then
rcond = zero
return
end if
end if
! compute the norm of the matrix a.
anorm = stdlib${ii}$_clanht( '1', n, d, e )
! compute the reciprocal of the condition number of a.
call stdlib${ii}$_cptcon( n, df, ef, anorm, rcond, rwork, info )
! compute the solution vectors x.
call stdlib${ii}$_clacpy( 'FULL', n, nrhs, b, ldb, x, ldx )
call stdlib${ii}$_cpttrs( 'LOWER', n, nrhs, df, ef, x, ldx, info )
! use iterative refinement to improve the computed solutions and
! compute error bounds and backward error estimates for them.
call stdlib${ii}$_cptrfs( 'LOWER', n, nrhs, d, e, df, ef, b, ldb, x, ldx, ferr,berr, work, &
rwork, info )
! set info = n+1 if the matrix is singular to working precision.
if( rcond<stdlib${ii}$_slamch( 'EPSILON' ) )info = n + 1_${ik}$
return
end subroutine stdlib${ii}$_cptsvx
pure module subroutine stdlib${ii}$_zptsvx( fact, n, nrhs, d, e, df, ef, b, ldb, x, ldx,rcond, ferr, berr,&
!! ZPTSVX uses the factorization A = L*D*L**H to compute the solution
!! to a complex system of linear equations A*X = B, where A is an
!! N-by-N Hermitian positive definite tridiagonal matrix and X and B
!! are N-by-NRHS matrices.
!! Error bounds on the solution and a condition estimate are also
!! provided.
work, rwork, info )
! -- lapack driver 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) :: fact
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: ldb, ldx, n, nrhs
real(dp), intent(out) :: rcond
! Array Arguments
real(dp), intent(out) :: berr(*), ferr(*), rwork(*)
real(dp), intent(in) :: d(*)
real(dp), intent(inout) :: df(*)
complex(dp), intent(in) :: b(ldb,*), e(*)
complex(dp), intent(inout) :: ef(*)
complex(dp), intent(out) :: work(*), x(ldx,*)
! =====================================================================
! Local Scalars
logical(lk) :: nofact
real(dp) :: anorm
! Intrinsic Functions
! Executable Statements
! test the input parameters.
info = 0_${ik}$
nofact = stdlib_lsame( fact, 'N' )
if( .not.nofact .and. .not.stdlib_lsame( fact, 'F' ) ) 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 = -9_${ik}$
else if( ldx<max( 1_${ik}$, n ) ) then
info = -11_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'ZPTSVX', -info )
return
end if
if( nofact ) then
! compute the l*d*l**h (or u**h*d*u) factorization of a.
call stdlib${ii}$_dcopy( n, d, 1_${ik}$, df, 1_${ik}$ )
if( n>1_${ik}$ )call stdlib${ii}$_zcopy( n-1, e, 1_${ik}$, ef, 1_${ik}$ )
call stdlib${ii}$_zpttrf( n, df, ef, info )
! return if info is non-zero.
if( info>0_${ik}$ )then
rcond = zero
return
end if
end if
! compute the norm of the matrix a.
anorm = stdlib${ii}$_zlanht( '1', n, d, e )
! compute the reciprocal of the condition number of a.
call stdlib${ii}$_zptcon( n, df, ef, anorm, rcond, rwork, info )
! compute the solution vectors x.
call stdlib${ii}$_zlacpy( 'FULL', n, nrhs, b, ldb, x, ldx )
call stdlib${ii}$_zpttrs( 'LOWER', n, nrhs, df, ef, x, ldx, info )
! use iterative refinement to improve the computed solutions and
! compute error bounds and backward error estimates for them.
call stdlib${ii}$_zptrfs( 'LOWER', n, nrhs, d, e, df, ef, b, ldb, x, ldx, ferr,berr, work, &
rwork, info )
! set info = n+1 if the matrix is singular to working precision.
if( rcond<stdlib${ii}$_dlamch( 'EPSILON' ) )info = n + 1_${ik}$
return
end subroutine stdlib${ii}$_zptsvx
#:for ck,ct,ci in CMPLX_KINDS_TYPES
#:if not ck in ["sp","dp"]
pure module subroutine stdlib${ii}$_${ci}$ptsvx( fact, n, nrhs, d, e, df, ef, b, ldb, x, ldx,rcond, ferr, berr,&
!! ZPTSVX: uses the factorization A = L*D*L**H to compute the solution
!! to a complex system of linear equations A*X = B, where A is an
!! N-by-N Hermitian positive definite tridiagonal matrix and X and B
!! are N-by-NRHS matrices.
!! Error bounds on the solution and a condition estimate are also
!! provided.
work, rwork, info )
! -- lapack driver 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) :: fact
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: ldb, ldx, n, nrhs
real(${ck}$), intent(out) :: rcond
! Array Arguments
real(${ck}$), intent(out) :: berr(*), ferr(*), rwork(*)
real(${ck}$), intent(in) :: d(*)
real(${ck}$), intent(inout) :: df(*)
complex(${ck}$), intent(in) :: b(ldb,*), e(*)
complex(${ck}$), intent(inout) :: ef(*)
complex(${ck}$), intent(out) :: work(*), x(ldx,*)
! =====================================================================
! Local Scalars
logical(lk) :: nofact
real(${ck}$) :: anorm
! Intrinsic Functions
! Executable Statements
! test the input parameters.
info = 0_${ik}$
nofact = stdlib_lsame( fact, 'N' )
if( .not.nofact .and. .not.stdlib_lsame( fact, 'F' ) ) 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 = -9_${ik}$
else if( ldx<max( 1_${ik}$, n ) ) then
info = -11_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'ZPTSVX', -info )
return
end if
if( nofact ) then
! compute the l*d*l**h (or u**h*d*u) factorization of a.
call stdlib${ii}$_${c2ri(ci)}$copy( n, d, 1_${ik}$, df, 1_${ik}$ )
if( n>1_${ik}$ )call stdlib${ii}$_${ci}$copy( n-1, e, 1_${ik}$, ef, 1_${ik}$ )
call stdlib${ii}$_${ci}$pttrf( n, df, ef, info )
! return if info is non-zero.
if( info>0_${ik}$ )then
rcond = zero
return
end if
end if
! compute the norm of the matrix a.
anorm = stdlib${ii}$_${ci}$lanht( '1', n, d, e )
! compute the reciprocal of the condition number of a.
call stdlib${ii}$_${ci}$ptcon( n, df, ef, anorm, rcond, rwork, info )
! compute the solution vectors x.
call stdlib${ii}$_${ci}$lacpy( 'FULL', n, nrhs, b, ldb, x, ldx )
call stdlib${ii}$_${ci}$pttrs( 'LOWER', n, nrhs, df, ef, x, ldx, info )
! use iterative refinement to improve the computed solutions and
! compute error bounds and backward error estimates for them.
call stdlib${ii}$_${ci}$ptrfs( 'LOWER', n, nrhs, d, e, df, ef, b, ldb, x, ldx, ferr,berr, work, &
rwork, info )
! set info = n+1 if the matrix is singular to working precision.
if( rcond<stdlib${ii}$_${c2ri(ci)}$lamch( 'EPSILON' ) )info = n + 1_${ik}$
return
end subroutine stdlib${ii}$_${ci}$ptsvx
#:endif
#:endfor
#:endfor
end submodule stdlib_lapack_solve_chol
