#:include "common.fypp"
submodule(stdlib_lapack_solve) stdlib_lapack_solve_lu
implicit none
contains
#:for ik,it,ii in LINALG_INT_KINDS_TYPES
pure module subroutine stdlib${ii}$_sgesv( n, nrhs, a, lda, ipiv, b, ldb, info )
!! SGESV computes the solution to a real system of linear equations
!! A * X = B,
!! where A is an N-by-N matrix and X and B are N-by-NRHS matrices.
!! The LU decomposition with partial pivoting and row interchanges is
!! used to factor A as
!! A = P * L * U,
!! where P is a permutation matrix, L is unit lower triangular, and U is
!! upper triangular. 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
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: lda, ldb, n, nrhs
! Array Arguments
integer(${ik}$), intent(out) :: ipiv(*)
real(sp), intent(inout) :: a(lda,*), b(ldb,*)
! =====================================================================
! 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( lda<max( 1_${ik}$, n ) ) then
info = -4_${ik}$
else if( ldb<max( 1_${ik}$, n ) ) then
info = -7_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'SGESV ', -info )
return
end if
! compute the lu factorization of a.
call stdlib${ii}$_sgetrf( n, n, a, lda, ipiv, info )
if( info==0_${ik}$ ) then
! solve the system a*x = b, overwriting b with x.
call stdlib${ii}$_sgetrs( 'NO TRANSPOSE', n, nrhs, a, lda, ipiv, b, ldb,info )
end if
return
end subroutine stdlib${ii}$_sgesv
pure module subroutine stdlib${ii}$_dgesv( n, nrhs, a, lda, ipiv, b, ldb, info )
!! DGESV computes the solution to a real system of linear equations
!! A * X = B,
!! where A is an N-by-N matrix and X and B are N-by-NRHS matrices.
!! The LU decomposition with partial pivoting and row interchanges is
!! used to factor A as
!! A = P * L * U,
!! where P is a permutation matrix, L is unit lower triangular, and U is
!! upper triangular. 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
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: lda, ldb, n, nrhs
! Array Arguments
integer(${ik}$), intent(out) :: ipiv(*)
real(dp), intent(inout) :: a(lda,*), b(ldb,*)
! =====================================================================
! 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( lda<max( 1_${ik}$, n ) ) then
info = -4_${ik}$
else if( ldb<max( 1_${ik}$, n ) ) then
info = -7_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'DGESV ', -info )
return
end if
! compute the lu factorization of a.
call stdlib${ii}$_dgetrf( n, n, a, lda, ipiv, info )
if( info==0_${ik}$ ) then
! solve the system a*x = b, overwriting b with x.
call stdlib${ii}$_dgetrs( 'NO TRANSPOSE', n, nrhs, a, lda, ipiv, b, ldb,info )
end if
return
end subroutine stdlib${ii}$_dgesv
#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
pure module subroutine stdlib${ii}$_${ri}$gesv( n, nrhs, a, lda, ipiv, b, ldb, info )
!! DGESV: computes the solution to a real system of linear equations
!! A * X = B,
!! where A is an N-by-N matrix and X and B are N-by-NRHS matrices.
!! The LU decomposition with partial pivoting and row interchanges is
!! used to factor A as
!! A = P * L * U,
!! where P is a permutation matrix, L is unit lower triangular, and U is
!! upper triangular. 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
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: lda, ldb, n, nrhs
! Array Arguments
integer(${ik}$), intent(out) :: ipiv(*)
real(${rk}$), intent(inout) :: a(lda,*), b(ldb,*)
! =====================================================================
! 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( lda<max( 1_${ik}$, n ) ) then
info = -4_${ik}$
else if( ldb<max( 1_${ik}$, n ) ) then
info = -7_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'DGESV ', -info )
return
end if
! compute the lu factorization of a.
call stdlib${ii}$_${ri}$getrf( n, n, a, lda, ipiv, info )
if( info==0_${ik}$ ) then
! solve the system a*x = b, overwriting b with x.
call stdlib${ii}$_${ri}$getrs( 'NO TRANSPOSE', n, nrhs, a, lda, ipiv, b, ldb,info )
end if
return
end subroutine stdlib${ii}$_${ri}$gesv
#:endif
#:endfor
pure module subroutine stdlib${ii}$_cgesv( n, nrhs, a, lda, ipiv, b, ldb, info )
!! CGESV computes the solution to a complex system of linear equations
!! A * X = B,
!! where A is an N-by-N matrix and X and B are N-by-NRHS matrices.
!! The LU decomposition with partial pivoting and row interchanges is
!! used to factor A as
!! A = P * L * U,
!! where P is a permutation matrix, L is unit lower triangular, and U is
!! upper triangular. 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
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: lda, ldb, n, nrhs
! Array Arguments
integer(${ik}$), intent(out) :: ipiv(*)
complex(sp), intent(inout) :: a(lda,*), b(ldb,*)
! =====================================================================
! 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( lda<max( 1_${ik}$, n ) ) then
info = -4_${ik}$
else if( ldb<max( 1_${ik}$, n ) ) then
info = -7_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'CGESV ', -info )
return
end if
! compute the lu factorization of a.
call stdlib${ii}$_cgetrf( n, n, a, lda, ipiv, info )
if( info==0_${ik}$ ) then
! solve the system a*x = b, overwriting b with x.
call stdlib${ii}$_cgetrs( 'NO TRANSPOSE', n, nrhs, a, lda, ipiv, b, ldb,info )
end if
return
end subroutine stdlib${ii}$_cgesv
pure module subroutine stdlib${ii}$_zgesv( n, nrhs, a, lda, ipiv, b, ldb, info )
!! ZGESV computes the solution to a complex system of linear equations
!! A * X = B,
!! where A is an N-by-N matrix and X and B are N-by-NRHS matrices.
!! The LU decomposition with partial pivoting and row interchanges is
!! used to factor A as
!! A = P * L * U,
!! where P is a permutation matrix, L is unit lower triangular, and U is
!! upper triangular. 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
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: lda, ldb, n, nrhs
! Array Arguments
integer(${ik}$), intent(out) :: ipiv(*)
complex(dp), intent(inout) :: a(lda,*), b(ldb,*)
! =====================================================================
! 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( lda<max( 1_${ik}$, n ) ) then
info = -4_${ik}$
else if( ldb<max( 1_${ik}$, n ) ) then
info = -7_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'ZGESV ', -info )
return
end if
! compute the lu factorization of a.
call stdlib${ii}$_zgetrf( n, n, a, lda, ipiv, info )
if( info==0_${ik}$ ) then
! solve the system a*x = b, overwriting b with x.
call stdlib${ii}$_zgetrs( 'NO TRANSPOSE', n, nrhs, a, lda, ipiv, b, ldb,info )
end if
return
end subroutine stdlib${ii}$_zgesv
#:for ck,ct,ci in CMPLX_KINDS_TYPES
#:if not ck in ["sp","dp"]
pure module subroutine stdlib${ii}$_${ci}$gesv( n, nrhs, a, lda, ipiv, b, ldb, info )
!! ZGESV: computes the solution to a complex system of linear equations
!! A * X = B,
!! where A is an N-by-N matrix and X and B are N-by-NRHS matrices.
!! The LU decomposition with partial pivoting and row interchanges is
!! used to factor A as
!! A = P * L * U,
!! where P is a permutation matrix, L is unit lower triangular, and U is
!! upper triangular. 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
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: lda, ldb, n, nrhs
! Array Arguments
integer(${ik}$), intent(out) :: ipiv(*)
complex(${ck}$), intent(inout) :: a(lda,*), b(ldb,*)
! =====================================================================
! 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( lda<max( 1_${ik}$, n ) ) then
info = -4_${ik}$
else if( ldb<max( 1_${ik}$, n ) ) then
info = -7_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'ZGESV ', -info )
return
end if
! compute the lu factorization of a.
call stdlib${ii}$_${ci}$getrf( n, n, a, lda, ipiv, info )
if( info==0_${ik}$ ) then
! solve the system a*x = b, overwriting b with x.
call stdlib${ii}$_${ci}$getrs( 'NO TRANSPOSE', n, nrhs, a, lda, ipiv, b, ldb,info )
end if
return
end subroutine stdlib${ii}$_${ci}$gesv
#:endif
#:endfor
module subroutine stdlib${ii}$_sgesvx( fact, trans, n, nrhs, a, lda, af, ldaf, ipiv,equed, r, c, b, ldb, &
!! SGESVX uses the LU factorization to compute the solution to a real
!! system of linear equations
!! A * X = B,
!! where A is an N-by-N matrix and X and B are N-by-NRHS matrices.
!! Error bounds on the solution and a condition estimate are also
!! provided.
x, 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, trans
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(inout) :: ipiv(*)
integer(${ik}$), intent(out) :: iwork(*)
real(sp), intent(inout) :: a(lda,*), af(ldaf,*), b(ldb,*), c(*), r(*)
real(sp), intent(out) :: berr(*), ferr(*), work(*), x(ldx,*)
! =====================================================================
! Local Scalars
logical(lk) :: colequ, equil, nofact, notran, rowequ
character :: norm
integer(${ik}$) :: i, infequ, j
real(sp) :: amax, anorm, bignum, colcnd, rcmax, rcmin, rowcnd, rpvgrw, smlnum
! Intrinsic Functions
! Executable Statements
info = 0_${ik}$
nofact = stdlib_lsame( fact, 'N' )
equil = stdlib_lsame( fact, 'E' )
notran = stdlib_lsame( trans, 'N' )
if( nofact .or. equil ) then
equed = 'N'
rowequ = .false.
colequ = .false.
else
rowequ = stdlib_lsame( equed, 'R' ) .or. stdlib_lsame( equed, 'B' )
colequ = stdlib_lsame( equed, 'C' ) .or. stdlib_lsame( equed, 'B' )
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.notran .and. .not.stdlib_lsame( trans, 'T' ) .and. .not.stdlib_lsame( &
trans, 'C' ) ) 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.( rowequ .or. colequ .or. stdlib_lsame( &
equed, 'N' ) ) ) then
info = -10_${ik}$
else
if( rowequ ) then
rcmin = bignum
rcmax = zero
do j = 1, n
rcmin = min( rcmin, r( j ) )
rcmax = max( rcmax, r( j ) )
end do
if( rcmin<=zero ) then
info = -11_${ik}$
else if( n>0_${ik}$ ) then
rowcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
else
rowcnd = one
end if
end if
if( colequ .and. info==0_${ik}$ ) then
rcmin = bignum
rcmax = zero
do j = 1, n
rcmin = min( rcmin, c( j ) )
rcmax = max( rcmax, c( j ) )
end do
if( rcmin<=zero ) then
info = -12_${ik}$
else if( n>0_${ik}$ ) then
colcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
else
colcnd = one
end if
end if
if( info==0_${ik}$ ) then
if( ldb<max( 1_${ik}$, n ) ) then
info = -14_${ik}$
else if( ldx<max( 1_${ik}$, n ) ) then
info = -16_${ik}$
end if
end if
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'SGESVX', -info )
return
end if
if( equil ) then
! compute row and column scalings to equilibrate the matrix a.
call stdlib${ii}$_sgeequ( n, n, a, lda, r, c, rowcnd, colcnd, amax, infequ )
if( infequ==0_${ik}$ ) then
! equilibrate the matrix.
call stdlib${ii}$_slaqge( n, n, a, lda, r, c, rowcnd, colcnd, amax,equed )
rowequ = stdlib_lsame( equed, 'R' ) .or. stdlib_lsame( equed, 'B' )
colequ = stdlib_lsame( equed, 'C' ) .or. stdlib_lsame( equed, 'B' )
end if
end if
! scale the right hand side.
if( notran ) then
if( rowequ ) then
do j = 1, nrhs
do i = 1, n
b( i, j ) = r( i )*b( i, j )
end do
end do
end if
else if( colequ ) then
do j = 1, nrhs
do i = 1, n
b( i, j ) = c( i )*b( i, j )
end do
end do
end if
if( nofact .or. equil ) then
! compute the lu factorization of a.
call stdlib${ii}$_slacpy( 'FULL', n, n, a, lda, af, ldaf )
call stdlib${ii}$_sgetrf( n, n, af, ldaf, ipiv, info )
! return if info is non-zero.
if( info>0_${ik}$ ) then
! compute the reciprocal pivot growth factor of the
! leading rank-deficient info columns of a.
rpvgrw = stdlib${ii}$_slantr( 'M', 'U', 'N', info, info, af, ldaf,work )
if( rpvgrw==zero ) then
rpvgrw = one
else
rpvgrw = stdlib${ii}$_slange( 'M', n, info, a, lda, work ) / rpvgrw
end if
work( 1_${ik}$ ) = rpvgrw
rcond = zero
return
end if
end if
! compute the norm of the matrix a and the
! reciprocal pivot growth factor rpvgrw.
if( notran ) then
norm = '1'
else
norm = 'I'
end if
anorm = stdlib${ii}$_slange( norm, n, n, a, lda, work )
rpvgrw = stdlib${ii}$_slantr( 'M', 'U', 'N', n, n, af, ldaf, work )
if( rpvgrw==zero ) then
rpvgrw = one
else
rpvgrw = stdlib${ii}$_slange( 'M', n, n, a, lda, work ) / rpvgrw
end if
! compute the reciprocal of the condition number of a.
call stdlib${ii}$_sgecon( norm, 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}$_sgetrs( trans, n, nrhs, af, ldaf, ipiv, x, ldx, info )
! use iterative refinement to improve the computed solution and
! compute error bounds and backward error estimates for it.
call stdlib${ii}$_sgerfs( trans, n, nrhs, a, lda, af, ldaf, ipiv, b, ldb, x,ldx, ferr, berr, &
work, iwork, info )
! transform the solution matrix x to a solution of the original
! system.
if( notran ) then
if( colequ ) then
do j = 1, nrhs
do i = 1, n
x( i, j ) = c( i )*x( i, j )
end do
end do
do j = 1, nrhs
ferr( j ) = ferr( j ) / colcnd
end do
end if
else if( rowequ ) then
do j = 1, nrhs
do i = 1, n
x( i, j ) = r( i )*x( i, j )
end do
end do
do j = 1, nrhs
ferr( j ) = ferr( j ) / rowcnd
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}$
work( 1_${ik}$ ) = rpvgrw
return
end subroutine stdlib${ii}$_sgesvx
module subroutine stdlib${ii}$_dgesvx( fact, trans, n, nrhs, a, lda, af, ldaf, ipiv,equed, r, c, b, ldb, &
!! DGESVX uses the LU factorization to compute the solution to a real
!! system of linear equations
!! A * X = B,
!! where A is an N-by-N matrix and X and B are N-by-NRHS matrices.
!! Error bounds on the solution and a condition estimate are also
!! provided.
x, 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, trans
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(inout) :: ipiv(*)
integer(${ik}$), intent(out) :: iwork(*)
real(dp), intent(inout) :: a(lda,*), af(ldaf,*), b(ldb,*), c(*), r(*)
real(dp), intent(out) :: berr(*), ferr(*), work(*), x(ldx,*)
! =====================================================================
! Local Scalars
logical(lk) :: colequ, equil, nofact, notran, rowequ
character :: norm
integer(${ik}$) :: i, infequ, j
real(dp) :: amax, anorm, bignum, colcnd, rcmax, rcmin, rowcnd, rpvgrw, smlnum
! Intrinsic Functions
! Executable Statements
info = 0_${ik}$
nofact = stdlib_lsame( fact, 'N' )
equil = stdlib_lsame( fact, 'E' )
notran = stdlib_lsame( trans, 'N' )
if( nofact .or. equil ) then
equed = 'N'
rowequ = .false.
colequ = .false.
else
rowequ = stdlib_lsame( equed, 'R' ) .or. stdlib_lsame( equed, 'B' )
colequ = stdlib_lsame( equed, 'C' ) .or. stdlib_lsame( equed, 'B' )
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.notran .and. .not.stdlib_lsame( trans, 'T' ) .and. .not.stdlib_lsame( &
trans, 'C' ) ) 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.( rowequ .or. colequ .or. stdlib_lsame( &
equed, 'N' ) ) ) then
info = -10_${ik}$
else
if( rowequ ) then
rcmin = bignum
rcmax = zero
do j = 1, n
rcmin = min( rcmin, r( j ) )
rcmax = max( rcmax, r( j ) )
end do
if( rcmin<=zero ) then
info = -11_${ik}$
else if( n>0_${ik}$ ) then
rowcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
else
rowcnd = one
end if
end if
if( colequ .and. info==0_${ik}$ ) then
rcmin = bignum
rcmax = zero
do j = 1, n
rcmin = min( rcmin, c( j ) )
rcmax = max( rcmax, c( j ) )
end do
if( rcmin<=zero ) then
info = -12_${ik}$
else if( n>0_${ik}$ ) then
colcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
else
colcnd = one
end if
end if
if( info==0_${ik}$ ) then
if( ldb<max( 1_${ik}$, n ) ) then
info = -14_${ik}$
else if( ldx<max( 1_${ik}$, n ) ) then
info = -16_${ik}$
end if
end if
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'DGESVX', -info )
return
end if
if( equil ) then
! compute row and column scalings to equilibrate the matrix a.
call stdlib${ii}$_dgeequ( n, n, a, lda, r, c, rowcnd, colcnd, amax, infequ )
if( infequ==0_${ik}$ ) then
! equilibrate the matrix.
call stdlib${ii}$_dlaqge( n, n, a, lda, r, c, rowcnd, colcnd, amax,equed )
rowequ = stdlib_lsame( equed, 'R' ) .or. stdlib_lsame( equed, 'B' )
colequ = stdlib_lsame( equed, 'C' ) .or. stdlib_lsame( equed, 'B' )
end if
end if
! scale the right hand side.
if( notran ) then
if( rowequ ) then
do j = 1, nrhs
do i = 1, n
b( i, j ) = r( i )*b( i, j )
end do
end do
end if
else if( colequ ) then
do j = 1, nrhs
do i = 1, n
b( i, j ) = c( i )*b( i, j )
end do
end do
end if
if( nofact .or. equil ) then
! compute the lu factorization of a.
call stdlib${ii}$_dlacpy( 'FULL', n, n, a, lda, af, ldaf )
call stdlib${ii}$_dgetrf( n, n, af, ldaf, ipiv, info )
! return if info is non-zero.
if( info>0_${ik}$ ) then
! compute the reciprocal pivot growth factor of the
! leading rank-deficient info columns of a.
rpvgrw = stdlib${ii}$_dlantr( 'M', 'U', 'N', info, info, af, ldaf,work )
if( rpvgrw==zero ) then
rpvgrw = one
else
rpvgrw = stdlib${ii}$_dlange( 'M', n, info, a, lda, work ) / rpvgrw
end if
work( 1_${ik}$ ) = rpvgrw
rcond = zero
return
end if
end if
! compute the norm of the matrix a and the
! reciprocal pivot growth factor rpvgrw.
if( notran ) then
norm = '1'
else
norm = 'I'
end if
anorm = stdlib${ii}$_dlange( norm, n, n, a, lda, work )
rpvgrw = stdlib${ii}$_dlantr( 'M', 'U', 'N', n, n, af, ldaf, work )
if( rpvgrw==zero ) then
rpvgrw = one
else
rpvgrw = stdlib${ii}$_dlange( 'M', n, n, a, lda, work ) / rpvgrw
end if
! compute the reciprocal of the condition number of a.
call stdlib${ii}$_dgecon( norm, 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}$_dgetrs( trans, n, nrhs, af, ldaf, ipiv, x, ldx, info )
! use iterative refinement to improve the computed solution and
! compute error bounds and backward error estimates for it.
call stdlib${ii}$_dgerfs( trans, n, nrhs, a, lda, af, ldaf, ipiv, b, ldb, x,ldx, ferr, berr, &
work, iwork, info )
! transform the solution matrix x to a solution of the original
! system.
if( notran ) then
if( colequ ) then
do j = 1, nrhs
do i = 1, n
x( i, j ) = c( i )*x( i, j )
end do
end do
do j = 1, nrhs
ferr( j ) = ferr( j ) / colcnd
end do
end if
else if( rowequ ) then
do j = 1, nrhs
do i = 1, n
x( i, j ) = r( i )*x( i, j )
end do
end do
do j = 1, nrhs
ferr( j ) = ferr( j ) / rowcnd
end do
end if
work( 1_${ik}$ ) = rpvgrw
! 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}$_dgesvx
#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
module subroutine stdlib${ii}$_${ri}$gesvx( fact, trans, n, nrhs, a, lda, af, ldaf, ipiv,equed, r, c, b, ldb, &
!! DGESVX: uses the LU factorization to compute the solution to a real
!! system of linear equations
!! A * X = B,
!! where A is an N-by-N matrix and X and B are N-by-NRHS matrices.
!! Error bounds on the solution and a condition estimate are also
!! provided.
x, 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, trans
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(inout) :: ipiv(*)
integer(${ik}$), intent(out) :: iwork(*)
real(${rk}$), intent(inout) :: a(lda,*), af(ldaf,*), b(ldb,*), c(*), r(*)
real(${rk}$), intent(out) :: berr(*), ferr(*), work(*), x(ldx,*)
! =====================================================================
! Local Scalars
logical(lk) :: colequ, equil, nofact, notran, rowequ
character :: norm
integer(${ik}$) :: i, infequ, j
real(${rk}$) :: amax, anorm, bignum, colcnd, rcmax, rcmin, rowcnd, rpvgrw, smlnum
! Intrinsic Functions
! Executable Statements
info = 0_${ik}$
nofact = stdlib_lsame( fact, 'N' )
equil = stdlib_lsame( fact, 'E' )
notran = stdlib_lsame( trans, 'N' )
if( nofact .or. equil ) then
equed = 'N'
rowequ = .false.
colequ = .false.
else
rowequ = stdlib_lsame( equed, 'R' ) .or. stdlib_lsame( equed, 'B' )
colequ = stdlib_lsame( equed, 'C' ) .or. stdlib_lsame( equed, 'B' )
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.notran .and. .not.stdlib_lsame( trans, 'T' ) .and. .not.stdlib_lsame( &
trans, 'C' ) ) 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.( rowequ .or. colequ .or. stdlib_lsame( &
equed, 'N' ) ) ) then
info = -10_${ik}$
else
if( rowequ ) then
rcmin = bignum
rcmax = zero
do j = 1, n
rcmin = min( rcmin, r( j ) )
rcmax = max( rcmax, r( j ) )
end do
if( rcmin<=zero ) then
info = -11_${ik}$
else if( n>0_${ik}$ ) then
rowcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
else
rowcnd = one
end if
end if
if( colequ .and. info==0_${ik}$ ) then
rcmin = bignum
rcmax = zero
do j = 1, n
rcmin = min( rcmin, c( j ) )
rcmax = max( rcmax, c( j ) )
end do
if( rcmin<=zero ) then
info = -12_${ik}$
else if( n>0_${ik}$ ) then
colcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
else
colcnd = one
end if
end if
if( info==0_${ik}$ ) then
if( ldb<max( 1_${ik}$, n ) ) then
info = -14_${ik}$
else if( ldx<max( 1_${ik}$, n ) ) then
info = -16_${ik}$
end if
end if
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'DGESVX', -info )
return
end if
if( equil ) then
! compute row and column scalings to equilibrate the matrix a.
call stdlib${ii}$_${ri}$geequ( n, n, a, lda, r, c, rowcnd, colcnd, amax, infequ )
if( infequ==0_${ik}$ ) then
! equilibrate the matrix.
call stdlib${ii}$_${ri}$laqge( n, n, a, lda, r, c, rowcnd, colcnd, amax,equed )
rowequ = stdlib_lsame( equed, 'R' ) .or. stdlib_lsame( equed, 'B' )
colequ = stdlib_lsame( equed, 'C' ) .or. stdlib_lsame( equed, 'B' )
end if
end if
! scale the right hand side.
if( notran ) then
if( rowequ ) then
do j = 1, nrhs
do i = 1, n
b( i, j ) = r( i )*b( i, j )
end do
end do
end if
else if( colequ ) then
do j = 1, nrhs
do i = 1, n
b( i, j ) = c( i )*b( i, j )
end do
end do
end if
if( nofact .or. equil ) then
! compute the lu factorization of a.
call stdlib${ii}$_${ri}$lacpy( 'FULL', n, n, a, lda, af, ldaf )
call stdlib${ii}$_${ri}$getrf( n, n, af, ldaf, ipiv, info )
! return if info is non-zero.
if( info>0_${ik}$ ) then
! compute the reciprocal pivot growth factor of the
! leading rank-deficient info columns of a.
rpvgrw = stdlib${ii}$_${ri}$lantr( 'M', 'U', 'N', info, info, af, ldaf,work )
if( rpvgrw==zero ) then
rpvgrw = one
else
rpvgrw = stdlib${ii}$_${ri}$lange( 'M', n, info, a, lda, work ) / rpvgrw
end if
work( 1_${ik}$ ) = rpvgrw
rcond = zero
return
end if
end if
! compute the norm of the matrix a and the
! reciprocal pivot growth factor rpvgrw.
if( notran ) then
norm = '1'
else
norm = 'I'
end if
anorm = stdlib${ii}$_${ri}$lange( norm, n, n, a, lda, work )
rpvgrw = stdlib${ii}$_${ri}$lantr( 'M', 'U', 'N', n, n, af, ldaf, work )
if( rpvgrw==zero ) then
rpvgrw = one
else
rpvgrw = stdlib${ii}$_${ri}$lange( 'M', n, n, a, lda, work ) / rpvgrw
end if
! compute the reciprocal of the condition number of a.
call stdlib${ii}$_${ri}$gecon( norm, 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}$getrs( trans, n, nrhs, af, ldaf, ipiv, 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}$gerfs( trans, n, nrhs, a, lda, af, ldaf, ipiv, b, ldb, x,ldx, ferr, berr, &
work, iwork, info )
! transform the solution matrix x to a solution of the original
! system.
if( notran ) then
if( colequ ) then
do j = 1, nrhs
do i = 1, n
x( i, j ) = c( i )*x( i, j )
end do
end do
do j = 1, nrhs
ferr( j ) = ferr( j ) / colcnd
end do
end if
else if( rowequ ) then
do j = 1, nrhs
do i = 1, n
x( i, j ) = r( i )*x( i, j )
end do
end do
do j = 1, nrhs
ferr( j ) = ferr( j ) / rowcnd
end do
end if
work( 1_${ik}$ ) = rpvgrw
! 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}$gesvx
#:endif
#:endfor
module subroutine stdlib${ii}$_cgesvx( fact, trans, n, nrhs, a, lda, af, ldaf, ipiv,equed, r, c, b, ldb, &
!! CGESVX uses the LU factorization to compute the solution to a complex
!! system of linear equations
!! A * X = B,
!! where A is an N-by-N matrix and X and B are N-by-NRHS matrices.
!! Error bounds on the solution and a condition estimate are also
!! provided.
x, 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, trans
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(inout) :: ipiv(*)
real(sp), intent(out) :: berr(*), ferr(*), rwork(*)
real(sp), intent(inout) :: c(*), r(*)
complex(sp), intent(inout) :: a(lda,*), af(ldaf,*), b(ldb,*)
complex(sp), intent(out) :: work(*), x(ldx,*)
! =====================================================================
! Local Scalars
logical(lk) :: colequ, equil, nofact, notran, rowequ
character :: norm
integer(${ik}$) :: i, infequ, j
real(sp) :: amax, anorm, bignum, colcnd, rcmax, rcmin, rowcnd, rpvgrw, smlnum
! Intrinsic Functions
! Executable Statements
info = 0_${ik}$
nofact = stdlib_lsame( fact, 'N' )
equil = stdlib_lsame( fact, 'E' )
notran = stdlib_lsame( trans, 'N' )
if( nofact .or. equil ) then
equed = 'N'
rowequ = .false.
colequ = .false.
else
rowequ = stdlib_lsame( equed, 'R' ) .or. stdlib_lsame( equed, 'B' )
colequ = stdlib_lsame( equed, 'C' ) .or. stdlib_lsame( equed, 'B' )
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.notran .and. .not.stdlib_lsame( trans, 'T' ) .and. .not.stdlib_lsame( &
trans, 'C' ) ) 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.( rowequ .or. colequ .or. stdlib_lsame( &
equed, 'N' ) ) ) then
info = -10_${ik}$
else
if( rowequ ) then
rcmin = bignum
rcmax = zero
do j = 1, n
rcmin = min( rcmin, r( j ) )
rcmax = max( rcmax, r( j ) )
end do
if( rcmin<=zero ) then
info = -11_${ik}$
else if( n>0_${ik}$ ) then
rowcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
else
rowcnd = one
end if
end if
if( colequ .and. info==0_${ik}$ ) then
rcmin = bignum
rcmax = zero
do j = 1, n
rcmin = min( rcmin, c( j ) )
rcmax = max( rcmax, c( j ) )
end do
if( rcmin<=zero ) then
info = -12_${ik}$
else if( n>0_${ik}$ ) then
colcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
else
colcnd = one
end if
end if
if( info==0_${ik}$ ) then
if( ldb<max( 1_${ik}$, n ) ) then
info = -14_${ik}$
else if( ldx<max( 1_${ik}$, n ) ) then
info = -16_${ik}$
end if
end if
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'CGESVX', -info )
return
end if
if( equil ) then
! compute row and column scalings to equilibrate the matrix a.
call stdlib${ii}$_cgeequ( n, n, a, lda, r, c, rowcnd, colcnd, amax, infequ )
if( infequ==0_${ik}$ ) then
! equilibrate the matrix.
call stdlib${ii}$_claqge( n, n, a, lda, r, c, rowcnd, colcnd, amax,equed )
rowequ = stdlib_lsame( equed, 'R' ) .or. stdlib_lsame( equed, 'B' )
colequ = stdlib_lsame( equed, 'C' ) .or. stdlib_lsame( equed, 'B' )
end if
end if
! scale the right hand side.
if( notran ) then
if( rowequ ) then
do j = 1, nrhs
do i = 1, n
b( i, j ) = r( i )*b( i, j )
end do
end do
end if
else if( colequ ) then
do j = 1, nrhs
do i = 1, n
b( i, j ) = c( i )*b( i, j )
end do
end do
end if
if( nofact .or. equil ) then
! compute the lu factorization of a.
call stdlib${ii}$_clacpy( 'FULL', n, n, a, lda, af, ldaf )
call stdlib${ii}$_cgetrf( n, n, af, ldaf, ipiv, info )
! return if info is non-zero.
if( info>0_${ik}$ ) then
! compute the reciprocal pivot growth factor of the
! leading rank-deficient info columns of a.
rpvgrw = stdlib${ii}$_clantr( 'M', 'U', 'N', info, info, af, ldaf,rwork )
if( rpvgrw==zero ) then
rpvgrw = one
else
rpvgrw = stdlib${ii}$_clange( 'M', n, info, a, lda, rwork ) /rpvgrw
end if
rwork( 1_${ik}$ ) = rpvgrw
rcond = zero
return
end if
end if
! compute the norm of the matrix a and the
! reciprocal pivot growth factor rpvgrw.
if( notran ) then
norm = '1'
else
norm = 'I'
end if
anorm = stdlib${ii}$_clange( norm, n, n, a, lda, rwork )
rpvgrw = stdlib${ii}$_clantr( 'M', 'U', 'N', n, n, af, ldaf, rwork )
if( rpvgrw==zero ) then
rpvgrw = one
else
rpvgrw = stdlib${ii}$_clange( 'M', n, n, a, lda, rwork ) / rpvgrw
end if
! compute the reciprocal of the condition number of a.
call stdlib${ii}$_cgecon( norm, 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}$_cgetrs( trans, n, nrhs, af, ldaf, ipiv, x, ldx, info )
! use iterative refinement to improve the computed solution and
! compute error bounds and backward error estimates for it.
call stdlib${ii}$_cgerfs( trans, n, nrhs, a, lda, af, ldaf, ipiv, b, ldb, x,ldx, ferr, berr, &
work, rwork, info )
! transform the solution matrix x to a solution of the original
! system.
if( notran ) then
if( colequ ) then
do j = 1, nrhs
do i = 1, n
x( i, j ) = c( i )*x( i, j )
end do
end do
do j = 1, nrhs
ferr( j ) = ferr( j ) / colcnd
end do
end if
else if( rowequ ) then
do j = 1, nrhs
do i = 1, n
x( i, j ) = r( i )*x( i, j )
end do
end do
do j = 1, nrhs
ferr( j ) = ferr( j ) / rowcnd
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}$
rwork( 1_${ik}$ ) = rpvgrw
return
end subroutine stdlib${ii}$_cgesvx
module subroutine stdlib${ii}$_zgesvx( fact, trans, n, nrhs, a, lda, af, ldaf, ipiv,equed, r, c, b, ldb, &
!! ZGESVX uses the LU factorization to compute the solution to a complex
!! system of linear equations
!! A * X = B,
!! where A is an N-by-N matrix and X and B are N-by-NRHS matrices.
!! Error bounds on the solution and a condition estimate are also
!! provided.
x, 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, trans
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(inout) :: ipiv(*)
real(dp), intent(out) :: berr(*), ferr(*), rwork(*)
real(dp), intent(inout) :: c(*), r(*)
complex(dp), intent(inout) :: a(lda,*), af(ldaf,*), b(ldb,*)
complex(dp), intent(out) :: work(*), x(ldx,*)
! =====================================================================
! Local Scalars
logical(lk) :: colequ, equil, nofact, notran, rowequ
character :: norm
integer(${ik}$) :: i, infequ, j
real(dp) :: amax, anorm, bignum, colcnd, rcmax, rcmin, rowcnd, rpvgrw, smlnum
! Intrinsic Functions
! Executable Statements
info = 0_${ik}$
nofact = stdlib_lsame( fact, 'N' )
equil = stdlib_lsame( fact, 'E' )
notran = stdlib_lsame( trans, 'N' )
if( nofact .or. equil ) then
equed = 'N'
rowequ = .false.
colequ = .false.
else
rowequ = stdlib_lsame( equed, 'R' ) .or. stdlib_lsame( equed, 'B' )
colequ = stdlib_lsame( equed, 'C' ) .or. stdlib_lsame( equed, 'B' )
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.notran .and. .not.stdlib_lsame( trans, 'T' ) .and. .not.stdlib_lsame( &
trans, 'C' ) ) 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.( rowequ .or. colequ .or. stdlib_lsame( &
equed, 'N' ) ) ) then
info = -10_${ik}$
else
if( rowequ ) then
rcmin = bignum
rcmax = zero
do j = 1, n
rcmin = min( rcmin, r( j ) )
rcmax = max( rcmax, r( j ) )
end do
if( rcmin<=zero ) then
info = -11_${ik}$
else if( n>0_${ik}$ ) then
rowcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
else
rowcnd = one
end if
end if
if( colequ .and. info==0_${ik}$ ) then
rcmin = bignum
rcmax = zero
do j = 1, n
rcmin = min( rcmin, c( j ) )
rcmax = max( rcmax, c( j ) )
end do
if( rcmin<=zero ) then
info = -12_${ik}$
else if( n>0_${ik}$ ) then
colcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
else
colcnd = one
end if
end if
if( info==0_${ik}$ ) then
if( ldb<max( 1_${ik}$, n ) ) then
info = -14_${ik}$
else if( ldx<max( 1_${ik}$, n ) ) then
info = -16_${ik}$
end if
end if
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'ZGESVX', -info )
return
end if
if( equil ) then
! compute row and column scalings to equilibrate the matrix a.
call stdlib${ii}$_zgeequ( n, n, a, lda, r, c, rowcnd, colcnd, amax, infequ )
if( infequ==0_${ik}$ ) then
! equilibrate the matrix.
call stdlib${ii}$_zlaqge( n, n, a, lda, r, c, rowcnd, colcnd, amax,equed )
rowequ = stdlib_lsame( equed, 'R' ) .or. stdlib_lsame( equed, 'B' )
colequ = stdlib_lsame( equed, 'C' ) .or. stdlib_lsame( equed, 'B' )
end if
end if
! scale the right hand side.
if( notran ) then
if( rowequ ) then
do j = 1, nrhs
do i = 1, n
b( i, j ) = r( i )*b( i, j )
end do
end do
end if
else if( colequ ) then
do j = 1, nrhs
do i = 1, n
b( i, j ) = c( i )*b( i, j )
end do
end do
end if
if( nofact .or. equil ) then
! compute the lu factorization of a.
call stdlib${ii}$_zlacpy( 'FULL', n, n, a, lda, af, ldaf )
call stdlib${ii}$_zgetrf( n, n, af, ldaf, ipiv, info )
! return if info is non-zero.
if( info>0_${ik}$ ) then
! compute the reciprocal pivot growth factor of the
! leading rank-deficient info columns of a.
rpvgrw = stdlib${ii}$_zlantr( 'M', 'U', 'N', info, info, af, ldaf,rwork )
if( rpvgrw==zero ) then
rpvgrw = one
else
rpvgrw = stdlib${ii}$_zlange( 'M', n, info, a, lda, rwork ) /rpvgrw
end if
rwork( 1_${ik}$ ) = rpvgrw
rcond = zero
return
end if
end if
! compute the norm of the matrix a and the
! reciprocal pivot growth factor rpvgrw.
if( notran ) then
norm = '1'
else
norm = 'I'
end if
anorm = stdlib${ii}$_zlange( norm, n, n, a, lda, rwork )
rpvgrw = stdlib${ii}$_zlantr( 'M', 'U', 'N', n, n, af, ldaf, rwork )
if( rpvgrw==zero ) then
rpvgrw = one
else
rpvgrw = stdlib${ii}$_zlange( 'M', n, n, a, lda, rwork ) / rpvgrw
end if
! compute the reciprocal of the condition number of a.
call stdlib${ii}$_zgecon( norm, 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}$_zgetrs( trans, n, nrhs, af, ldaf, ipiv, x, ldx, info )
! use iterative refinement to improve the computed solution and
! compute error bounds and backward error estimates for it.
call stdlib${ii}$_zgerfs( trans, n, nrhs, a, lda, af, ldaf, ipiv, b, ldb, x,ldx, ferr, berr, &
work, rwork, info )
! transform the solution matrix x to a solution of the original
! system.
if( notran ) then
if( colequ ) then
do j = 1, nrhs
do i = 1, n
x( i, j ) = c( i )*x( i, j )
end do
end do
do j = 1, nrhs
ferr( j ) = ferr( j ) / colcnd
end do
end if
else if( rowequ ) then
do j = 1, nrhs
do i = 1, n
x( i, j ) = r( i )*x( i, j )
end do
end do
do j = 1, nrhs
ferr( j ) = ferr( j ) / rowcnd
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}$
rwork( 1_${ik}$ ) = rpvgrw
return
end subroutine stdlib${ii}$_zgesvx
#:for ck,ct,ci in CMPLX_KINDS_TYPES
#:if not ck in ["sp","dp"]
module subroutine stdlib${ii}$_${ci}$gesvx( fact, trans, n, nrhs, a, lda, af, ldaf, ipiv,equed, r, c, b, ldb, &
!! ZGESVX: uses the LU factorization to compute the solution to a complex
!! system of linear equations
!! A * X = B,
!! where A is an N-by-N matrix and X and B are N-by-NRHS matrices.
!! Error bounds on the solution and a condition estimate are also
!! provided.
x, 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, trans
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: lda, ldaf, ldb, ldx, n, nrhs
real(${ck}$), intent(out) :: rcond
! Array Arguments
integer(${ik}$), intent(inout) :: ipiv(*)
real(${ck}$), intent(out) :: berr(*), ferr(*), rwork(*)
real(${ck}$), intent(inout) :: c(*), r(*)
complex(${ck}$), intent(inout) :: a(lda,*), af(ldaf,*), b(ldb,*)
complex(${ck}$), intent(out) :: work(*), x(ldx,*)
! =====================================================================
! Local Scalars
logical(lk) :: colequ, equil, nofact, notran, rowequ
character :: norm
integer(${ik}$) :: i, infequ, j
real(${ck}$) :: amax, anorm, bignum, colcnd, rcmax, rcmin, rowcnd, rpvgrw, smlnum
! Intrinsic Functions
! Executable Statements
info = 0_${ik}$
nofact = stdlib_lsame( fact, 'N' )
equil = stdlib_lsame( fact, 'E' )
notran = stdlib_lsame( trans, 'N' )
if( nofact .or. equil ) then
equed = 'N'
rowequ = .false.
colequ = .false.
else
rowequ = stdlib_lsame( equed, 'R' ) .or. stdlib_lsame( equed, 'B' )
colequ = stdlib_lsame( equed, 'C' ) .or. stdlib_lsame( equed, 'B' )
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.notran .and. .not.stdlib_lsame( trans, 'T' ) .and. .not.stdlib_lsame( &
trans, 'C' ) ) 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.( rowequ .or. colequ .or. stdlib_lsame( &
equed, 'N' ) ) ) then
info = -10_${ik}$
else
if( rowequ ) then
rcmin = bignum
rcmax = zero
do j = 1, n
rcmin = min( rcmin, r( j ) )
rcmax = max( rcmax, r( j ) )
end do
if( rcmin<=zero ) then
info = -11_${ik}$
else if( n>0_${ik}$ ) then
rowcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
else
rowcnd = one
end if
end if
if( colequ .and. info==0_${ik}$ ) then
rcmin = bignum
rcmax = zero
do j = 1, n
rcmin = min( rcmin, c( j ) )
rcmax = max( rcmax, c( j ) )
end do
if( rcmin<=zero ) then
info = -12_${ik}$
else if( n>0_${ik}$ ) then
colcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
else
colcnd = one
end if
end if
if( info==0_${ik}$ ) then
if( ldb<max( 1_${ik}$, n ) ) then
info = -14_${ik}$
else if( ldx<max( 1_${ik}$, n ) ) then
info = -16_${ik}$
end if
end if
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'ZGESVX', -info )
return
end if
if( equil ) then
! compute row and column scalings to equilibrate the matrix a.
call stdlib${ii}$_${ci}$geequ( n, n, a, lda, r, c, rowcnd, colcnd, amax, infequ )
if( infequ==0_${ik}$ ) then
! equilibrate the matrix.
call stdlib${ii}$_${ci}$laqge( n, n, a, lda, r, c, rowcnd, colcnd, amax,equed )
rowequ = stdlib_lsame( equed, 'R' ) .or. stdlib_lsame( equed, 'B' )
colequ = stdlib_lsame( equed, 'C' ) .or. stdlib_lsame( equed, 'B' )
end if
end if
! scale the right hand side.
if( notran ) then
if( rowequ ) then
do j = 1, nrhs
do i = 1, n
b( i, j ) = r( i )*b( i, j )
end do
end do
end if
else if( colequ ) then
do j = 1, nrhs
do i = 1, n
b( i, j ) = c( i )*b( i, j )
end do
end do
end if
if( nofact .or. equil ) then
! compute the lu factorization of a.
call stdlib${ii}$_${ci}$lacpy( 'FULL', n, n, a, lda, af, ldaf )
call stdlib${ii}$_${ci}$getrf( n, n, af, ldaf, ipiv, info )
! return if info is non-zero.
if( info>0_${ik}$ ) then
! compute the reciprocal pivot growth factor of the
! leading rank-deficient info columns of a.
rpvgrw = stdlib${ii}$_${ci}$lantr( 'M', 'U', 'N', info, info, af, ldaf,rwork )
if( rpvgrw==zero ) then
rpvgrw = one
else
rpvgrw = stdlib${ii}$_${ci}$lange( 'M', n, info, a, lda, rwork ) /rpvgrw
end if
rwork( 1_${ik}$ ) = rpvgrw
rcond = zero
return
end if
end if
! compute the norm of the matrix a and the
! reciprocal pivot growth factor rpvgrw.
if( notran ) then
norm = '1'
else
norm = 'I'
end if
anorm = stdlib${ii}$_${ci}$lange( norm, n, n, a, lda, rwork )
rpvgrw = stdlib${ii}$_${ci}$lantr( 'M', 'U', 'N', n, n, af, ldaf, rwork )
if( rpvgrw==zero ) then
rpvgrw = one
else
rpvgrw = stdlib${ii}$_${ci}$lange( 'M', n, n, a, lda, rwork ) / rpvgrw
end if
! compute the reciprocal of the condition number of a.
call stdlib${ii}$_${ci}$gecon( norm, 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}$getrs( trans, n, nrhs, af, ldaf, ipiv, 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}$gerfs( trans, n, nrhs, a, lda, af, ldaf, ipiv, b, ldb, x,ldx, ferr, berr, &
work, rwork, info )
! transform the solution matrix x to a solution of the original
! system.
if( notran ) then
if( colequ ) then
do j = 1, nrhs
do i = 1, n
x( i, j ) = c( i )*x( i, j )
end do
end do
do j = 1, nrhs
ferr( j ) = ferr( j ) / colcnd
end do
end if
else if( rowequ ) then
do j = 1, nrhs
do i = 1, n
x( i, j ) = r( i )*x( i, j )
end do
end do
do j = 1, nrhs
ferr( j ) = ferr( j ) / rowcnd
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}$
rwork( 1_${ik}$ ) = rpvgrw
return
end subroutine stdlib${ii}$_${ci}$gesvx
#:endif
#:endfor
pure module subroutine stdlib${ii}$_sgbsv( n, kl, ku, nrhs, ab, ldab, ipiv, b, ldb, info )
!! SGBSV computes the solution to a real system of linear equations
!! A * X = B, where A is a band matrix of order N with KL subdiagonals
!! and KU superdiagonals, and X and B are N-by-NRHS matrices.
!! The LU decomposition with partial pivoting and row interchanges is
!! used to factor A as A = L * U, where L is a product of permutation
!! and unit lower triangular matrices with KL subdiagonals, and U is
!! upper triangular with KL+KU superdiagonals. 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
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: kl, ku, ldab, ldb, n, nrhs
! Array Arguments
integer(${ik}$), intent(out) :: ipiv(*)
real(sp), intent(inout) :: ab(ldab,*), b(ldb,*)
! =====================================================================
! Intrinsic Functions
! Executable Statements
! test the input parameters.
info = 0_${ik}$
if( n<0_${ik}$ ) then
info = -1_${ik}$
else if( kl<0_${ik}$ ) then
info = -2_${ik}$
else if( ku<0_${ik}$ ) then
info = -3_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -4_${ik}$
else if( ldab<2_${ik}$*kl+ku+1 ) then
info = -6_${ik}$
else if( ldb<max( n, 1_${ik}$ ) ) then
info = -9_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'SGBSV ', -info )
return
end if
! compute the lu factorization of the band matrix a.
call stdlib${ii}$_sgbtrf( n, n, kl, ku, ab, ldab, ipiv, info )
if( info==0_${ik}$ ) then
! solve the system a*x = b, overwriting b with x.
call stdlib${ii}$_sgbtrs( 'NO TRANSPOSE', n, kl, ku, nrhs, ab, ldab, ipiv,b, ldb, info )
end if
return
end subroutine stdlib${ii}$_sgbsv
pure module subroutine stdlib${ii}$_dgbsv( n, kl, ku, nrhs, ab, ldab, ipiv, b, ldb, info )
!! DGBSV computes the solution to a real system of linear equations
!! A * X = B, where A is a band matrix of order N with KL subdiagonals
!! and KU superdiagonals, and X and B are N-by-NRHS matrices.
!! The LU decomposition with partial pivoting and row interchanges is
!! used to factor A as A = L * U, where L is a product of permutation
!! and unit lower triangular matrices with KL subdiagonals, and U is
!! upper triangular with KL+KU superdiagonals. 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
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: kl, ku, ldab, ldb, n, nrhs
! Array Arguments
integer(${ik}$), intent(out) :: ipiv(*)
real(dp), intent(inout) :: ab(ldab,*), b(ldb,*)
! =====================================================================
! Intrinsic Functions
! Executable Statements
! test the input parameters.
info = 0_${ik}$
if( n<0_${ik}$ ) then
info = -1_${ik}$
else if( kl<0_${ik}$ ) then
info = -2_${ik}$
else if( ku<0_${ik}$ ) then
info = -3_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -4_${ik}$
else if( ldab<2_${ik}$*kl+ku+1 ) then
info = -6_${ik}$
else if( ldb<max( n, 1_${ik}$ ) ) then
info = -9_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'DGBSV ', -info )
return
end if
! compute the lu factorization of the band matrix a.
call stdlib${ii}$_dgbtrf( n, n, kl, ku, ab, ldab, ipiv, info )
if( info==0_${ik}$ ) then
! solve the system a*x = b, overwriting b with x.
call stdlib${ii}$_dgbtrs( 'NO TRANSPOSE', n, kl, ku, nrhs, ab, ldab, ipiv,b, ldb, info )
end if
return
end subroutine stdlib${ii}$_dgbsv
#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
pure module subroutine stdlib${ii}$_${ri}$gbsv( n, kl, ku, nrhs, ab, ldab, ipiv, b, ldb, info )
!! DGBSV: computes the solution to a real system of linear equations
!! A * X = B, where A is a band matrix of order N with KL subdiagonals
!! and KU superdiagonals, and X and B are N-by-NRHS matrices.
!! The LU decomposition with partial pivoting and row interchanges is
!! used to factor A as A = L * U, where L is a product of permutation
!! and unit lower triangular matrices with KL subdiagonals, and U is
!! upper triangular with KL+KU superdiagonals. 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
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: kl, ku, ldab, ldb, n, nrhs
! Array Arguments
integer(${ik}$), intent(out) :: ipiv(*)
real(${rk}$), intent(inout) :: ab(ldab,*), b(ldb,*)
! =====================================================================
! Intrinsic Functions
! Executable Statements
! test the input parameters.
info = 0_${ik}$
if( n<0_${ik}$ ) then
info = -1_${ik}$
else if( kl<0_${ik}$ ) then
info = -2_${ik}$
else if( ku<0_${ik}$ ) then
info = -3_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -4_${ik}$
else if( ldab<2_${ik}$*kl+ku+1 ) then
info = -6_${ik}$
else if( ldb<max( n, 1_${ik}$ ) ) then
info = -9_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'DGBSV ', -info )
return
end if
! compute the lu factorization of the band matrix a.
call stdlib${ii}$_${ri}$gbtrf( n, n, kl, ku, ab, ldab, ipiv, info )
if( info==0_${ik}$ ) then
! solve the system a*x = b, overwriting b with x.
call stdlib${ii}$_${ri}$gbtrs( 'NO TRANSPOSE', n, kl, ku, nrhs, ab, ldab, ipiv,b, ldb, info )
end if
return
end subroutine stdlib${ii}$_${ri}$gbsv
#:endif
#:endfor
pure module subroutine stdlib${ii}$_cgbsv( n, kl, ku, nrhs, ab, ldab, ipiv, b, ldb, info )
!! CGBSV computes the solution to a complex system of linear equations
!! A * X = B, where A is a band matrix of order N with KL subdiagonals
!! and KU superdiagonals, and X and B are N-by-NRHS matrices.
!! The LU decomposition with partial pivoting and row interchanges is
!! used to factor A as A = L * U, where L is a product of permutation
!! and unit lower triangular matrices with KL subdiagonals, and U is
!! upper triangular with KL+KU superdiagonals. 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
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: kl, ku, ldab, ldb, n, nrhs
! Array Arguments
integer(${ik}$), intent(out) :: ipiv(*)
complex(sp), intent(inout) :: ab(ldab,*), b(ldb,*)
! =====================================================================
! Intrinsic Functions
! Executable Statements
! test the input parameters.
info = 0_${ik}$
if( n<0_${ik}$ ) then
info = -1_${ik}$
else if( kl<0_${ik}$ ) then
info = -2_${ik}$
else if( ku<0_${ik}$ ) then
info = -3_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -4_${ik}$
else if( ldab<2_${ik}$*kl+ku+1 ) then
info = -6_${ik}$
else if( ldb<max( n, 1_${ik}$ ) ) then
info = -9_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'CGBSV ', -info )
return
end if
! compute the lu factorization of the band matrix a.
call stdlib${ii}$_cgbtrf( n, n, kl, ku, ab, ldab, ipiv, info )
if( info==0_${ik}$ ) then
! solve the system a*x = b, overwriting b with x.
call stdlib${ii}$_cgbtrs( 'NO TRANSPOSE', n, kl, ku, nrhs, ab, ldab, ipiv,b, ldb, info )
end if
return
end subroutine stdlib${ii}$_cgbsv
pure module subroutine stdlib${ii}$_zgbsv( n, kl, ku, nrhs, ab, ldab, ipiv, b, ldb, info )
!! ZGBSV computes the solution to a complex system of linear equations
!! A * X = B, where A is a band matrix of order N with KL subdiagonals
!! and KU superdiagonals, and X and B are N-by-NRHS matrices.
!! The LU decomposition with partial pivoting and row interchanges is
!! used to factor A as A = L * U, where L is a product of permutation
!! and unit lower triangular matrices with KL subdiagonals, and U is
!! upper triangular with KL+KU superdiagonals. 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
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: kl, ku, ldab, ldb, n, nrhs
! Array Arguments
integer(${ik}$), intent(out) :: ipiv(*)
complex(dp), intent(inout) :: ab(ldab,*), b(ldb,*)
! =====================================================================
! Intrinsic Functions
! Executable Statements
! test the input parameters.
info = 0_${ik}$
if( n<0_${ik}$ ) then
info = -1_${ik}$
else if( kl<0_${ik}$ ) then
info = -2_${ik}$
else if( ku<0_${ik}$ ) then
info = -3_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -4_${ik}$
else if( ldab<2_${ik}$*kl+ku+1 ) then
info = -6_${ik}$
else if( ldb<max( n, 1_${ik}$ ) ) then
info = -9_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'ZGBSV ', -info )
return
end if
! compute the lu factorization of the band matrix a.
call stdlib${ii}$_zgbtrf( n, n, kl, ku, ab, ldab, ipiv, info )
if( info==0_${ik}$ ) then
! solve the system a*x = b, overwriting b with x.
call stdlib${ii}$_zgbtrs( 'NO TRANSPOSE', n, kl, ku, nrhs, ab, ldab, ipiv,b, ldb, info )
end if
return
end subroutine stdlib${ii}$_zgbsv
#:for ck,ct,ci in CMPLX_KINDS_TYPES
#:if not ck in ["sp","dp"]
pure module subroutine stdlib${ii}$_${ci}$gbsv( n, kl, ku, nrhs, ab, ldab, ipiv, b, ldb, info )
!! ZGBSV: computes the solution to a complex system of linear equations
!! A * X = B, where A is a band matrix of order N with KL subdiagonals
!! and KU superdiagonals, and X and B are N-by-NRHS matrices.
!! The LU decomposition with partial pivoting and row interchanges is
!! used to factor A as A = L * U, where L is a product of permutation
!! and unit lower triangular matrices with KL subdiagonals, and U is
!! upper triangular with KL+KU superdiagonals. 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
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: kl, ku, ldab, ldb, n, nrhs
! Array Arguments
integer(${ik}$), intent(out) :: ipiv(*)
complex(${ck}$), intent(inout) :: ab(ldab,*), b(ldb,*)
! =====================================================================
! Intrinsic Functions
! Executable Statements
! test the input parameters.
info = 0_${ik}$
if( n<0_${ik}$ ) then
info = -1_${ik}$
else if( kl<0_${ik}$ ) then
info = -2_${ik}$
else if( ku<0_${ik}$ ) then
info = -3_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -4_${ik}$
else if( ldab<2_${ik}$*kl+ku+1 ) then
info = -6_${ik}$
else if( ldb<max( n, 1_${ik}$ ) ) then
info = -9_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'ZGBSV ', -info )
return
end if
! compute the lu factorization of the band matrix a.
call stdlib${ii}$_${ci}$gbtrf( n, n, kl, ku, ab, ldab, ipiv, info )
if( info==0_${ik}$ ) then
! solve the system a*x = b, overwriting b with x.
call stdlib${ii}$_${ci}$gbtrs( 'NO TRANSPOSE', n, kl, ku, nrhs, ab, ldab, ipiv,b, ldb, info )
end if
return
end subroutine stdlib${ii}$_${ci}$gbsv
#:endif
#:endfor
module subroutine stdlib${ii}$_sgbsvx( fact, trans, n, kl, ku, nrhs, ab, ldab, afb,ldafb, ipiv, equed, r, &
!! SGBSVX uses the LU factorization to compute the solution to a real
!! system of linear equations A * X = B, A**T * X = B, or A**H * X = B,
!! where A is a band matrix of order N with KL subdiagonals and KU
!! superdiagonals, and X and B are N-by-NRHS matrices.
!! Error bounds on the solution and a condition estimate are also
!! provided.
c, b, ldb, x, 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, trans
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: kl, ku, ldab, ldafb, ldb, ldx, n, nrhs
real(sp), intent(out) :: rcond
! Array Arguments
integer(${ik}$), intent(inout) :: ipiv(*)
integer(${ik}$), intent(out) :: iwork(*)
real(sp), intent(inout) :: ab(ldab,*), afb(ldafb,*), b(ldb,*), c(*), r(*)
real(sp), intent(out) :: berr(*), ferr(*), work(*), x(ldx,*)
! =====================================================================
! moved setting of info = n+1 so info does not subsequently get
! overwritten. sven, 17 mar 05.
! =====================================================================
! Local Scalars
logical(lk) :: colequ, equil, nofact, notran, rowequ
character :: norm
integer(${ik}$) :: i, infequ, j, j1, j2
real(sp) :: amax, anorm, bignum, colcnd, rcmax, rcmin, rowcnd, rpvgrw, smlnum
! Intrinsic Functions
! Executable Statements
info = 0_${ik}$
nofact = stdlib_lsame( fact, 'N' )
equil = stdlib_lsame( fact, 'E' )
notran = stdlib_lsame( trans, 'N' )
if( nofact .or. equil ) then
equed = 'N'
rowequ = .false.
colequ = .false.
else
rowequ = stdlib_lsame( equed, 'R' ) .or. stdlib_lsame( equed, 'B' )
colequ = stdlib_lsame( equed, 'C' ) .or. stdlib_lsame( equed, 'B' )
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.notran .and. .not.stdlib_lsame( trans, 'T' ) .and. .not.stdlib_lsame( &
trans, 'C' ) ) then
info = -2_${ik}$
else if( n<0_${ik}$ ) then
info = -3_${ik}$
else if( kl<0_${ik}$ ) then
info = -4_${ik}$
else if( ku<0_${ik}$ ) then
info = -5_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -6_${ik}$
else if( ldab<kl+ku+1 ) then
info = -8_${ik}$
else if( ldafb<2_${ik}$*kl+ku+1 ) then
info = -10_${ik}$
else if( stdlib_lsame( fact, 'F' ) .and. .not.( rowequ .or. colequ .or. stdlib_lsame( &
equed, 'N' ) ) ) then
info = -12_${ik}$
else
if( rowequ ) then
rcmin = bignum
rcmax = zero
do j = 1, n
rcmin = min( rcmin, r( j ) )
rcmax = max( rcmax, r( j ) )
end do
if( rcmin<=zero ) then
info = -13_${ik}$
else if( n>0_${ik}$ ) then
rowcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
else
rowcnd = one
end if
end if
if( colequ .and. info==0_${ik}$ ) then
rcmin = bignum
rcmax = zero
do j = 1, n
rcmin = min( rcmin, c( j ) )
rcmax = max( rcmax, c( j ) )
end do
if( rcmin<=zero ) then
info = -14_${ik}$
else if( n>0_${ik}$ ) then
colcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
else
colcnd = one
end if
end if
if( info==0_${ik}$ ) then
if( ldb<max( 1_${ik}$, n ) ) then
info = -16_${ik}$
else if( ldx<max( 1_${ik}$, n ) ) then
info = -18_${ik}$
end if
end if
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'SGBSVX', -info )
return
end if
if( equil ) then
! compute row and column scalings to equilibrate the matrix a.
call stdlib${ii}$_sgbequ( n, n, kl, ku, ab, ldab, r, c, rowcnd, colcnd,amax, infequ )
if( infequ==0_${ik}$ ) then
! equilibrate the matrix.
call stdlib${ii}$_slaqgb( n, n, kl, ku, ab, ldab, r, c, rowcnd, colcnd,amax, equed )
rowequ = stdlib_lsame( equed, 'R' ) .or. stdlib_lsame( equed, 'B' )
colequ = stdlib_lsame( equed, 'C' ) .or. stdlib_lsame( equed, 'B' )
end if
end if
! scale the right hand side.
if( notran ) then
if( rowequ ) then
do j = 1, nrhs
do i = 1, n
b( i, j ) = r( i )*b( i, j )
end do
end do
end if
else if( colequ ) then
do j = 1, nrhs
do i = 1, n
b( i, j ) = c( i )*b( i, j )
end do
end do
end if
if( nofact .or. equil ) then
! compute the lu factorization of the band matrix a.
do j = 1, n
j1 = max( j-ku, 1_${ik}$ )
j2 = min( j+kl, n )
call stdlib${ii}$_scopy( j2-j1+1, ab( ku+1-j+j1, j ), 1_${ik}$,afb( kl+ku+1-j+j1, j ), 1_${ik}$ )
end do
call stdlib${ii}$_sgbtrf( n, n, kl, ku, afb, ldafb, ipiv, info )
! return if info is non-zero.
if( info>0_${ik}$ ) then
! compute the reciprocal pivot growth factor of the
! leading rank-deficient info columns of a.
anorm = zero
do j = 1, info
do i = max( ku+2-j, 1 ), min( n+ku+1-j, kl+ku+1 )
anorm = max( anorm, abs( ab( i, j ) ) )
end do
end do
rpvgrw = stdlib${ii}$_slantb( 'M', 'U', 'N', info, min( info-1, kl+ku ),afb( max( 1_${ik}$, &
kl+ku+2-info ), 1_${ik}$ ), ldafb,work )
if( rpvgrw==zero ) then
rpvgrw = one
else
rpvgrw = anorm / rpvgrw
end if
work( 1_${ik}$ ) = rpvgrw
rcond = zero
return
end if
end if
! compute the norm of the matrix a and the
! reciprocal pivot growth factor rpvgrw.
if( notran ) then
norm = '1'
else
norm = 'I'
end if
anorm = stdlib${ii}$_slangb( norm, n, kl, ku, ab, ldab, work )
rpvgrw = stdlib${ii}$_slantb( 'M', 'U', 'N', n, kl+ku, afb, ldafb, work )
if( rpvgrw==zero ) then
rpvgrw = one
else
rpvgrw = stdlib${ii}$_slangb( 'M', n, kl, ku, ab, ldab, work ) / rpvgrw
end if
! compute the reciprocal of the condition number of a.
call stdlib${ii}$_sgbcon( norm, n, kl, ku, afb, ldafb, ipiv, anorm, rcond,work, iwork, info )
! compute the solution matrix x.
call stdlib${ii}$_slacpy( 'FULL', n, nrhs, b, ldb, x, ldx )
call stdlib${ii}$_sgbtrs( trans, n, kl, ku, nrhs, afb, ldafb, ipiv, x, ldx,info )
! use iterative refinement to improve the computed solution and
! compute error bounds and backward error estimates for it.
call stdlib${ii}$_sgbrfs( trans, n, kl, ku, nrhs, ab, ldab, afb, ldafb, ipiv,b, ldb, x, ldx, &
ferr, berr, work, iwork, info )
! transform the solution matrix x to a solution of the original
! system.
if( notran ) then
if( colequ ) then
do j = 1, nrhs
do i = 1, n
x( i, j ) = c( i )*x( i, j )
end do
end do
do j = 1, nrhs
ferr( j ) = ferr( j ) / colcnd
end do
end if
else if( rowequ ) then
do j = 1, nrhs
do i = 1, n
x( i, j ) = r( i )*x( i, j )
end do
end do
do j = 1, nrhs
ferr( j ) = ferr( j ) / rowcnd
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}$
work( 1_${ik}$ ) = rpvgrw
return
end subroutine stdlib${ii}$_sgbsvx
module subroutine stdlib${ii}$_dgbsvx( fact, trans, n, kl, ku, nrhs, ab, ldab, afb,ldafb, ipiv, equed, r, &
!! DGBSVX uses the LU factorization to compute the solution to a real
!! system of linear equations A * X = B, A**T * X = B, or A**H * X = B,
!! where A is a band matrix of order N with KL subdiagonals and KU
!! superdiagonals, and X and B are N-by-NRHS matrices.
!! Error bounds on the solution and a condition estimate are also
!! provided.
c, b, ldb, x, 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, trans
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: kl, ku, ldab, ldafb, ldb, ldx, n, nrhs
real(dp), intent(out) :: rcond
! Array Arguments
integer(${ik}$), intent(inout) :: ipiv(*)
integer(${ik}$), intent(out) :: iwork(*)
real(dp), intent(inout) :: ab(ldab,*), afb(ldafb,*), b(ldb,*), c(*), r(*)
real(dp), intent(out) :: berr(*), ferr(*), work(*), x(ldx,*)
! =====================================================================
! Local Scalars
logical(lk) :: colequ, equil, nofact, notran, rowequ
character :: norm
integer(${ik}$) :: i, infequ, j, j1, j2
real(dp) :: amax, anorm, bignum, colcnd, rcmax, rcmin, rowcnd, rpvgrw, smlnum
! Intrinsic Functions
! Executable Statements
info = 0_${ik}$
nofact = stdlib_lsame( fact, 'N' )
equil = stdlib_lsame( fact, 'E' )
notran = stdlib_lsame( trans, 'N' )
if( nofact .or. equil ) then
equed = 'N'
rowequ = .false.
colequ = .false.
else
rowequ = stdlib_lsame( equed, 'R' ) .or. stdlib_lsame( equed, 'B' )
colequ = stdlib_lsame( equed, 'C' ) .or. stdlib_lsame( equed, 'B' )
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.notran .and. .not.stdlib_lsame( trans, 'T' ) .and. .not.stdlib_lsame( &
trans, 'C' ) ) then
info = -2_${ik}$
else if( n<0_${ik}$ ) then
info = -3_${ik}$
else if( kl<0_${ik}$ ) then
info = -4_${ik}$
else if( ku<0_${ik}$ ) then
info = -5_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -6_${ik}$
else if( ldab<kl+ku+1 ) then
info = -8_${ik}$
else if( ldafb<2_${ik}$*kl+ku+1 ) then
info = -10_${ik}$
else if( stdlib_lsame( fact, 'F' ) .and. .not.( rowequ .or. colequ .or. stdlib_lsame( &
equed, 'N' ) ) ) then
info = -12_${ik}$
else
if( rowequ ) then
rcmin = bignum
rcmax = zero
do j = 1, n
rcmin = min( rcmin, r( j ) )
rcmax = max( rcmax, r( j ) )
end do
if( rcmin<=zero ) then
info = -13_${ik}$
else if( n>0_${ik}$ ) then
rowcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
else
rowcnd = one
end if
end if
if( colequ .and. info==0_${ik}$ ) then
rcmin = bignum
rcmax = zero
do j = 1, n
rcmin = min( rcmin, c( j ) )
rcmax = max( rcmax, c( j ) )
end do
if( rcmin<=zero ) then
info = -14_${ik}$
else if( n>0_${ik}$ ) then
colcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
else
colcnd = one
end if
end if
if( info==0_${ik}$ ) then
if( ldb<max( 1_${ik}$, n ) ) then
info = -16_${ik}$
else if( ldx<max( 1_${ik}$, n ) ) then
info = -18_${ik}$
end if
end if
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'DGBSVX', -info )
return
end if
if( equil ) then
! compute row and column scalings to equilibrate the matrix a.
call stdlib${ii}$_dgbequ( n, n, kl, ku, ab, ldab, r, c, rowcnd, colcnd,amax, infequ )
if( infequ==0_${ik}$ ) then
! equilibrate the matrix.
call stdlib${ii}$_dlaqgb( n, n, kl, ku, ab, ldab, r, c, rowcnd, colcnd,amax, equed )
rowequ = stdlib_lsame( equed, 'R' ) .or. stdlib_lsame( equed, 'B' )
colequ = stdlib_lsame( equed, 'C' ) .or. stdlib_lsame( equed, 'B' )
end if
end if
! scale the right hand side.
if( notran ) then
if( rowequ ) then
do j = 1, nrhs
do i = 1, n
b( i, j ) = r( i )*b( i, j )
end do
end do
end if
else if( colequ ) then
do j = 1, nrhs
do i = 1, n
b( i, j ) = c( i )*b( i, j )
end do
end do
end if
if( nofact .or. equil ) then
! compute the lu factorization of the band matrix a.
do j = 1, n
j1 = max( j-ku, 1_${ik}$ )
j2 = min( j+kl, n )
call stdlib${ii}$_dcopy( j2-j1+1, ab( ku+1-j+j1, j ), 1_${ik}$,afb( kl+ku+1-j+j1, j ), 1_${ik}$ )
end do
call stdlib${ii}$_dgbtrf( n, n, kl, ku, afb, ldafb, ipiv, info )
! return if info is non-zero.
if( info>0_${ik}$ ) then
! compute the reciprocal pivot growth factor of the
! leading rank-deficient info columns of a.
anorm = zero
do j = 1, info
do i = max( ku+2-j, 1 ), min( n+ku+1-j, kl+ku+1 )
anorm = max( anorm, abs( ab( i, j ) ) )
end do
end do
rpvgrw = stdlib${ii}$_dlantb( 'M', 'U', 'N', info, min( info-1, kl+ku ),afb( max( 1_${ik}$, &
kl+ku+2-info ), 1_${ik}$ ), ldafb,work )
if( rpvgrw==zero ) then
rpvgrw = one
else
rpvgrw = anorm / rpvgrw
end if
work( 1_${ik}$ ) = rpvgrw
rcond = zero
return
end if
end if
! compute the norm of the matrix a and the
! reciprocal pivot growth factor rpvgrw.
if( notran ) then
norm = '1'
else
norm = 'I'
end if
anorm = stdlib${ii}$_dlangb( norm, n, kl, ku, ab, ldab, work )
rpvgrw = stdlib${ii}$_dlantb( 'M', 'U', 'N', n, kl+ku, afb, ldafb, work )
if( rpvgrw==zero ) then
rpvgrw = one
else
rpvgrw = stdlib${ii}$_dlangb( 'M', n, kl, ku, ab, ldab, work ) / rpvgrw
end if
! compute the reciprocal of the condition number of a.
call stdlib${ii}$_dgbcon( norm, n, kl, ku, afb, ldafb, ipiv, anorm, rcond,work, iwork, info )
! compute the solution matrix x.
call stdlib${ii}$_dlacpy( 'FULL', n, nrhs, b, ldb, x, ldx )
call stdlib${ii}$_dgbtrs( trans, n, kl, ku, nrhs, afb, ldafb, ipiv, x, ldx,info )
! use iterative refinement to improve the computed solution and
! compute error bounds and backward error estimates for it.
call stdlib${ii}$_dgbrfs( trans, n, kl, ku, nrhs, ab, ldab, afb, ldafb, ipiv,b, ldb, x, ldx, &
ferr, berr, work, iwork, info )
! transform the solution matrix x to a solution of the original
! system.
if( notran ) then
if( colequ ) then
do j = 1, nrhs
do i = 1, n
x( i, j ) = c( i )*x( i, j )
end do
end do
do j = 1, nrhs
ferr( j ) = ferr( j ) / colcnd
end do
end if
else if( rowequ ) then
do j = 1, nrhs
do i = 1, n
x( i, j ) = r( i )*x( i, j )
end do
end do
do j = 1, nrhs
ferr( j ) = ferr( j ) / rowcnd
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}$
work( 1_${ik}$ ) = rpvgrw
return
end subroutine stdlib${ii}$_dgbsvx
#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
module subroutine stdlib${ii}$_${ri}$gbsvx( fact, trans, n, kl, ku, nrhs, ab, ldab, afb,ldafb, ipiv, equed, r, &
!! DGBSVX: uses the LU factorization to compute the solution to a real
!! system of linear equations A * X = B, A**T * X = B, or A**H * X = B,
!! where A is a band matrix of order N with KL subdiagonals and KU
!! superdiagonals, and X and B are N-by-NRHS matrices.
!! Error bounds on the solution and a condition estimate are also
!! provided.
c, b, ldb, x, 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, trans
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: kl, ku, ldab, ldafb, ldb, ldx, n, nrhs
real(${rk}$), intent(out) :: rcond
! Array Arguments
integer(${ik}$), intent(inout) :: ipiv(*)
integer(${ik}$), intent(out) :: iwork(*)
real(${rk}$), intent(inout) :: ab(ldab,*), afb(ldafb,*), b(ldb,*), c(*), r(*)
real(${rk}$), intent(out) :: berr(*), ferr(*), work(*), x(ldx,*)
! =====================================================================
! Local Scalars
logical(lk) :: colequ, equil, nofact, notran, rowequ
character :: norm
integer(${ik}$) :: i, infequ, j, j1, j2
real(${rk}$) :: amax, anorm, bignum, colcnd, rcmax, rcmin, rowcnd, rpvgrw, smlnum
! Intrinsic Functions
! Executable Statements
info = 0_${ik}$
nofact = stdlib_lsame( fact, 'N' )
equil = stdlib_lsame( fact, 'E' )
notran = stdlib_lsame( trans, 'N' )
if( nofact .or. equil ) then
equed = 'N'
rowequ = .false.
colequ = .false.
else
rowequ = stdlib_lsame( equed, 'R' ) .or. stdlib_lsame( equed, 'B' )
colequ = stdlib_lsame( equed, 'C' ) .or. stdlib_lsame( equed, 'B' )
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.notran .and. .not.stdlib_lsame( trans, 'T' ) .and. .not.stdlib_lsame( &
trans, 'C' ) ) then
info = -2_${ik}$
else if( n<0_${ik}$ ) then
info = -3_${ik}$
else if( kl<0_${ik}$ ) then
info = -4_${ik}$
else if( ku<0_${ik}$ ) then
info = -5_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -6_${ik}$
else if( ldab<kl+ku+1 ) then
info = -8_${ik}$
else if( ldafb<2_${ik}$*kl+ku+1 ) then
info = -10_${ik}$
else if( stdlib_lsame( fact, 'F' ) .and. .not.( rowequ .or. colequ .or. stdlib_lsame( &
equed, 'N' ) ) ) then
info = -12_${ik}$
else
if( rowequ ) then
rcmin = bignum
rcmax = zero
do j = 1, n
rcmin = min( rcmin, r( j ) )
rcmax = max( rcmax, r( j ) )
end do
if( rcmin<=zero ) then
info = -13_${ik}$
else if( n>0_${ik}$ ) then
rowcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
else
rowcnd = one
end if
end if
if( colequ .and. info==0_${ik}$ ) then
rcmin = bignum
rcmax = zero
do j = 1, n
rcmin = min( rcmin, c( j ) )
rcmax = max( rcmax, c( j ) )
end do
if( rcmin<=zero ) then
info = -14_${ik}$
else if( n>0_${ik}$ ) then
colcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
else
colcnd = one
end if
end if
if( info==0_${ik}$ ) then
if( ldb<max( 1_${ik}$, n ) ) then
info = -16_${ik}$
else if( ldx<max( 1_${ik}$, n ) ) then
info = -18_${ik}$
end if
end if
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'DGBSVX', -info )
return
end if
if( equil ) then
! compute row and column scalings to equilibrate the matrix a.
call stdlib${ii}$_${ri}$gbequ( n, n, kl, ku, ab, ldab, r, c, rowcnd, colcnd,amax, infequ )
if( infequ==0_${ik}$ ) then
! equilibrate the matrix.
call stdlib${ii}$_${ri}$laqgb( n, n, kl, ku, ab, ldab, r, c, rowcnd, colcnd,amax, equed )
rowequ = stdlib_lsame( equed, 'R' ) .or. stdlib_lsame( equed, 'B' )
colequ = stdlib_lsame( equed, 'C' ) .or. stdlib_lsame( equed, 'B' )
end if
end if
! scale the right hand side.
if( notran ) then
if( rowequ ) then
do j = 1, nrhs
do i = 1, n
b( i, j ) = r( i )*b( i, j )
end do
end do
end if
else if( colequ ) then
do j = 1, nrhs
do i = 1, n
b( i, j ) = c( i )*b( i, j )
end do
end do
end if
if( nofact .or. equil ) then
! compute the lu factorization of the band matrix a.
do j = 1, n
j1 = max( j-ku, 1_${ik}$ )
j2 = min( j+kl, n )
call stdlib${ii}$_${ri}$copy( j2-j1+1, ab( ku+1-j+j1, j ), 1_${ik}$,afb( kl+ku+1-j+j1, j ), 1_${ik}$ )
end do
call stdlib${ii}$_${ri}$gbtrf( n, n, kl, ku, afb, ldafb, ipiv, info )
! return if info is non-zero.
if( info>0_${ik}$ ) then
! compute the reciprocal pivot growth factor of the
! leading rank-deficient info columns of a.
anorm = zero
do j = 1, info
do i = max( ku+2-j, 1 ), min( n+ku+1-j, kl+ku+1 )
anorm = max( anorm, abs( ab( i, j ) ) )
end do
end do
rpvgrw = stdlib${ii}$_${ri}$lantb( 'M', 'U', 'N', info, min( info-1, kl+ku ),afb( max( 1_${ik}$, &
kl+ku+2-info ), 1_${ik}$ ), ldafb,work )
if( rpvgrw==zero ) then
rpvgrw = one
else
rpvgrw = anorm / rpvgrw
end if
work( 1_${ik}$ ) = rpvgrw
rcond = zero
return
end if
end if
! compute the norm of the matrix a and the
! reciprocal pivot growth factor rpvgrw.
if( notran ) then
norm = '1'
else
norm = 'I'
end if
anorm = stdlib${ii}$_${ri}$langb( norm, n, kl, ku, ab, ldab, work )
rpvgrw = stdlib${ii}$_${ri}$lantb( 'M', 'U', 'N', n, kl+ku, afb, ldafb, work )
if( rpvgrw==zero ) then
rpvgrw = one
else
rpvgrw = stdlib${ii}$_${ri}$langb( 'M', n, kl, ku, ab, ldab, work ) / rpvgrw
end if
! compute the reciprocal of the condition number of a.
call stdlib${ii}$_${ri}$gbcon( norm, n, kl, ku, afb, ldafb, ipiv, 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}$gbtrs( trans, n, kl, ku, nrhs, afb, ldafb, ipiv, 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}$gbrfs( trans, n, kl, ku, nrhs, ab, ldab, afb, ldafb, ipiv,b, ldb, x, ldx, &
ferr, berr, work, iwork, info )
! transform the solution matrix x to a solution of the original
! system.
if( notran ) then
if( colequ ) then
do j = 1, nrhs
do i = 1, n
x( i, j ) = c( i )*x( i, j )
end do
end do
do j = 1, nrhs
ferr( j ) = ferr( j ) / colcnd
end do
end if
else if( rowequ ) then
do j = 1, nrhs
do i = 1, n
x( i, j ) = r( i )*x( i, j )
end do
end do
do j = 1, nrhs
ferr( j ) = ferr( j ) / rowcnd
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}$
work( 1_${ik}$ ) = rpvgrw
return
end subroutine stdlib${ii}$_${ri}$gbsvx
#:endif
#:endfor
module subroutine stdlib${ii}$_cgbsvx( fact, trans, n, kl, ku, nrhs, ab, ldab, afb,ldafb, ipiv, equed, r, &
!! CGBSVX uses the LU factorization to compute the solution to a complex
!! system of linear equations A * X = B, A**T * X = B, or A**H * X = B,
!! where A is a band matrix of order N with KL subdiagonals and KU
!! superdiagonals, and X and B are N-by-NRHS matrices.
!! Error bounds on the solution and a condition estimate are also
!! provided.
c, b, ldb, x, 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, trans
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: kl, ku, ldab, ldafb, ldb, ldx, n, nrhs
real(sp), intent(out) :: rcond
! Array Arguments
integer(${ik}$), intent(inout) :: ipiv(*)
real(sp), intent(out) :: berr(*), ferr(*), rwork(*)
real(sp), intent(inout) :: c(*), r(*)
complex(sp), intent(inout) :: ab(ldab,*), afb(ldafb,*), b(ldb,*)
complex(sp), intent(out) :: work(*), x(ldx,*)
! =====================================================================
! moved setting of info = n+1 so info does not subsequently get
! overwritten. sven, 17 mar 05.
! =====================================================================
! Local Scalars
logical(lk) :: colequ, equil, nofact, notran, rowequ
character :: norm
integer(${ik}$) :: i, infequ, j, j1, j2
real(sp) :: amax, anorm, bignum, colcnd, rcmax, rcmin, rowcnd, rpvgrw, smlnum
! Intrinsic Functions
! Executable Statements
info = 0_${ik}$
nofact = stdlib_lsame( fact, 'N' )
equil = stdlib_lsame( fact, 'E' )
notran = stdlib_lsame( trans, 'N' )
if( nofact .or. equil ) then
equed = 'N'
rowequ = .false.
colequ = .false.
else
rowequ = stdlib_lsame( equed, 'R' ) .or. stdlib_lsame( equed, 'B' )
colequ = stdlib_lsame( equed, 'C' ) .or. stdlib_lsame( equed, 'B' )
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.notran .and. .not.stdlib_lsame( trans, 'T' ) .and. .not.stdlib_lsame( &
trans, 'C' ) ) then
info = -2_${ik}$
else if( n<0_${ik}$ ) then
info = -3_${ik}$
else if( kl<0_${ik}$ ) then
info = -4_${ik}$
else if( ku<0_${ik}$ ) then
info = -5_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -6_${ik}$
else if( ldab<kl+ku+1 ) then
info = -8_${ik}$
else if( ldafb<2_${ik}$*kl+ku+1 ) then
info = -10_${ik}$
else if( stdlib_lsame( fact, 'F' ) .and. .not.( rowequ .or. colequ .or. stdlib_lsame( &
equed, 'N' ) ) ) then
info = -12_${ik}$
else
if( rowequ ) then
rcmin = bignum
rcmax = zero
do j = 1, n
rcmin = min( rcmin, r( j ) )
rcmax = max( rcmax, r( j ) )
end do
if( rcmin<=zero ) then
info = -13_${ik}$
else if( n>0_${ik}$ ) then
rowcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
else
rowcnd = one
end if
end if
if( colequ .and. info==0_${ik}$ ) then
rcmin = bignum
rcmax = zero
do j = 1, n
rcmin = min( rcmin, c( j ) )
rcmax = max( rcmax, c( j ) )
end do
if( rcmin<=zero ) then
info = -14_${ik}$
else if( n>0_${ik}$ ) then
colcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
else
colcnd = one
end if
end if
if( info==0_${ik}$ ) then
if( ldb<max( 1_${ik}$, n ) ) then
info = -16_${ik}$
else if( ldx<max( 1_${ik}$, n ) ) then
info = -18_${ik}$
end if
end if
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'CGBSVX', -info )
return
end if
if( equil ) then
! compute row and column scalings to equilibrate the matrix a.
call stdlib${ii}$_cgbequ( n, n, kl, ku, ab, ldab, r, c, rowcnd, colcnd,amax, infequ )
if( infequ==0_${ik}$ ) then
! equilibrate the matrix.
call stdlib${ii}$_claqgb( n, n, kl, ku, ab, ldab, r, c, rowcnd, colcnd,amax, equed )
rowequ = stdlib_lsame( equed, 'R' ) .or. stdlib_lsame( equed, 'B' )
colequ = stdlib_lsame( equed, 'C' ) .or. stdlib_lsame( equed, 'B' )
end if
end if
! scale the right hand side.
if( notran ) then
if( rowequ ) then
do j = 1, nrhs
do i = 1, n
b( i, j ) = r( i )*b( i, j )
end do
end do
end if
else if( colequ ) then
do j = 1, nrhs
do i = 1, n
b( i, j ) = c( i )*b( i, j )
end do
end do
end if
if( nofact .or. equil ) then
! compute the lu factorization of the band matrix a.
do j = 1, n
j1 = max( j-ku, 1_${ik}$ )
j2 = min( j+kl, n )
call stdlib${ii}$_ccopy( j2-j1+1, ab( ku+1-j+j1, j ), 1_${ik}$,afb( kl+ku+1-j+j1, j ), 1_${ik}$ )
end do
call stdlib${ii}$_cgbtrf( n, n, kl, ku, afb, ldafb, ipiv, info )
! return if info is non-zero.
if( info>0_${ik}$ ) then
! compute the reciprocal pivot growth factor of the
! leading rank-deficient info columns of a.
anorm = zero
do j = 1, info
do i = max( ku+2-j, 1 ), min( n+ku+1-j, kl+ku+1 )
anorm = max( anorm, abs( ab( i, j ) ) )
end do
end do
rpvgrw = stdlib${ii}$_clantb( 'M', 'U', 'N', info, min( info-1, kl+ku ),afb( max( 1_${ik}$, &
kl+ku+2-info ), 1_${ik}$ ), ldafb,rwork )
if( rpvgrw==zero ) then
rpvgrw = one
else
rpvgrw = anorm / rpvgrw
end if
rwork( 1_${ik}$ ) = rpvgrw
rcond = zero
return
end if
end if
! compute the norm of the matrix a and the
! reciprocal pivot growth factor rpvgrw.
if( notran ) then
norm = '1'
else
norm = 'I'
end if
anorm = stdlib${ii}$_clangb( norm, n, kl, ku, ab, ldab, rwork )
rpvgrw = stdlib${ii}$_clantb( 'M', 'U', 'N', n, kl+ku, afb, ldafb, rwork )
if( rpvgrw==zero ) then
rpvgrw = one
else
rpvgrw = stdlib${ii}$_clangb( 'M', n, kl, ku, ab, ldab, rwork ) / rpvgrw
end if
! compute the reciprocal of the condition number of a.
call stdlib${ii}$_cgbcon( norm, n, kl, ku, afb, ldafb, ipiv, anorm, rcond,work, rwork, info )
! compute the solution matrix x.
call stdlib${ii}$_clacpy( 'FULL', n, nrhs, b, ldb, x, ldx )
call stdlib${ii}$_cgbtrs( trans, n, kl, ku, nrhs, afb, ldafb, ipiv, x, ldx,info )
! use iterative refinement to improve the computed solution and
! compute error bounds and backward error estimates for it.
call stdlib${ii}$_cgbrfs( trans, n, kl, ku, nrhs, ab, ldab, afb, ldafb, ipiv,b, ldb, x, ldx, &
ferr, berr, work, rwork, info )
! transform the solution matrix x to a solution of the original
! system.
if( notran ) then
if( colequ ) then
do j = 1, nrhs
do i = 1, n
x( i, j ) = c( i )*x( i, j )
end do
end do
do j = 1, nrhs
ferr( j ) = ferr( j ) / colcnd
end do
end if
else if( rowequ ) then
do j = 1, nrhs
do i = 1, n
x( i, j ) = r( i )*x( i, j )
end do
end do
do j = 1, nrhs
ferr( j ) = ferr( j ) / rowcnd
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}$
rwork( 1_${ik}$ ) = rpvgrw
return
end subroutine stdlib${ii}$_cgbsvx
module subroutine stdlib${ii}$_zgbsvx( fact, trans, n, kl, ku, nrhs, ab, ldab, afb,ldafb, ipiv, equed, r, &
!! ZGBSVX uses the LU factorization to compute the solution to a complex
!! system of linear equations A * X = B, A**T * X = B, or A**H * X = B,
!! where A is a band matrix of order N with KL subdiagonals and KU
!! superdiagonals, and X and B are N-by-NRHS matrices.
!! Error bounds on the solution and a condition estimate are also
!! provided.
c, b, ldb, x, 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, trans
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: kl, ku, ldab, ldafb, ldb, ldx, n, nrhs
real(dp), intent(out) :: rcond
! Array Arguments
integer(${ik}$), intent(inout) :: ipiv(*)
real(dp), intent(out) :: berr(*), ferr(*), rwork(*)
real(dp), intent(inout) :: c(*), r(*)
complex(dp), intent(inout) :: ab(ldab,*), afb(ldafb,*), b(ldb,*)
complex(dp), intent(out) :: work(*), x(ldx,*)
! =====================================================================
! moved setting of info = n+1 so info does not subsequently get
! overwritten. sven, 17 mar 05.
! =====================================================================
! Local Scalars
logical(lk) :: colequ, equil, nofact, notran, rowequ
character :: norm
integer(${ik}$) :: i, infequ, j, j1, j2
real(dp) :: amax, anorm, bignum, colcnd, rcmax, rcmin, rowcnd, rpvgrw, smlnum
! Intrinsic Functions
! Executable Statements
info = 0_${ik}$
nofact = stdlib_lsame( fact, 'N' )
equil = stdlib_lsame( fact, 'E' )
notran = stdlib_lsame( trans, 'N' )
if( nofact .or. equil ) then
equed = 'N'
rowequ = .false.
colequ = .false.
else
rowequ = stdlib_lsame( equed, 'R' ) .or. stdlib_lsame( equed, 'B' )
colequ = stdlib_lsame( equed, 'C' ) .or. stdlib_lsame( equed, 'B' )
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.notran .and. .not.stdlib_lsame( trans, 'T' ) .and. .not.stdlib_lsame( &
trans, 'C' ) ) then
info = -2_${ik}$
else if( n<0_${ik}$ ) then
info = -3_${ik}$
else if( kl<0_${ik}$ ) then
info = -4_${ik}$
else if( ku<0_${ik}$ ) then
info = -5_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -6_${ik}$
else if( ldab<kl+ku+1 ) then
info = -8_${ik}$
else if( ldafb<2_${ik}$*kl+ku+1 ) then
info = -10_${ik}$
else if( stdlib_lsame( fact, 'F' ) .and. .not.( rowequ .or. colequ .or. stdlib_lsame( &
equed, 'N' ) ) ) then
info = -12_${ik}$
else
if( rowequ ) then
rcmin = bignum
rcmax = zero
do j = 1, n
rcmin = min( rcmin, r( j ) )
rcmax = max( rcmax, r( j ) )
end do
if( rcmin<=zero ) then
info = -13_${ik}$
else if( n>0_${ik}$ ) then
rowcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
else
rowcnd = one
end if
end if
if( colequ .and. info==0_${ik}$ ) then
rcmin = bignum
rcmax = zero
do j = 1, n
rcmin = min( rcmin, c( j ) )
rcmax = max( rcmax, c( j ) )
end do
if( rcmin<=zero ) then
info = -14_${ik}$
else if( n>0_${ik}$ ) then
colcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
else
colcnd = one
end if
end if
if( info==0_${ik}$ ) then
if( ldb<max( 1_${ik}$, n ) ) then
info = -16_${ik}$
else if( ldx<max( 1_${ik}$, n ) ) then
info = -18_${ik}$
end if
end if
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'ZGBSVX', -info )
return
end if
if( equil ) then
! compute row and column scalings to equilibrate the matrix a.
call stdlib${ii}$_zgbequ( n, n, kl, ku, ab, ldab, r, c, rowcnd, colcnd,amax, infequ )
if( infequ==0_${ik}$ ) then
! equilibrate the matrix.
call stdlib${ii}$_zlaqgb( n, n, kl, ku, ab, ldab, r, c, rowcnd, colcnd,amax, equed )
rowequ = stdlib_lsame( equed, 'R' ) .or. stdlib_lsame( equed, 'B' )
colequ = stdlib_lsame( equed, 'C' ) .or. stdlib_lsame( equed, 'B' )
end if
end if
! scale the right hand side.
if( notran ) then
if( rowequ ) then
do j = 1, nrhs
do i = 1, n
b( i, j ) = r( i )*b( i, j )
end do
end do
end if
else if( colequ ) then
do j = 1, nrhs
do i = 1, n
b( i, j ) = c( i )*b( i, j )
end do
end do
end if
if( nofact .or. equil ) then
! compute the lu factorization of the band matrix a.
do j = 1, n
j1 = max( j-ku, 1_${ik}$ )
j2 = min( j+kl, n )
call stdlib${ii}$_zcopy( j2-j1+1, ab( ku+1-j+j1, j ), 1_${ik}$,afb( kl+ku+1-j+j1, j ), 1_${ik}$ )
end do
call stdlib${ii}$_zgbtrf( n, n, kl, ku, afb, ldafb, ipiv, info )
! return if info is non-zero.
if( info>0_${ik}$ ) then
! compute the reciprocal pivot growth factor of the
! leading rank-deficient info columns of a.
anorm = zero
do j = 1, info
do i = max( ku+2-j, 1 ), min( n+ku+1-j, kl+ku+1 )
anorm = max( anorm, abs( ab( i, j ) ) )
end do
end do
rpvgrw = stdlib${ii}$_zlantb( 'M', 'U', 'N', info, min( info-1, kl+ku ),afb( max( 1_${ik}$, &
kl+ku+2-info ), 1_${ik}$ ), ldafb,rwork )
if( rpvgrw==zero ) then
rpvgrw = one
else
rpvgrw = anorm / rpvgrw
end if
rwork( 1_${ik}$ ) = rpvgrw
rcond = zero
return
end if
end if
! compute the norm of the matrix a and the
! reciprocal pivot growth factor rpvgrw.
if( notran ) then
norm = '1'
else
norm = 'I'
end if
anorm = stdlib${ii}$_zlangb( norm, n, kl, ku, ab, ldab, rwork )
rpvgrw = stdlib${ii}$_zlantb( 'M', 'U', 'N', n, kl+ku, afb, ldafb, rwork )
if( rpvgrw==zero ) then
rpvgrw = one
else
rpvgrw = stdlib${ii}$_zlangb( 'M', n, kl, ku, ab, ldab, rwork ) / rpvgrw
end if
! compute the reciprocal of the condition number of a.
call stdlib${ii}$_zgbcon( norm, n, kl, ku, afb, ldafb, ipiv, anorm, rcond,work, rwork, info )
! compute the solution matrix x.
call stdlib${ii}$_zlacpy( 'FULL', n, nrhs, b, ldb, x, ldx )
call stdlib${ii}$_zgbtrs( trans, n, kl, ku, nrhs, afb, ldafb, ipiv, x, ldx,info )
! use iterative refinement to improve the computed solution and
! compute error bounds and backward error estimates for it.
call stdlib${ii}$_zgbrfs( trans, n, kl, ku, nrhs, ab, ldab, afb, ldafb, ipiv,b, ldb, x, ldx, &
ferr, berr, work, rwork, info )
! transform the solution matrix x to a solution of the original
! system.
if( notran ) then
if( colequ ) then
do j = 1, nrhs
do i = 1, n
x( i, j ) = c( i )*x( i, j )
end do
end do
do j = 1, nrhs
ferr( j ) = ferr( j ) / colcnd
end do
end if
else if( rowequ ) then
do j = 1, nrhs
do i = 1, n
x( i, j ) = r( i )*x( i, j )
end do
end do
do j = 1, nrhs
ferr( j ) = ferr( j ) / rowcnd
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}$
rwork( 1_${ik}$ ) = rpvgrw
return
end subroutine stdlib${ii}$_zgbsvx
#:for ck,ct,ci in CMPLX_KINDS_TYPES
#:if not ck in ["sp","dp"]
module subroutine stdlib${ii}$_${ci}$gbsvx( fact, trans, n, kl, ku, nrhs, ab, ldab, afb,ldafb, ipiv, equed, r, &
!! ZGBSVX: uses the LU factorization to compute the solution to a complex
!! system of linear equations A * X = B, A**T * X = B, or A**H * X = B,
!! where A is a band matrix of order N with KL subdiagonals and KU
!! superdiagonals, and X and B are N-by-NRHS matrices.
!! Error bounds on the solution and a condition estimate are also
!! provided.
c, b, ldb, x, 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, trans
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: kl, ku, ldab, ldafb, ldb, ldx, n, nrhs
real(${ck}$), intent(out) :: rcond
! Array Arguments
integer(${ik}$), intent(inout) :: ipiv(*)
real(${ck}$), intent(out) :: berr(*), ferr(*), rwork(*)
real(${ck}$), intent(inout) :: c(*), r(*)
complex(${ck}$), intent(inout) :: ab(ldab,*), afb(ldafb,*), b(ldb,*)
complex(${ck}$), intent(out) :: work(*), x(ldx,*)
! =====================================================================
! moved setting of info = n+1 so info does not subsequently get
! overwritten. sven, 17 mar 05.
! =====================================================================
! Local Scalars
logical(lk) :: colequ, equil, nofact, notran, rowequ
character :: norm
integer(${ik}$) :: i, infequ, j, j1, j2
real(${ck}$) :: amax, anorm, bignum, colcnd, rcmax, rcmin, rowcnd, rpvgrw, smlnum
! Intrinsic Functions
! Executable Statements
info = 0_${ik}$
nofact = stdlib_lsame( fact, 'N' )
equil = stdlib_lsame( fact, 'E' )
notran = stdlib_lsame( trans, 'N' )
if( nofact .or. equil ) then
equed = 'N'
rowequ = .false.
colequ = .false.
else
rowequ = stdlib_lsame( equed, 'R' ) .or. stdlib_lsame( equed, 'B' )
colequ = stdlib_lsame( equed, 'C' ) .or. stdlib_lsame( equed, 'B' )
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.notran .and. .not.stdlib_lsame( trans, 'T' ) .and. .not.stdlib_lsame( &
trans, 'C' ) ) then
info = -2_${ik}$
else if( n<0_${ik}$ ) then
info = -3_${ik}$
else if( kl<0_${ik}$ ) then
info = -4_${ik}$
else if( ku<0_${ik}$ ) then
info = -5_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -6_${ik}$
else if( ldab<kl+ku+1 ) then
info = -8_${ik}$
else if( ldafb<2_${ik}$*kl+ku+1 ) then
info = -10_${ik}$
else if( stdlib_lsame( fact, 'F' ) .and. .not.( rowequ .or. colequ .or. stdlib_lsame( &
equed, 'N' ) ) ) then
info = -12_${ik}$
else
if( rowequ ) then
rcmin = bignum
rcmax = zero
do j = 1, n
rcmin = min( rcmin, r( j ) )
rcmax = max( rcmax, r( j ) )
end do
if( rcmin<=zero ) then
info = -13_${ik}$
else if( n>0_${ik}$ ) then
rowcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
else
rowcnd = one
end if
end if
if( colequ .and. info==0_${ik}$ ) then
rcmin = bignum
rcmax = zero
do j = 1, n
rcmin = min( rcmin, c( j ) )
rcmax = max( rcmax, c( j ) )
end do
if( rcmin<=zero ) then
info = -14_${ik}$
else if( n>0_${ik}$ ) then
colcnd = max( rcmin, smlnum ) / min( rcmax, bignum )
else
colcnd = one
end if
end if
if( info==0_${ik}$ ) then
if( ldb<max( 1_${ik}$, n ) ) then
info = -16_${ik}$
else if( ldx<max( 1_${ik}$, n ) ) then
info = -18_${ik}$
end if
end if
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'ZGBSVX', -info )
return
end if
if( equil ) then
! compute row and column scalings to equilibrate the matrix a.
call stdlib${ii}$_${ci}$gbequ( n, n, kl, ku, ab, ldab, r, c, rowcnd, colcnd,amax, infequ )
if( infequ==0_${ik}$ ) then
! equilibrate the matrix.
call stdlib${ii}$_${ci}$laqgb( n, n, kl, ku, ab, ldab, r, c, rowcnd, colcnd,amax, equed )
rowequ = stdlib_lsame( equed, 'R' ) .or. stdlib_lsame( equed, 'B' )
colequ = stdlib_lsame( equed, 'C' ) .or. stdlib_lsame( equed, 'B' )
end if
end if
! scale the right hand side.
if( notran ) then
if( rowequ ) then
do j = 1, nrhs
do i = 1, n
b( i, j ) = r( i )*b( i, j )
end do
end do
end if
else if( colequ ) then
do j = 1, nrhs
do i = 1, n
b( i, j ) = c( i )*b( i, j )
end do
end do
end if
if( nofact .or. equil ) then
! compute the lu factorization of the band matrix a.
do j = 1, n
j1 = max( j-ku, 1_${ik}$ )
j2 = min( j+kl, n )
call stdlib${ii}$_${ci}$copy( j2-j1+1, ab( ku+1-j+j1, j ), 1_${ik}$,afb( kl+ku+1-j+j1, j ), 1_${ik}$ )
end do
call stdlib${ii}$_${ci}$gbtrf( n, n, kl, ku, afb, ldafb, ipiv, info )
! return if info is non-zero.
if( info>0_${ik}$ ) then
! compute the reciprocal pivot growth factor of the
! leading rank-deficient info columns of a.
anorm = zero
do j = 1, info
do i = max( ku+2-j, 1 ), min( n+ku+1-j, kl+ku+1 )
anorm = max( anorm, abs( ab( i, j ) ) )
end do
end do
rpvgrw = stdlib${ii}$_${ci}$lantb( 'M', 'U', 'N', info, min( info-1, kl+ku ),afb( max( 1_${ik}$, &
kl+ku+2-info ), 1_${ik}$ ), ldafb,rwork )
if( rpvgrw==zero ) then
rpvgrw = one
else
rpvgrw = anorm / rpvgrw
end if
rwork( 1_${ik}$ ) = rpvgrw
rcond = zero
return
end if
end if
! compute the norm of the matrix a and the
! reciprocal pivot growth factor rpvgrw.
if( notran ) then
norm = '1'
else
norm = 'I'
end if
anorm = stdlib${ii}$_${ci}$langb( norm, n, kl, ku, ab, ldab, rwork )
rpvgrw = stdlib${ii}$_${ci}$lantb( 'M', 'U', 'N', n, kl+ku, afb, ldafb, rwork )
if( rpvgrw==zero ) then
rpvgrw = one
else
rpvgrw = stdlib${ii}$_${ci}$langb( 'M', n, kl, ku, ab, ldab, rwork ) / rpvgrw
end if
! compute the reciprocal of the condition number of a.
call stdlib${ii}$_${ci}$gbcon( norm, n, kl, ku, afb, ldafb, ipiv, 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}$gbtrs( trans, n, kl, ku, nrhs, afb, ldafb, ipiv, 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}$gbrfs( trans, n, kl, ku, nrhs, ab, ldab, afb, ldafb, ipiv,b, ldb, x, ldx, &
ferr, berr, work, rwork, info )
! transform the solution matrix x to a solution of the original
! system.
if( notran ) then
if( colequ ) then
do j = 1, nrhs
do i = 1, n
x( i, j ) = c( i )*x( i, j )
end do
end do
do j = 1, nrhs
ferr( j ) = ferr( j ) / colcnd
end do
end if
else if( rowequ ) then
do j = 1, nrhs
do i = 1, n
x( i, j ) = r( i )*x( i, j )
end do
end do
do j = 1, nrhs
ferr( j ) = ferr( j ) / rowcnd
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}$
rwork( 1_${ik}$ ) = rpvgrw
return
end subroutine stdlib${ii}$_${ci}$gbsvx
#:endif
#:endfor
pure module subroutine stdlib${ii}$_sgtsv( n, nrhs, dl, d, du, b, ldb, info )
!! SGTSV solves the equation
!! A*X = B,
!! where A is an n by n tridiagonal matrix, by Gaussian elimination with
!! partial pivoting.
!! Note that the equation A**T*X = B may be solved by interchanging the
!! order of the arguments DU and DL.
! -- 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(*), dl(*), du(*)
! =====================================================================
! Local Scalars
integer(${ik}$) :: i, j
real(sp) :: fact, temp
! Intrinsic Functions
! Executable Statements
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 = -7_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'SGTSV ', -info )
return
end if
if( n==0 )return
if( nrhs==1_${ik}$ ) then
loop_10: do i = 1, n - 2
if( abs( d( i ) )>=abs( dl( i ) ) ) then
! no row interchange required
if( d( i )/=zero ) then
fact = dl( i ) / d( i )
d( i+1 ) = d( i+1 ) - fact*du( i )
b( i+1, 1_${ik}$ ) = b( i+1, 1_${ik}$ ) - fact*b( i, 1_${ik}$ )
else
info = i
return
end if
dl( i ) = zero
else
! interchange rows i and i+1
fact = d( i ) / dl( i )
d( i ) = dl( i )
temp = d( i+1 )
d( i+1 ) = du( i ) - fact*temp
dl( i ) = du( i+1 )
du( i+1 ) = -fact*dl( i )
du( i ) = temp
temp = b( i, 1_${ik}$ )
b( i, 1_${ik}$ ) = b( i+1, 1_${ik}$ )
b( i+1, 1_${ik}$ ) = temp - fact*b( i+1, 1_${ik}$ )
end if
end do loop_10
if( n>1_${ik}$ ) then
i = n - 1_${ik}$
if( abs( d( i ) )>=abs( dl( i ) ) ) then
if( d( i )/=zero ) then
fact = dl( i ) / d( i )
d( i+1 ) = d( i+1 ) - fact*du( i )
b( i+1, 1_${ik}$ ) = b( i+1, 1_${ik}$ ) - fact*b( i, 1_${ik}$ )
else
info = i
return
end if
else
fact = d( i ) / dl( i )
d( i ) = dl( i )
temp = d( i+1 )
d( i+1 ) = du( i ) - fact*temp
du( i ) = temp
temp = b( i, 1_${ik}$ )
b( i, 1_${ik}$ ) = b( i+1, 1_${ik}$ )
b( i+1, 1_${ik}$ ) = temp - fact*b( i+1, 1_${ik}$ )
end if
end if
if( d( n )==zero ) then
info = n
return
end if
else
loop_40: do i = 1, n - 2
if( abs( d( i ) )>=abs( dl( i ) ) ) then
! no row interchange required
if( d( i )/=zero ) then
fact = dl( i ) / d( i )
d( i+1 ) = d( i+1 ) - fact*du( i )
do j = 1, nrhs
b( i+1, j ) = b( i+1, j ) - fact*b( i, j )
end do
else
info = i
return
end if
dl( i ) = zero
else
! interchange rows i and i+1
fact = d( i ) / dl( i )
d( i ) = dl( i )
temp = d( i+1 )
d( i+1 ) = du( i ) - fact*temp
dl( i ) = du( i+1 )
du( i+1 ) = -fact*dl( i )
du( i ) = temp
do j = 1, nrhs
temp = b( i, j )
b( i, j ) = b( i+1, j )
b( i+1, j ) = temp - fact*b( i+1, j )
end do
end if
end do loop_40
if( n>1_${ik}$ ) then
i = n - 1_${ik}$
if( abs( d( i ) )>=abs( dl( i ) ) ) then
if( d( i )/=zero ) then
fact = dl( i ) / d( i )
d( i+1 ) = d( i+1 ) - fact*du( i )
do j = 1, nrhs
b( i+1, j ) = b( i+1, j ) - fact*b( i, j )
end do
else
info = i
return
end if
else
fact = d( i ) / dl( i )
d( i ) = dl( i )
temp = d( i+1 )
d( i+1 ) = du( i ) - fact*temp
du( i ) = temp
do j = 1, nrhs
temp = b( i, j )
b( i, j ) = b( i+1, j )
b( i+1, j ) = temp - fact*b( i+1, j )
end do
end if
end if
if( d( n )==zero ) then
info = n
return
end if
end if
! back solve with the matrix u from the factorization.
if( nrhs<=2_${ik}$ ) then
j = 1_${ik}$
70 continue
b( n, j ) = b( n, j ) / d( n )
if( n>1_${ik}$ )b( n-1, j ) = ( b( n-1, j )-du( n-1 )*b( n, j ) ) / d( n-1 )
do i = n - 2, 1, -1
b( i, j ) = ( b( i, j )-du( i )*b( i+1, j )-dl( i )*b( i+2, j ) ) / d( i )
end do
if( j<nrhs ) then
j = j + 1_${ik}$
go to 70
end if
else
do j = 1, nrhs
b( n, j ) = b( n, j ) / d( n )
if( n>1_${ik}$ )b( n-1, j ) = ( b( n-1, j )-du( n-1 )*b( n, j ) ) /d( n-1 )
do i = n - 2, 1, -1
b( i, j ) = ( b( i, j )-du( i )*b( i+1, j )-dl( i )*b( i+2, j ) ) / d( i )
end do
end do
end if
return
end subroutine stdlib${ii}$_sgtsv
pure module subroutine stdlib${ii}$_dgtsv( n, nrhs, dl, d, du, b, ldb, info )
!! DGTSV solves the equation
!! A*X = B,
!! where A is an n by n tridiagonal matrix, by Gaussian elimination with
!! partial pivoting.
!! Note that the equation A**T*X = B may be solved by interchanging the
!! order of the arguments DU and DL.
! -- 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(*), dl(*), du(*)
! =====================================================================
! Local Scalars
integer(${ik}$) :: i, j
real(dp) :: fact, temp
! Intrinsic Functions
! Executable Statements
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 = -7_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'DGTSV ', -info )
return
end if
if( n==0 )return
if( nrhs==1_${ik}$ ) then
loop_10: do i = 1, n - 2
if( abs( d( i ) )>=abs( dl( i ) ) ) then
! no row interchange required
if( d( i )/=zero ) then
fact = dl( i ) / d( i )
d( i+1 ) = d( i+1 ) - fact*du( i )
b( i+1, 1_${ik}$ ) = b( i+1, 1_${ik}$ ) - fact*b( i, 1_${ik}$ )
else
info = i
return
end if
dl( i ) = zero
else
! interchange rows i and i+1
fact = d( i ) / dl( i )
d( i ) = dl( i )
temp = d( i+1 )
d( i+1 ) = du( i ) - fact*temp
dl( i ) = du( i+1 )
du( i+1 ) = -fact*dl( i )
du( i ) = temp
temp = b( i, 1_${ik}$ )
b( i, 1_${ik}$ ) = b( i+1, 1_${ik}$ )
b( i+1, 1_${ik}$ ) = temp - fact*b( i+1, 1_${ik}$ )
end if
end do loop_10
if( n>1_${ik}$ ) then
i = n - 1_${ik}$
if( abs( d( i ) )>=abs( dl( i ) ) ) then
if( d( i )/=zero ) then
fact = dl( i ) / d( i )
d( i+1 ) = d( i+1 ) - fact*du( i )
b( i+1, 1_${ik}$ ) = b( i+1, 1_${ik}$ ) - fact*b( i, 1_${ik}$ )
else
info = i
return
end if
else
fact = d( i ) / dl( i )
d( i ) = dl( i )
temp = d( i+1 )
d( i+1 ) = du( i ) - fact*temp
du( i ) = temp
temp = b( i, 1_${ik}$ )
b( i, 1_${ik}$ ) = b( i+1, 1_${ik}$ )
b( i+1, 1_${ik}$ ) = temp - fact*b( i+1, 1_${ik}$ )
end if
end if
if( d( n )==zero ) then
info = n
return
end if
else
loop_40: do i = 1, n - 2
if( abs( d( i ) )>=abs( dl( i ) ) ) then
! no row interchange required
if( d( i )/=zero ) then
fact = dl( i ) / d( i )
d( i+1 ) = d( i+1 ) - fact*du( i )
do j = 1, nrhs
b( i+1, j ) = b( i+1, j ) - fact*b( i, j )
end do
else
info = i
return
end if
dl( i ) = zero
else
! interchange rows i and i+1
fact = d( i ) / dl( i )
d( i ) = dl( i )
temp = d( i+1 )
d( i+1 ) = du( i ) - fact*temp
dl( i ) = du( i+1 )
du( i+1 ) = -fact*dl( i )
du( i ) = temp
do j = 1, nrhs
temp = b( i, j )
b( i, j ) = b( i+1, j )
b( i+1, j ) = temp - fact*b( i+1, j )
end do
end if
end do loop_40
if( n>1_${ik}$ ) then
i = n - 1_${ik}$
if( abs( d( i ) )>=abs( dl( i ) ) ) then
if( d( i )/=zero ) then
fact = dl( i ) / d( i )
d( i+1 ) = d( i+1 ) - fact*du( i )
do j = 1, nrhs
b( i+1, j ) = b( i+1, j ) - fact*b( i, j )
end do
else
info = i
return
end if
else
fact = d( i ) / dl( i )
d( i ) = dl( i )
temp = d( i+1 )
d( i+1 ) = du( i ) - fact*temp
du( i ) = temp
do j = 1, nrhs
temp = b( i, j )
b( i, j ) = b( i+1, j )
b( i+1, j ) = temp - fact*b( i+1, j )
end do
end if
end if
if( d( n )==zero ) then
info = n
return
end if
end if
! back solve with the matrix u from the factorization.
if( nrhs<=2_${ik}$ ) then
j = 1_${ik}$
70 continue
b( n, j ) = b( n, j ) / d( n )
if( n>1_${ik}$ )b( n-1, j ) = ( b( n-1, j )-du( n-1 )*b( n, j ) ) / d( n-1 )
do i = n - 2, 1, -1
b( i, j ) = ( b( i, j )-du( i )*b( i+1, j )-dl( i )*b( i+2, j ) ) / d( i )
end do
if( j<nrhs ) then
j = j + 1_${ik}$
go to 70
end if
else
do j = 1, nrhs
b( n, j ) = b( n, j ) / d( n )
if( n>1_${ik}$ )b( n-1, j ) = ( b( n-1, j )-du( n-1 )*b( n, j ) ) /d( n-1 )
do i = n - 2, 1, -1
b( i, j ) = ( b( i, j )-du( i )*b( i+1, j )-dl( i )*b( i+2, j ) ) / d( i )
end do
end do
end if
return
end subroutine stdlib${ii}$_dgtsv
#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
pure module subroutine stdlib${ii}$_${ri}$gtsv( n, nrhs, dl, d, du, b, ldb, info )
!! DGTSV: solves the equation
!! A*X = B,
!! where A is an n by n tridiagonal matrix, by Gaussian elimination with
!! partial pivoting.
!! Note that the equation A**T*X = B may be solved by interchanging the
!! order of the arguments DU and DL.
! -- 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(*), dl(*), du(*)
! =====================================================================
! Local Scalars
integer(${ik}$) :: i, j
real(${rk}$) :: fact, temp
! Intrinsic Functions
! Executable Statements
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 = -7_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'DGTSV ', -info )
return
end if
if( n==0 )return
if( nrhs==1_${ik}$ ) then
loop_10: do i = 1, n - 2
if( abs( d( i ) )>=abs( dl( i ) ) ) then
! no row interchange required
if( d( i )/=zero ) then
fact = dl( i ) / d( i )
d( i+1 ) = d( i+1 ) - fact*du( i )
b( i+1, 1_${ik}$ ) = b( i+1, 1_${ik}$ ) - fact*b( i, 1_${ik}$ )
else
info = i
return
end if
dl( i ) = zero
else
! interchange rows i and i+1
fact = d( i ) / dl( i )
d( i ) = dl( i )
temp = d( i+1 )
d( i+1 ) = du( i ) - fact*temp
dl( i ) = du( i+1 )
du( i+1 ) = -fact*dl( i )
du( i ) = temp
temp = b( i, 1_${ik}$ )
b( i, 1_${ik}$ ) = b( i+1, 1_${ik}$ )
b( i+1, 1_${ik}$ ) = temp - fact*b( i+1, 1_${ik}$ )
end if
end do loop_10
if( n>1_${ik}$ ) then
i = n - 1_${ik}$
if( abs( d( i ) )>=abs( dl( i ) ) ) then
if( d( i )/=zero ) then
fact = dl( i ) / d( i )
d( i+1 ) = d( i+1 ) - fact*du( i )
b( i+1, 1_${ik}$ ) = b( i+1, 1_${ik}$ ) - fact*b( i, 1_${ik}$ )
else
info = i
return
end if
else
fact = d( i ) / dl( i )
d( i ) = dl( i )
temp = d( i+1 )
d( i+1 ) = du( i ) - fact*temp
du( i ) = temp
temp = b( i, 1_${ik}$ )
b( i, 1_${ik}$ ) = b( i+1, 1_${ik}$ )
b( i+1, 1_${ik}$ ) = temp - fact*b( i+1, 1_${ik}$ )
end if
end if
if( d( n )==zero ) then
info = n
return
end if
else
loop_40: do i = 1, n - 2
if( abs( d( i ) )>=abs( dl( i ) ) ) then
! no row interchange required
if( d( i )/=zero ) then
fact = dl( i ) / d( i )
d( i+1 ) = d( i+1 ) - fact*du( i )
do j = 1, nrhs
b( i+1, j ) = b( i+1, j ) - fact*b( i, j )
end do
else
info = i
return
end if
dl( i ) = zero
else
! interchange rows i and i+1
fact = d( i ) / dl( i )
d( i ) = dl( i )
temp = d( i+1 )
d( i+1 ) = du( i ) - fact*temp
dl( i ) = du( i+1 )
du( i+1 ) = -fact*dl( i )
du( i ) = temp
do j = 1, nrhs
temp = b( i, j )
b( i, j ) = b( i+1, j )
b( i+1, j ) = temp - fact*b( i+1, j )
end do
end if
end do loop_40
if( n>1_${ik}$ ) then
i = n - 1_${ik}$
if( abs( d( i ) )>=abs( dl( i ) ) ) then
if( d( i )/=zero ) then
fact = dl( i ) / d( i )
d( i+1 ) = d( i+1 ) - fact*du( i )
do j = 1, nrhs
b( i+1, j ) = b( i+1, j ) - fact*b( i, j )
end do
else
info = i
return
end if
else
fact = d( i ) / dl( i )
d( i ) = dl( i )
temp = d( i+1 )
d( i+1 ) = du( i ) - fact*temp
du( i ) = temp
do j = 1, nrhs
temp = b( i, j )
b( i, j ) = b( i+1, j )
b( i+1, j ) = temp - fact*b( i+1, j )
end do
end if
end if
if( d( n )==zero ) then
info = n
return
end if
end if
! back solve with the matrix u from the factorization.
if( nrhs<=2_${ik}$ ) then
j = 1_${ik}$
70 continue
b( n, j ) = b( n, j ) / d( n )
if( n>1_${ik}$ )b( n-1, j ) = ( b( n-1, j )-du( n-1 )*b( n, j ) ) / d( n-1 )
do i = n - 2, 1, -1
b( i, j ) = ( b( i, j )-du( i )*b( i+1, j )-dl( i )*b( i+2, j ) ) / d( i )
end do
if( j<nrhs ) then
j = j + 1_${ik}$
go to 70
end if
else
do j = 1, nrhs
b( n, j ) = b( n, j ) / d( n )
if( n>1_${ik}$ )b( n-1, j ) = ( b( n-1, j )-du( n-1 )*b( n, j ) ) /d( n-1 )
do i = n - 2, 1, -1
b( i, j ) = ( b( i, j )-du( i )*b( i+1, j )-dl( i )*b( i+2, j ) ) / d( i )
end do
end do
end if
return
end subroutine stdlib${ii}$_${ri}$gtsv
#:endif
#:endfor
pure module subroutine stdlib${ii}$_cgtsv( n, nrhs, dl, d, du, b, ldb, info )
!! CGTSV solves the equation
!! A*X = B,
!! where A is an N-by-N tridiagonal matrix, by Gaussian elimination with
!! partial pivoting.
!! Note that the equation A**T *X = B may be solved by interchanging the
!! order of the arguments DU and DL.
! -- 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
complex(sp), intent(inout) :: b(ldb,*), d(*), dl(*), du(*)
! =====================================================================
! Local Scalars
integer(${ik}$) :: j, k
complex(sp) :: mult, temp, zdum
! Intrinsic Functions
! Statement Functions
real(sp) :: cabs1
! Statement Function Definitions
cabs1( zdum ) = abs( real( zdum,KIND=sp) ) + abs( aimag( zdum ) )
! Executable Statements
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 = -7_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'CGTSV ', -info )
return
end if
if( n==0 )return
loop_30: do k = 1, n - 1
if( dl( k )==czero ) then
! subdiagonal is czero, no elimination is required.
if( d( k )==czero ) then
! diagonal is czero: set info = k and return; a unique
! solution can not be found.
info = k
return
end if
else if( cabs1( d( k ) )>=cabs1( dl( k ) ) ) then
! no row interchange required
mult = dl( k ) / d( k )
d( k+1 ) = d( k+1 ) - mult*du( k )
do j = 1, nrhs
b( k+1, j ) = b( k+1, j ) - mult*b( k, j )
end do
if( k<( n-1 ) )dl( k ) = czero
else
! interchange rows k and k+1
mult = d( k ) / dl( k )
d( k ) = dl( k )
temp = d( k+1 )
d( k+1 ) = du( k ) - mult*temp
if( k<( n-1 ) ) then
dl( k ) = du( k+1 )
du( k+1 ) = -mult*dl( k )
end if
du( k ) = temp
do j = 1, nrhs
temp = b( k, j )
b( k, j ) = b( k+1, j )
b( k+1, j ) = temp - mult*b( k+1, j )
end do
end if
end do loop_30
if( d( n )==czero ) then
info = n
return
end if
! back solve with the matrix u from the factorization.
do j = 1, nrhs
b( n, j ) = b( n, j ) / d( n )
if( n>1_${ik}$ )b( n-1, j ) = ( b( n-1, j )-du( n-1 )*b( n, j ) ) / d( n-1 )
do k = n - 2, 1, -1
b( k, j ) = ( b( k, j )-du( k )*b( k+1, j )-dl( k )*b( k+2, j ) ) / d( k )
end do
end do
return
end subroutine stdlib${ii}$_cgtsv
pure module subroutine stdlib${ii}$_zgtsv( n, nrhs, dl, d, du, b, ldb, info )
!! ZGTSV solves the equation
!! A*X = B,
!! where A is an N-by-N tridiagonal matrix, by Gaussian elimination with
!! partial pivoting.
!! Note that the equation A**T *X = B may be solved by interchanging the
!! order of the arguments DU and DL.
! -- 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
complex(dp), intent(inout) :: b(ldb,*), d(*), dl(*), du(*)
! =====================================================================
! Local Scalars
integer(${ik}$) :: j, k
complex(dp) :: mult, temp, zdum
! Intrinsic Functions
! Statement Functions
real(dp) :: cabs1
! Statement Function Definitions
cabs1( zdum ) = abs( real( zdum,KIND=dp) ) + abs( aimag( zdum ) )
! Executable Statements
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 = -7_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'ZGTSV ', -info )
return
end if
if( n==0 )return
loop_30: do k = 1, n - 1
if( dl( k )==czero ) then
! subdiagonal is czero, no elimination is required.
if( d( k )==czero ) then
! diagonal is czero: set info = k and return; a unique
! solution can not be found.
info = k
return
end if
else if( cabs1( d( k ) )>=cabs1( dl( k ) ) ) then
! no row interchange required
mult = dl( k ) / d( k )
d( k+1 ) = d( k+1 ) - mult*du( k )
do j = 1, nrhs
b( k+1, j ) = b( k+1, j ) - mult*b( k, j )
end do
if( k<( n-1 ) )dl( k ) = czero
else
! interchange rows k and k+1
mult = d( k ) / dl( k )
d( k ) = dl( k )
temp = d( k+1 )
d( k+1 ) = du( k ) - mult*temp
if( k<( n-1 ) ) then
dl( k ) = du( k+1 )
du( k+1 ) = -mult*dl( k )
end if
du( k ) = temp
do j = 1, nrhs
temp = b( k, j )
b( k, j ) = b( k+1, j )
b( k+1, j ) = temp - mult*b( k+1, j )
end do
end if
end do loop_30
if( d( n )==czero ) then
info = n
return
end if
! back solve with the matrix u from the factorization.
do j = 1, nrhs
b( n, j ) = b( n, j ) / d( n )
if( n>1_${ik}$ )b( n-1, j ) = ( b( n-1, j )-du( n-1 )*b( n, j ) ) / d( n-1 )
do k = n - 2, 1, -1
b( k, j ) = ( b( k, j )-du( k )*b( k+1, j )-dl( k )*b( k+2, j ) ) / d( k )
end do
end do
return
end subroutine stdlib${ii}$_zgtsv
#:for ck,ct,ci in CMPLX_KINDS_TYPES
#:if not ck in ["sp","dp"]
pure module subroutine stdlib${ii}$_${ci}$gtsv( n, nrhs, dl, d, du, b, ldb, info )
!! ZGTSV: solves the equation
!! A*X = B,
!! where A is an N-by-N tridiagonal matrix, by Gaussian elimination with
!! partial pivoting.
!! Note that the equation A**T *X = B may be solved by interchanging the
!! order of the arguments DU and DL.
! -- 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
complex(${ck}$), intent(inout) :: b(ldb,*), d(*), dl(*), du(*)
! =====================================================================
! Local Scalars
integer(${ik}$) :: j, k
complex(${ck}$) :: mult, temp, zdum
! Intrinsic Functions
! Statement Functions
real(${ck}$) :: cabs1
! Statement Function Definitions
cabs1( zdum ) = abs( real( zdum,KIND=${ck}$) ) + abs( aimag( zdum ) )
! Executable Statements
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 = -7_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'ZGTSV ', -info )
return
end if
if( n==0 )return
loop_30: do k = 1, n - 1
if( dl( k )==czero ) then
! subdiagonal is czero, no elimination is required.
if( d( k )==czero ) then
! diagonal is czero: set info = k and return; a unique
! solution can not be found.
info = k
return
end if
else if( cabs1( d( k ) )>=cabs1( dl( k ) ) ) then
! no row interchange required
mult = dl( k ) / d( k )
d( k+1 ) = d( k+1 ) - mult*du( k )
do j = 1, nrhs
b( k+1, j ) = b( k+1, j ) - mult*b( k, j )
end do
if( k<( n-1 ) )dl( k ) = czero
else
! interchange rows k and k+1
mult = d( k ) / dl( k )
d( k ) = dl( k )
temp = d( k+1 )
d( k+1 ) = du( k ) - mult*temp
if( k<( n-1 ) ) then
dl( k ) = du( k+1 )
du( k+1 ) = -mult*dl( k )
end if
du( k ) = temp
do j = 1, nrhs
temp = b( k, j )
b( k, j ) = b( k+1, j )
b( k+1, j ) = temp - mult*b( k+1, j )
end do
end if
end do loop_30
if( d( n )==czero ) then
info = n
return
end if
! back solve with the matrix u from the factorization.
do j = 1, nrhs
b( n, j ) = b( n, j ) / d( n )
if( n>1_${ik}$ )b( n-1, j ) = ( b( n-1, j )-du( n-1 )*b( n, j ) ) / d( n-1 )
do k = n - 2, 1, -1
b( k, j ) = ( b( k, j )-du( k )*b( k+1, j )-dl( k )*b( k+2, j ) ) / d( k )
end do
end do
return
end subroutine stdlib${ii}$_${ci}$gtsv
#:endif
#:endfor
pure module subroutine stdlib${ii}$_sgtsvx( fact, trans, n, nrhs, dl, d, du, dlf, df, duf,du2, ipiv, b, &
!! SGTSVX uses the LU factorization to compute the solution to a real
!! system of linear equations A * X = B or A**T * X = B,
!! where A is a tridiagonal matrix of order N and X and B are N-by-NRHS
!! matrices.
!! Error bounds on the solution and a condition estimate are also
!! provided.
ldb, x, 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(in) :: fact, trans
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: ldb, ldx, n, nrhs
real(sp), intent(out) :: rcond
! Array Arguments
integer(${ik}$), intent(inout) :: ipiv(*)
integer(${ik}$), intent(out) :: iwork(*)
real(sp), intent(in) :: b(ldb,*), d(*), dl(*), du(*)
real(sp), intent(out) :: berr(*), ferr(*), work(*), x(ldx,*)
real(sp), intent(inout) :: df(*), dlf(*), du2(*), duf(*)
! =====================================================================
! Local Scalars
logical(lk) :: nofact, notran
character :: norm
real(sp) :: anorm
! Intrinsic Functions
! Executable Statements
info = 0_${ik}$
nofact = stdlib_lsame( fact, 'N' )
notran = stdlib_lsame( trans, 'N' )
if( .not.nofact .and. .not.stdlib_lsame( fact, 'F' ) ) then
info = -1_${ik}$
else if( .not.notran .and. .not.stdlib_lsame( trans, 'T' ) .and. .not.stdlib_lsame( &
trans, 'C' ) ) then
info = -2_${ik}$
else if( n<0_${ik}$ ) then
info = -3_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -4_${ik}$
else if( ldb<max( 1_${ik}$, n ) ) then
info = -14_${ik}$
else if( ldx<max( 1_${ik}$, n ) ) then
info = -16_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'SGTSVX', -info )
return
end if
if( nofact ) then
! compute the lu factorization of a.
call stdlib${ii}$_scopy( n, d, 1_${ik}$, df, 1_${ik}$ )
if( n>1_${ik}$ ) then
call stdlib${ii}$_scopy( n-1, dl, 1_${ik}$, dlf, 1_${ik}$ )
call stdlib${ii}$_scopy( n-1, du, 1_${ik}$, duf, 1_${ik}$ )
end if
call stdlib${ii}$_sgttrf( n, dlf, df, duf, du2, ipiv, 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.
if( notran ) then
norm = '1'
else
norm = 'I'
end if
anorm = stdlib${ii}$_slangt( norm, n, dl, d, du )
! compute the reciprocal of the condition number of a.
call stdlib${ii}$_sgtcon( norm, n, dlf, df, duf, du2, ipiv, anorm, rcond, work,iwork, info )
! compute the solution vectors x.
call stdlib${ii}$_slacpy( 'FULL', n, nrhs, b, ldb, x, ldx )
call stdlib${ii}$_sgttrs( trans, n, nrhs, dlf, df, duf, du2, ipiv, x, ldx,info )
! use iterative refinement to improve the computed solutions and
! compute error bounds and backward error estimates for them.
call stdlib${ii}$_sgtrfs( trans, n, nrhs, dl, d, du, dlf, df, duf, du2, ipiv,b, ldb, x, ldx, &
ferr, berr, work, iwork, 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}$_sgtsvx
pure module subroutine stdlib${ii}$_dgtsvx( fact, trans, n, nrhs, dl, d, du, dlf, df, duf,du2, ipiv, b, &
!! DGTSVX uses the LU factorization to compute the solution to a real
!! system of linear equations A * X = B or A**T * X = B,
!! where A is a tridiagonal matrix of order N and X and B are N-by-NRHS
!! matrices.
!! Error bounds on the solution and a condition estimate are also
!! provided.
ldb, x, 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(in) :: fact, trans
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: ldb, ldx, n, nrhs
real(dp), intent(out) :: rcond
! Array Arguments
integer(${ik}$), intent(inout) :: ipiv(*)
integer(${ik}$), intent(out) :: iwork(*)
real(dp), intent(in) :: b(ldb,*), d(*), dl(*), du(*)
real(dp), intent(out) :: berr(*), ferr(*), work(*), x(ldx,*)
real(dp), intent(inout) :: df(*), dlf(*), du2(*), duf(*)
! =====================================================================
! Local Scalars
logical(lk) :: nofact, notran
character :: norm
real(dp) :: anorm
! Intrinsic Functions
! Executable Statements
info = 0_${ik}$
nofact = stdlib_lsame( fact, 'N' )
notran = stdlib_lsame( trans, 'N' )
if( .not.nofact .and. .not.stdlib_lsame( fact, 'F' ) ) then
info = -1_${ik}$
else if( .not.notran .and. .not.stdlib_lsame( trans, 'T' ) .and. .not.stdlib_lsame( &
trans, 'C' ) ) then
info = -2_${ik}$
else if( n<0_${ik}$ ) then
info = -3_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -4_${ik}$
else if( ldb<max( 1_${ik}$, n ) ) then
info = -14_${ik}$
else if( ldx<max( 1_${ik}$, n ) ) then
info = -16_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'DGTSVX', -info )
return
end if
if( nofact ) then
! compute the lu factorization of a.
call stdlib${ii}$_dcopy( n, d, 1_${ik}$, df, 1_${ik}$ )
if( n>1_${ik}$ ) then
call stdlib${ii}$_dcopy( n-1, dl, 1_${ik}$, dlf, 1_${ik}$ )
call stdlib${ii}$_dcopy( n-1, du, 1_${ik}$, duf, 1_${ik}$ )
end if
call stdlib${ii}$_dgttrf( n, dlf, df, duf, du2, ipiv, 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.
if( notran ) then
norm = '1'
else
norm = 'I'
end if
anorm = stdlib${ii}$_dlangt( norm, n, dl, d, du )
! compute the reciprocal of the condition number of a.
call stdlib${ii}$_dgtcon( norm, n, dlf, df, duf, du2, ipiv, anorm, rcond, work,iwork, info )
! compute the solution vectors x.
call stdlib${ii}$_dlacpy( 'FULL', n, nrhs, b, ldb, x, ldx )
call stdlib${ii}$_dgttrs( trans, n, nrhs, dlf, df, duf, du2, ipiv, x, ldx,info )
! use iterative refinement to improve the computed solutions and
! compute error bounds and backward error estimates for them.
call stdlib${ii}$_dgtrfs( trans, n, nrhs, dl, d, du, dlf, df, duf, du2, ipiv,b, ldb, x, ldx, &
ferr, berr, work, iwork, 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}$_dgtsvx
#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
pure module subroutine stdlib${ii}$_${ri}$gtsvx( fact, trans, n, nrhs, dl, d, du, dlf, df, duf,du2, ipiv, b, &
!! DGTSVX: uses the LU factorization to compute the solution to a real
!! system of linear equations A * X = B or A**T * X = B,
!! where A is a tridiagonal matrix of order N and X and B are N-by-NRHS
!! matrices.
!! Error bounds on the solution and a condition estimate are also
!! provided.
ldb, x, 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(in) :: fact, trans
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: ldb, ldx, n, nrhs
real(${rk}$), intent(out) :: rcond
! Array Arguments
integer(${ik}$), intent(inout) :: ipiv(*)
integer(${ik}$), intent(out) :: iwork(*)
real(${rk}$), intent(in) :: b(ldb,*), d(*), dl(*), du(*)
real(${rk}$), intent(out) :: berr(*), ferr(*), work(*), x(ldx,*)
real(${rk}$), intent(inout) :: df(*), dlf(*), du2(*), duf(*)
! =====================================================================
! Local Scalars
logical(lk) :: nofact, notran
character :: norm
real(${rk}$) :: anorm
! Intrinsic Functions
! Executable Statements
info = 0_${ik}$
nofact = stdlib_lsame( fact, 'N' )
notran = stdlib_lsame( trans, 'N' )
if( .not.nofact .and. .not.stdlib_lsame( fact, 'F' ) ) then
info = -1_${ik}$
else if( .not.notran .and. .not.stdlib_lsame( trans, 'T' ) .and. .not.stdlib_lsame( &
trans, 'C' ) ) then
info = -2_${ik}$
else if( n<0_${ik}$ ) then
info = -3_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -4_${ik}$
else if( ldb<max( 1_${ik}$, n ) ) then
info = -14_${ik}$
else if( ldx<max( 1_${ik}$, n ) ) then
info = -16_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'DGTSVX', -info )
return
end if
if( nofact ) then
! compute the lu factorization of a.
call stdlib${ii}$_${ri}$copy( n, d, 1_${ik}$, df, 1_${ik}$ )
if( n>1_${ik}$ ) then
call stdlib${ii}$_${ri}$copy( n-1, dl, 1_${ik}$, dlf, 1_${ik}$ )
call stdlib${ii}$_${ri}$copy( n-1, du, 1_${ik}$, duf, 1_${ik}$ )
end if
call stdlib${ii}$_${ri}$gttrf( n, dlf, df, duf, du2, ipiv, 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.
if( notran ) then
norm = '1'
else
norm = 'I'
end if
anorm = stdlib${ii}$_${ri}$langt( norm, n, dl, d, du )
! compute the reciprocal of the condition number of a.
call stdlib${ii}$_${ri}$gtcon( norm, n, dlf, df, duf, du2, ipiv, anorm, rcond, work,iwork, info )
! compute the solution vectors x.
call stdlib${ii}$_${ri}$lacpy( 'FULL', n, nrhs, b, ldb, x, ldx )
call stdlib${ii}$_${ri}$gttrs( trans, n, nrhs, dlf, df, duf, du2, ipiv, 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}$gtrfs( trans, n, nrhs, dl, d, du, dlf, df, duf, du2, ipiv,b, ldb, x, ldx, &
ferr, berr, work, iwork, 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}$gtsvx
#:endif
#:endfor
pure module subroutine stdlib${ii}$_cgtsvx( fact, trans, n, nrhs, dl, d, du, dlf, df, duf,du2, ipiv, b, &
!! CGTSVX uses the LU factorization to compute the solution to a complex
!! system of linear equations A * X = B, A**T * X = B, or A**H * X = B,
!! where A is a tridiagonal matrix of order N and X and B are N-by-NRHS
!! matrices.
!! Error bounds on the solution and a condition estimate are also
!! provided.
ldb, x, 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(in) :: fact, trans
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: ldb, ldx, n, nrhs
real(sp), intent(out) :: rcond
! Array Arguments
integer(${ik}$), intent(inout) :: ipiv(*)
real(sp), intent(out) :: berr(*), ferr(*), rwork(*)
complex(sp), intent(in) :: b(ldb,*), d(*), dl(*), du(*)
complex(sp), intent(inout) :: df(*), dlf(*), du2(*), duf(*)
complex(sp), intent(out) :: work(*), x(ldx,*)
! =====================================================================
! Local Scalars
logical(lk) :: nofact, notran
character :: norm
real(sp) :: anorm
! Intrinsic Functions
! Executable Statements
info = 0_${ik}$
nofact = stdlib_lsame( fact, 'N' )
notran = stdlib_lsame( trans, 'N' )
if( .not.nofact .and. .not.stdlib_lsame( fact, 'F' ) ) then
info = -1_${ik}$
else if( .not.notran .and. .not.stdlib_lsame( trans, 'T' ) .and. .not.stdlib_lsame( &
trans, 'C' ) ) then
info = -2_${ik}$
else if( n<0_${ik}$ ) then
info = -3_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -4_${ik}$
else if( ldb<max( 1_${ik}$, n ) ) then
info = -14_${ik}$
else if( ldx<max( 1_${ik}$, n ) ) then
info = -16_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'CGTSVX', -info )
return
end if
if( nofact ) then
! compute the lu factorization of a.
call stdlib${ii}$_ccopy( n, d, 1_${ik}$, df, 1_${ik}$ )
if( n>1_${ik}$ ) then
call stdlib${ii}$_ccopy( n-1, dl, 1_${ik}$, dlf, 1_${ik}$ )
call stdlib${ii}$_ccopy( n-1, du, 1_${ik}$, duf, 1_${ik}$ )
end if
call stdlib${ii}$_cgttrf( n, dlf, df, duf, du2, ipiv, 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.
if( notran ) then
norm = '1'
else
norm = 'I'
end if
anorm = stdlib${ii}$_clangt( norm, n, dl, d, du )
! compute the reciprocal of the condition number of a.
call stdlib${ii}$_cgtcon( norm, n, dlf, df, duf, du2, ipiv, anorm, rcond, work,info )
! compute the solution vectors x.
call stdlib${ii}$_clacpy( 'FULL', n, nrhs, b, ldb, x, ldx )
call stdlib${ii}$_cgttrs( trans, n, nrhs, dlf, df, duf, du2, ipiv, x, ldx,info )
! use iterative refinement to improve the computed solutions and
! compute error bounds and backward error estimates for them.
call stdlib${ii}$_cgtrfs( trans, n, nrhs, dl, d, du, dlf, df, duf, du2, ipiv,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}$_cgtsvx
pure module subroutine stdlib${ii}$_zgtsvx( fact, trans, n, nrhs, dl, d, du, dlf, df, duf,du2, ipiv, b, &
!! ZGTSVX uses the LU factorization to compute the solution to a complex
!! system of linear equations A * X = B, A**T * X = B, or A**H * X = B,
!! where A is a tridiagonal matrix of order N and X and B are N-by-NRHS
!! matrices.
!! Error bounds on the solution and a condition estimate are also
!! provided.
ldb, x, 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(in) :: fact, trans
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: ldb, ldx, n, nrhs
real(dp), intent(out) :: rcond
! Array Arguments
integer(${ik}$), intent(inout) :: ipiv(*)
real(dp), intent(out) :: berr(*), ferr(*), rwork(*)
complex(dp), intent(in) :: b(ldb,*), d(*), dl(*), du(*)
complex(dp), intent(inout) :: df(*), dlf(*), du2(*), duf(*)
complex(dp), intent(out) :: work(*), x(ldx,*)
! =====================================================================
! Local Scalars
logical(lk) :: nofact, notran
character :: norm
real(dp) :: anorm
! Intrinsic Functions
! Executable Statements
info = 0_${ik}$
nofact = stdlib_lsame( fact, 'N' )
notran = stdlib_lsame( trans, 'N' )
if( .not.nofact .and. .not.stdlib_lsame( fact, 'F' ) ) then
info = -1_${ik}$
else if( .not.notran .and. .not.stdlib_lsame( trans, 'T' ) .and. .not.stdlib_lsame( &
trans, 'C' ) ) then
info = -2_${ik}$
else if( n<0_${ik}$ ) then
info = -3_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -4_${ik}$
else if( ldb<max( 1_${ik}$, n ) ) then
info = -14_${ik}$
else if( ldx<max( 1_${ik}$, n ) ) then
info = -16_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'ZGTSVX', -info )
return
end if
if( nofact ) then
! compute the lu factorization of a.
call stdlib${ii}$_zcopy( n, d, 1_${ik}$, df, 1_${ik}$ )
if( n>1_${ik}$ ) then
call stdlib${ii}$_zcopy( n-1, dl, 1_${ik}$, dlf, 1_${ik}$ )
call stdlib${ii}$_zcopy( n-1, du, 1_${ik}$, duf, 1_${ik}$ )
end if
call stdlib${ii}$_zgttrf( n, dlf, df, duf, du2, ipiv, 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.
if( notran ) then
norm = '1'
else
norm = 'I'
end if
anorm = stdlib${ii}$_zlangt( norm, n, dl, d, du )
! compute the reciprocal of the condition number of a.
call stdlib${ii}$_zgtcon( norm, n, dlf, df, duf, du2, ipiv, anorm, rcond, work,info )
! compute the solution vectors x.
call stdlib${ii}$_zlacpy( 'FULL', n, nrhs, b, ldb, x, ldx )
call stdlib${ii}$_zgttrs( trans, n, nrhs, dlf, df, duf, du2, ipiv, x, ldx,info )
! use iterative refinement to improve the computed solutions and
! compute error bounds and backward error estimates for them.
call stdlib${ii}$_zgtrfs( trans, n, nrhs, dl, d, du, dlf, df, duf, du2, ipiv,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}$_zgtsvx
#:for ck,ct,ci in CMPLX_KINDS_TYPES
#:if not ck in ["sp","dp"]
pure module subroutine stdlib${ii}$_${ci}$gtsvx( fact, trans, n, nrhs, dl, d, du, dlf, df, duf,du2, ipiv, b, &
!! ZGTSVX: uses the LU factorization to compute the solution to a complex
!! system of linear equations A * X = B, A**T * X = B, or A**H * X = B,
!! where A is a tridiagonal matrix of order N and X and B are N-by-NRHS
!! matrices.
!! Error bounds on the solution and a condition estimate are also
!! provided.
ldb, x, 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(in) :: fact, trans
integer(${ik}$), intent(out) :: info
integer(${ik}$), intent(in) :: ldb, ldx, n, nrhs
real(${ck}$), intent(out) :: rcond
! Array Arguments
integer(${ik}$), intent(inout) :: ipiv(*)
real(${ck}$), intent(out) :: berr(*), ferr(*), rwork(*)
complex(${ck}$), intent(in) :: b(ldb,*), d(*), dl(*), du(*)
complex(${ck}$), intent(inout) :: df(*), dlf(*), du2(*), duf(*)
complex(${ck}$), intent(out) :: work(*), x(ldx,*)
! =====================================================================
! Local Scalars
logical(lk) :: nofact, notran
character :: norm
real(${ck}$) :: anorm
! Intrinsic Functions
! Executable Statements
info = 0_${ik}$
nofact = stdlib_lsame( fact, 'N' )
notran = stdlib_lsame( trans, 'N' )
if( .not.nofact .and. .not.stdlib_lsame( fact, 'F' ) ) then
info = -1_${ik}$
else if( .not.notran .and. .not.stdlib_lsame( trans, 'T' ) .and. .not.stdlib_lsame( &
trans, 'C' ) ) then
info = -2_${ik}$
else if( n<0_${ik}$ ) then
info = -3_${ik}$
else if( nrhs<0_${ik}$ ) then
info = -4_${ik}$
else if( ldb<max( 1_${ik}$, n ) ) then
info = -14_${ik}$
else if( ldx<max( 1_${ik}$, n ) ) then
info = -16_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'ZGTSVX', -info )
return
end if
if( nofact ) then
! compute the lu factorization of a.
call stdlib${ii}$_${ci}$copy( n, d, 1_${ik}$, df, 1_${ik}$ )
if( n>1_${ik}$ ) then
call stdlib${ii}$_${ci}$copy( n-1, dl, 1_${ik}$, dlf, 1_${ik}$ )
call stdlib${ii}$_${ci}$copy( n-1, du, 1_${ik}$, duf, 1_${ik}$ )
end if
call stdlib${ii}$_${ci}$gttrf( n, dlf, df, duf, du2, ipiv, 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.
if( notran ) then
norm = '1'
else
norm = 'I'
end if
anorm = stdlib${ii}$_${ci}$langt( norm, n, dl, d, du )
! compute the reciprocal of the condition number of a.
call stdlib${ii}$_${ci}$gtcon( norm, n, dlf, df, duf, du2, ipiv, anorm, rcond, work,info )
! compute the solution vectors x.
call stdlib${ii}$_${ci}$lacpy( 'FULL', n, nrhs, b, ldb, x, ldx )
call stdlib${ii}$_${ci}$gttrs( trans, n, nrhs, dlf, df, duf, du2, ipiv, 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}$gtrfs( trans, n, nrhs, dl, d, du, dlf, df, duf, du2, ipiv,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}$gtsvx
#:endif
#:endfor
#:endfor
end submodule stdlib_lapack_solve_lu
