#:include "common.fypp"
submodule(stdlib_lapack_others) stdlib_lapack_others_sm
implicit none
contains
#:for ik,it,ii in LINALG_INT_KINDS_TYPES
real(sp) module function stdlib${ii}$_sla_syrpvgrw( uplo, n, info, a, lda, af, ldaf, ipiv,work )
!! SLA_SYRPVGRW computes the reciprocal pivot growth factor
!! norm(A)/norm(U). The "max absolute element" norm is used. If this is
!! much less than 1, the stability of the LU factorization of the
!! (equilibrated) matrix A could be poor. This also means that the
!! solution X, estimated condition numbers, and error bounds could be
!! unreliable.
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: uplo
integer(${ik}$), intent(in) :: n, info, lda, ldaf
! Array Arguments
integer(${ik}$), intent(in) :: ipiv(*)
real(sp), intent(in) :: a(lda,*), af(ldaf,*)
real(sp), intent(out) :: work(*)
! =====================================================================
! Local Scalars
integer(${ik}$) :: ncols, i, j, k, kp
real(sp) :: amax, umax, rpvgrw, tmp
logical(lk) :: upper
! Intrinsic Functions
! Executable Statements
upper = stdlib_lsame( 'UPPER', uplo )
if ( info==0_${ik}$ ) then
if ( upper ) then
ncols = 1_${ik}$
else
ncols = n
end if
else
ncols = info
end if
rpvgrw = one
do i = 1, 2*n
work( i ) = zero
end do
! find the max magnitude entry of each column of a. compute the max
! for all n columns so we can apply the pivot permutation while
! looping below. assume a full factorization is the common case.
if ( upper ) then
do j = 1, n
do i = 1, j
work( n+i ) = max( abs( a( i, j ) ), work( n+i ) )
work( n+j ) = max( abs( a( i, j ) ), work( n+j ) )
end do
end do
else
do j = 1, n
do i = j, n
work( n+i ) = max( abs( a( i, j ) ), work( n+i ) )
work( n+j ) = max( abs( a( i, j ) ), work( n+j ) )
end do
end do
end if
! now find the max magnitude entry of each column of u or l. also
! permute the magnitudes of a above so they're in the same order as
! the factor.
! the iteration orders and permutations were copied from stdlib${ii}$_ssytrs.
! calls to stdlib${ii}$_sswap would be severe overkill.
if ( upper ) then
k = n
do while ( k < ncols .and. k>0 )
if ( ipiv( k )>0_${ik}$ ) then
! 1x1 pivot
kp = ipiv( k )
if ( kp /= k ) then
tmp = work( n+k )
work( n+k ) = work( n+kp )
work( n+kp ) = tmp
end if
do i = 1, k
work( k ) = max( abs( af( i, k ) ), work( k ) )
end do
k = k - 1_${ik}$
else
! 2x2 pivot
kp = -ipiv( k )
tmp = work( n+k-1 )
work( n+k-1 ) = work( n+kp )
work( n+kp ) = tmp
do i = 1, k-1
work( k ) = max( abs( af( i, k ) ), work( k ) )
work( k-1 ) = max( abs( af( i, k-1 ) ), work( k-1 ) )
end do
work( k ) = max( abs( af( k, k ) ), work( k ) )
k = k - 2_${ik}$
end if
end do
k = ncols
do while ( k <= n )
if ( ipiv( k )>0_${ik}$ ) then
kp = ipiv( k )
if ( kp /= k ) then
tmp = work( n+k )
work( n+k ) = work( n+kp )
work( n+kp ) = tmp
end if
k = k + 1_${ik}$
else
kp = -ipiv( k )
tmp = work( n+k )
work( n+k ) = work( n+kp )
work( n+kp ) = tmp
k = k + 2_${ik}$
end if
end do
else
k = 1_${ik}$
do while ( k <= ncols )
if ( ipiv( k )>0_${ik}$ ) then
! 1x1 pivot
kp = ipiv( k )
if ( kp /= k ) then
tmp = work( n+k )
work( n+k ) = work( n+kp )
work( n+kp ) = tmp
end if
do i = k, n
work( k ) = max( abs( af( i, k ) ), work( k ) )
end do
k = k + 1_${ik}$
else
! 2x2 pivot
kp = -ipiv( k )
tmp = work( n+k+1 )
work( n+k+1 ) = work( n+kp )
work( n+kp ) = tmp
do i = k+1, n
work( k ) = max( abs( af( i, k ) ), work( k ) )
work( k+1 ) = max( abs( af(i, k+1 ) ), work( k+1 ) )
end do
work( k ) = max( abs( af( k, k ) ), work( k ) )
k = k + 2_${ik}$
end if
end do
k = ncols
do while ( k >= 1 )
if ( ipiv( k )>0_${ik}$ ) then
kp = ipiv( k )
if ( kp /= k ) then
tmp = work( n+k )
work( n+k ) = work( n+kp )
work( n+kp ) = tmp
end if
k = k - 1_${ik}$
else
kp = -ipiv( k )
tmp = work( n+k )
work( n+k ) = work( n+kp )
work( n+kp ) = tmp
k = k - 2_${ik}$
endif
end do
end if
! compute the *inverse* of the max element growth factor. dividing
! by zero would imply the largest entry of the factor's column is
! zero. than can happen when either the column of a is zero or
! massive pivots made the factor underflow to zero. neither counts
! as growth in itself, so simply ignore terms with zero
! denominators.
if ( upper ) then
do i = ncols, n
umax = work( i )
amax = work( n+i )
if ( umax /= 0.0_sp ) then
rpvgrw = min( amax / umax, rpvgrw )
end if
end do
else
do i = 1, ncols
umax = work( i )
amax = work( n+i )
if ( umax /= 0.0_sp ) then
rpvgrw = min( amax / umax, rpvgrw )
end if
end do
end if
stdlib${ii}$_sla_syrpvgrw = rpvgrw
end function stdlib${ii}$_sla_syrpvgrw
real(dp) module function stdlib${ii}$_dla_syrpvgrw( uplo, n, info, a, lda, af,ldaf, ipiv, work )
!! DLA_SYRPVGRW computes the reciprocal pivot growth factor
!! norm(A)/norm(U). The "max absolute element" norm is used. If this is
!! much less than 1, the stability of the LU factorization of the
!! (equilibrated) matrix A could be poor. This also means that the
!! solution X, estimated condition numbers, and error bounds could be
!! unreliable.
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: uplo
integer(${ik}$), intent(in) :: n, info, lda, ldaf
! Array Arguments
integer(${ik}$), intent(in) :: ipiv(*)
real(dp), intent(in) :: a(lda,*), af(ldaf,*)
real(dp), intent(out) :: work(*)
! =====================================================================
! Local Scalars
integer(${ik}$) :: ncols, i, j, k, kp
real(dp) :: amax, umax, rpvgrw, tmp
logical(lk) :: upper
! Intrinsic Functions
! Executable Statements
upper = stdlib_lsame( 'UPPER', uplo )
if ( info==0_${ik}$ ) then
if ( upper ) then
ncols = 1_${ik}$
else
ncols = n
end if
else
ncols = info
end if
rpvgrw = one
do i = 1, 2*n
work( i ) = zero
end do
! find the max magnitude entry of each column of a. compute the max
! for all n columns so we can apply the pivot permutation while
! looping below. assume a full factorization is the common case.
if ( upper ) then
do j = 1, n
do i = 1, j
work( n+i ) = max( abs( a( i, j ) ), work( n+i ) )
work( n+j ) = max( abs( a( i, j ) ), work( n+j ) )
end do
end do
else
do j = 1, n
do i = j, n
work( n+i ) = max( abs( a( i, j ) ), work( n+i ) )
work( n+j ) = max( abs( a( i, j ) ), work( n+j ) )
end do
end do
end if
! now find the max magnitude entry of each column of u or l. also
! permute the magnitudes of a above so they're in the same order as
! the factor.
! the iteration orders and permutations were copied from stdlib${ii}$_dsytrs.
! calls to stdlib${ii}$_sswap would be severe overkill.
if ( upper ) then
k = n
do while ( k < ncols .and. k>0 )
if ( ipiv( k )>0_${ik}$ ) then
! 1x1 pivot
kp = ipiv( k )
if ( kp /= k ) then
tmp = work( n+k )
work( n+k ) = work( n+kp )
work( n+kp ) = tmp
end if
do i = 1, k
work( k ) = max( abs( af( i, k ) ), work( k ) )
end do
k = k - 1_${ik}$
else
! 2x2 pivot
kp = -ipiv( k )
tmp = work( n+k-1 )
work( n+k-1 ) = work( n+kp )
work( n+kp ) = tmp
do i = 1, k-1
work( k ) = max( abs( af( i, k ) ), work( k ) )
work( k-1 ) = max( abs( af( i, k-1 ) ), work( k-1 ) )
end do
work( k ) = max( abs( af( k, k ) ), work( k ) )
k = k - 2_${ik}$
end if
end do
k = ncols
do while ( k <= n )
if ( ipiv( k )>0_${ik}$ ) then
kp = ipiv( k )
if ( kp /= k ) then
tmp = work( n+k )
work( n+k ) = work( n+kp )
work( n+kp ) = tmp
end if
k = k + 1_${ik}$
else
kp = -ipiv( k )
tmp = work( n+k )
work( n+k ) = work( n+kp )
work( n+kp ) = tmp
k = k + 2_${ik}$
end if
end do
else
k = 1_${ik}$
do while ( k <= ncols )
if ( ipiv( k )>0_${ik}$ ) then
! 1x1 pivot
kp = ipiv( k )
if ( kp /= k ) then
tmp = work( n+k )
work( n+k ) = work( n+kp )
work( n+kp ) = tmp
end if
do i = k, n
work( k ) = max( abs( af( i, k ) ), work( k ) )
end do
k = k + 1_${ik}$
else
! 2x2 pivot
kp = -ipiv( k )
tmp = work( n+k+1 )
work( n+k+1 ) = work( n+kp )
work( n+kp ) = tmp
do i = k+1, n
work( k ) = max( abs( af( i, k ) ), work( k ) )
work( k+1 ) = max( abs( af(i, k+1 ) ), work( k+1 ) )
end do
work( k ) = max( abs( af( k, k ) ), work( k ) )
k = k + 2_${ik}$
end if
end do
k = ncols
do while ( k >= 1 )
if ( ipiv( k )>0_${ik}$ ) then
kp = ipiv( k )
if ( kp /= k ) then
tmp = work( n+k )
work( n+k ) = work( n+kp )
work( n+kp ) = tmp
end if
k = k - 1_${ik}$
else
kp = -ipiv( k )
tmp = work( n+k )
work( n+k ) = work( n+kp )
work( n+kp ) = tmp
k = k - 2_${ik}$
endif
end do
end if
! compute the *inverse* of the max element growth factor. dividing
! by zero would imply the largest entry of the factor's column is
! zero. than can happen when either the column of a is zero or
! massive pivots made the factor underflow to zero. neither counts
! as growth in itself, so simply ignore terms with zero
! denominators.
if ( upper ) then
do i = ncols, n
umax = work( i )
amax = work( n+i )
if ( umax /= zero ) then
rpvgrw = min( amax / umax, rpvgrw )
end if
end do
else
do i = 1, ncols
umax = work( i )
amax = work( n+i )
if ( umax /= zero ) then
rpvgrw = min( amax / umax, rpvgrw )
end if
end do
end if
stdlib${ii}$_dla_syrpvgrw = rpvgrw
end function stdlib${ii}$_dla_syrpvgrw
#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
real(${rk}$) module function stdlib${ii}$_${ri}$la_syrpvgrw( uplo, n, info, a, lda, af,ldaf, ipiv, work )
!! DLA_SYRPVGRW: computes the reciprocal pivot growth factor
!! norm(A)/norm(U). The "max absolute element" norm is used. If this is
!! much less than 1, the stability of the LU factorization of the
!! (equilibrated) matrix A could be poor. This also means that the
!! solution X, estimated condition numbers, and error bounds could be
!! unreliable.
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_${rk}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: uplo
integer(${ik}$), intent(in) :: n, info, lda, ldaf
! Array Arguments
integer(${ik}$), intent(in) :: ipiv(*)
real(${rk}$), intent(in) :: a(lda,*), af(ldaf,*)
real(${rk}$), intent(out) :: work(*)
! =====================================================================
! Local Scalars
integer(${ik}$) :: ncols, i, j, k, kp
real(${rk}$) :: amax, umax, rpvgrw, tmp
logical(lk) :: upper
! Intrinsic Functions
! Executable Statements
upper = stdlib_lsame( 'UPPER', uplo )
if ( info==0_${ik}$ ) then
if ( upper ) then
ncols = 1_${ik}$
else
ncols = n
end if
else
ncols = info
end if
rpvgrw = one
do i = 1, 2*n
work( i ) = zero
end do
! find the max magnitude entry of each column of a. compute the max
! for all n columns so we can apply the pivot permutation while
! looping below. assume a full factorization is the common case.
if ( upper ) then
do j = 1, n
do i = 1, j
work( n+i ) = max( abs( a( i, j ) ), work( n+i ) )
work( n+j ) = max( abs( a( i, j ) ), work( n+j ) )
end do
end do
else
do j = 1, n
do i = j, n
work( n+i ) = max( abs( a( i, j ) ), work( n+i ) )
work( n+j ) = max( abs( a( i, j ) ), work( n+j ) )
end do
end do
end if
! now find the max magnitude entry of each column of u or l. also
! permute the magnitudes of a above so they're in the same order as
! the factor.
! the iteration orders and permutations were copied from stdlib${ii}$_${ri}$sytrs.
! calls to stdlib${ii}$_dswap would be severe overkill.
if ( upper ) then
k = n
do while ( k < ncols .and. k>0 )
if ( ipiv( k )>0_${ik}$ ) then
! 1x1 pivot
kp = ipiv( k )
if ( kp /= k ) then
tmp = work( n+k )
work( n+k ) = work( n+kp )
work( n+kp ) = tmp
end if
do i = 1, k
work( k ) = max( abs( af( i, k ) ), work( k ) )
end do
k = k - 1_${ik}$
else
! 2x2 pivot
kp = -ipiv( k )
tmp = work( n+k-1 )
work( n+k-1 ) = work( n+kp )
work( n+kp ) = tmp
do i = 1, k-1
work( k ) = max( abs( af( i, k ) ), work( k ) )
work( k-1 ) = max( abs( af( i, k-1 ) ), work( k-1 ) )
end do
work( k ) = max( abs( af( k, k ) ), work( k ) )
k = k - 2_${ik}$
end if
end do
k = ncols
do while ( k <= n )
if ( ipiv( k )>0_${ik}$ ) then
kp = ipiv( k )
if ( kp /= k ) then
tmp = work( n+k )
work( n+k ) = work( n+kp )
work( n+kp ) = tmp
end if
k = k + 1_${ik}$
else
kp = -ipiv( k )
tmp = work( n+k )
work( n+k ) = work( n+kp )
work( n+kp ) = tmp
k = k + 2_${ik}$
end if
end do
else
k = 1_${ik}$
do while ( k <= ncols )
if ( ipiv( k )>0_${ik}$ ) then
! 1x1 pivot
kp = ipiv( k )
if ( kp /= k ) then
tmp = work( n+k )
work( n+k ) = work( n+kp )
work( n+kp ) = tmp
end if
do i = k, n
work( k ) = max( abs( af( i, k ) ), work( k ) )
end do
k = k + 1_${ik}$
else
! 2x2 pivot
kp = -ipiv( k )
tmp = work( n+k+1 )
work( n+k+1 ) = work( n+kp )
work( n+kp ) = tmp
do i = k+1, n
work( k ) = max( abs( af( i, k ) ), work( k ) )
work( k+1 ) = max( abs( af(i, k+1 ) ), work( k+1 ) )
end do
work( k ) = max( abs( af( k, k ) ), work( k ) )
k = k + 2_${ik}$
end if
end do
k = ncols
do while ( k >= 1 )
if ( ipiv( k )>0_${ik}$ ) then
kp = ipiv( k )
if ( kp /= k ) then
tmp = work( n+k )
work( n+k ) = work( n+kp )
work( n+kp ) = tmp
end if
k = k - 1_${ik}$
else
kp = -ipiv( k )
tmp = work( n+k )
work( n+k ) = work( n+kp )
work( n+kp ) = tmp
k = k - 2_${ik}$
endif
end do
end if
! compute the *inverse* of the max element growth factor. dividing
! by zero would imply the largest entry of the factor's column is
! zero. than can happen when either the column of a is zero or
! massive pivots made the factor underflow to zero. neither counts
! as growth in itself, so simply ignore terms with zero
! denominators.
if ( upper ) then
do i = ncols, n
umax = work( i )
amax = work( n+i )
if ( umax /= zero ) then
rpvgrw = min( amax / umax, rpvgrw )
end if
end do
else
do i = 1, ncols
umax = work( i )
amax = work( n+i )
if ( umax /= zero ) then
rpvgrw = min( amax / umax, rpvgrw )
end if
end do
end if
stdlib${ii}$_${ri}$la_syrpvgrw = rpvgrw
end function stdlib${ii}$_${ri}$la_syrpvgrw
#:endif
#:endfor
real(sp) module function stdlib${ii}$_cla_syrpvgrw( uplo, n, info, a, lda, af, ldaf, ipiv,work )
!! CLA_SYRPVGRW computes the reciprocal pivot growth factor
!! norm(A)/norm(U). The "max absolute element" norm is used. If this is
!! much less than 1, the stability of the LU factorization of the
!! (equilibrated) matrix A could be poor. This also means that the
!! solution X, estimated condition numbers, and error bounds could be
!! unreliable.
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: uplo
integer(${ik}$), intent(in) :: n, info, lda, ldaf
! Array Arguments
complex(sp), intent(in) :: a(lda,*), af(ldaf,*)
real(sp), intent(out) :: work(*)
integer(${ik}$), intent(in) :: ipiv(*)
! =====================================================================
! Local Scalars
integer(${ik}$) :: ncols, i, j, k, kp
real(sp) :: amax, umax, rpvgrw, tmp
logical(lk) :: upper
complex(sp) :: zdum
! Intrinsic Functions
! Statement Functions
real(sp) :: cabs1
! Statement Function Definitions
cabs1( zdum ) = abs( real( zdum,KIND=sp) ) + abs( aimag ( zdum ) )
! Executable Statements
upper = stdlib_lsame( 'UPPER', uplo )
if ( info==0_${ik}$ ) then
if ( upper ) then
ncols = 1_${ik}$
else
ncols = n
end if
else
ncols = info
end if
rpvgrw = one
do i = 1, 2*n
work( i ) = zero
end do
! find the max magnitude entry of each column of a. compute the max
! for all n columns so we can apply the pivot permutation while
! looping below. assume a full factorization is the common case.
if ( upper ) then
do j = 1, n
do i = 1, j
work( n+i ) = max( cabs1( a( i, j ) ), work( n+i ) )
work( n+j ) = max( cabs1( a( i, j ) ), work( n+j ) )
end do
end do
else
do j = 1, n
do i = j, n
work( n+i ) = max( cabs1( a( i, j ) ), work( n+i ) )
work( n+j ) = max( cabs1( a( i, j ) ), work( n+j ) )
end do
end do
end if
! now find the max magnitude entry of each column of u or l. also
! permute the magnitudes of a above so they're in the same order as
! the factor.
! the iteration orders and permutations were copied from stdlib${ii}$_csytrs.
! calls to stdlib${ii}$_sswap would be severe overkill.
if ( upper ) then
k = n
do while ( k < ncols .and. k>0 )
if ( ipiv( k )>0_${ik}$ ) then
! 1x1 pivot
kp = ipiv( k )
if ( kp /= k ) then
tmp = work( n+k )
work( n+k ) = work( n+kp )
work( n+kp ) = tmp
end if
do i = 1, k
work( k ) = max( cabs1( af( i, k ) ), work( k ) )
end do
k = k - 1_${ik}$
else
! 2x2 pivot
kp = -ipiv( k )
tmp = work( n+k-1 )
work( n+k-1 ) = work( n+kp )
work( n+kp ) = tmp
do i = 1, k-1
work( k ) = max( cabs1( af( i, k ) ), work( k ) )
work( k-1 ) =max( cabs1( af( i, k-1 ) ), work( k-1 ) )
end do
work( k ) = max( cabs1( af( k, k ) ), work( k ) )
k = k - 2_${ik}$
end if
end do
k = ncols
do while ( k <= n )
if ( ipiv( k )>0_${ik}$ ) then
kp = ipiv( k )
if ( kp /= k ) then
tmp = work( n+k )
work( n+k ) = work( n+kp )
work( n+kp ) = tmp
end if
k = k + 1_${ik}$
else
kp = -ipiv( k )
tmp = work( n+k )
work( n+k ) = work( n+kp )
work( n+kp ) = tmp
k = k + 2_${ik}$
end if
end do
else
k = 1_${ik}$
do while ( k <= ncols )
if ( ipiv( k )>0_${ik}$ ) then
! 1x1 pivot
kp = ipiv( k )
if ( kp /= k ) then
tmp = work( n+k )
work( n+k ) = work( n+kp )
work( n+kp ) = tmp
end if
do i = k, n
work( k ) = max( cabs1( af( i, k ) ), work( k ) )
end do
k = k + 1_${ik}$
else
! 2x2 pivot
kp = -ipiv( k )
tmp = work( n+k+1 )
work( n+k+1 ) = work( n+kp )
work( n+kp ) = tmp
do i = k+1, n
work( k ) = max( cabs1( af( i, k ) ), work( k ) )
work( k+1 ) =max( cabs1( af( i, k+1 ) ), work( k+1 ) )
end do
work( k ) = max( cabs1( af( k, k ) ), work( k ) )
k = k + 2_${ik}$
end if
end do
k = ncols
do while ( k >= 1 )
if ( ipiv( k )>0_${ik}$ ) then
kp = ipiv( k )
if ( kp /= k ) then
tmp = work( n+k )
work( n+k ) = work( n+kp )
work( n+kp ) = tmp
end if
k = k - 1_${ik}$
else
kp = -ipiv( k )
tmp = work( n+k )
work( n+k ) = work( n+kp )
work( n+kp ) = tmp
k = k - 2_${ik}$
endif
end do
end if
! compute the *inverse* of the max element growth factor. dividing
! by zero would imply the largest entry of the factor's column is
! zero. than can happen when either the column of a is zero or
! massive pivots made the factor underflow to zero. neither counts
! as growth in itself, so simply ignore terms with zero
! denominators.
if ( upper ) then
do i = ncols, n
umax = work( i )
amax = work( n+i )
if ( umax /= 0.0_sp ) then
rpvgrw = min( amax / umax, rpvgrw )
end if
end do
else
do i = 1, ncols
umax = work( i )
amax = work( n+i )
if ( umax /= 0.0_sp ) then
rpvgrw = min( amax / umax, rpvgrw )
end if
end do
end if
stdlib${ii}$_cla_syrpvgrw = rpvgrw
end function stdlib${ii}$_cla_syrpvgrw
real(dp) module function stdlib${ii}$_zla_syrpvgrw( uplo, n, info, a, lda, af,ldaf, ipiv, work )
!! ZLA_SYRPVGRW computes the reciprocal pivot growth factor
!! norm(A)/norm(U). The "max absolute element" norm is used. If this is
!! much less than 1, the stability of the LU factorization of the
!! (equilibrated) matrix A could be poor. This also means that the
!! solution X, estimated condition numbers, and error bounds could be
!! unreliable.
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: uplo
integer(${ik}$), intent(in) :: n, info, lda, ldaf
! Array Arguments
complex(dp), intent(in) :: a(lda,*), af(ldaf,*)
real(dp), intent(out) :: work(*)
integer(${ik}$), intent(in) :: ipiv(*)
! =====================================================================
! Local Scalars
integer(${ik}$) :: ncols, i, j, k, kp
real(dp) :: amax, umax, rpvgrw, tmp
logical(lk) :: upper
complex(dp) :: zdum
! Intrinsic Functions
! Statement Functions
real(dp) :: cabs1
! Statement Function Definitions
cabs1( zdum ) = abs( real( zdum,KIND=dp) ) + abs( aimag ( zdum ) )
! Executable Statements
upper = stdlib_lsame( 'UPPER', uplo )
if ( info==0_${ik}$ ) then
if ( upper ) then
ncols = 1_${ik}$
else
ncols = n
end if
else
ncols = info
end if
rpvgrw = one
do i = 1, 2*n
work( i ) = zero
end do
! find the max magnitude entry of each column of a. compute the max
! for all n columns so we can apply the pivot permutation while
! looping below. assume a full factorization is the common case.
if ( upper ) then
do j = 1, n
do i = 1, j
work( n+i ) = max( cabs1( a( i, j ) ), work( n+i ) )
work( n+j ) = max( cabs1( a( i, j ) ), work( n+j ) )
end do
end do
else
do j = 1, n
do i = j, n
work( n+i ) = max( cabs1( a( i, j ) ), work( n+i ) )
work( n+j ) = max( cabs1( a( i, j ) ), work( n+j ) )
end do
end do
end if
! now find the max magnitude entry of each column of u or l. also
! permute the magnitudes of a above so they're in the same order as
! the factor.
! the iteration orders and permutations were copied from stdlib${ii}$_zsytrs.
! calls to stdlib${ii}$_sswap would be severe overkill.
if ( upper ) then
k = n
do while ( k < ncols .and. k>0 )
if ( ipiv( k )>0_${ik}$ ) then
! 1x1 pivot
kp = ipiv( k )
if ( kp /= k ) then
tmp = work( n+k )
work( n+k ) = work( n+kp )
work( n+kp ) = tmp
end if
do i = 1, k
work( k ) = max( cabs1( af( i, k ) ), work( k ) )
end do
k = k - 1_${ik}$
else
! 2x2 pivot
kp = -ipiv( k )
tmp = work( n+k-1 )
work( n+k-1 ) = work( n+kp )
work( n+kp ) = tmp
do i = 1, k-1
work( k ) = max( cabs1( af( i, k ) ), work( k ) )
work( k-1 ) =max( cabs1( af( i, k-1 ) ), work( k-1 ) )
end do
work( k ) = max( cabs1( af( k, k ) ), work( k ) )
k = k - 2_${ik}$
end if
end do
k = ncols
do while ( k <= n )
if ( ipiv( k )>0_${ik}$ ) then
kp = ipiv( k )
if ( kp /= k ) then
tmp = work( n+k )
work( n+k ) = work( n+kp )
work( n+kp ) = tmp
end if
k = k + 1_${ik}$
else
kp = -ipiv( k )
tmp = work( n+k )
work( n+k ) = work( n+kp )
work( n+kp ) = tmp
k = k + 2_${ik}$
end if
end do
else
k = 1_${ik}$
do while ( k <= ncols )
if ( ipiv( k )>0_${ik}$ ) then
! 1x1 pivot
kp = ipiv( k )
if ( kp /= k ) then
tmp = work( n+k )
work( n+k ) = work( n+kp )
work( n+kp ) = tmp
end if
do i = k, n
work( k ) = max( cabs1( af( i, k ) ), work( k ) )
end do
k = k + 1_${ik}$
else
! 2x2 pivot
kp = -ipiv( k )
tmp = work( n+k+1 )
work( n+k+1 ) = work( n+kp )
work( n+kp ) = tmp
do i = k+1, n
work( k ) = max( cabs1( af( i, k ) ), work( k ) )
work( k+1 ) =max( cabs1( af( i, k+1 ) ), work( k+1 ) )
end do
work( k ) = max( cabs1( af( k, k ) ), work( k ) )
k = k + 2_${ik}$
end if
end do
k = ncols
do while ( k >= 1 )
if ( ipiv( k )>0_${ik}$ ) then
kp = ipiv( k )
if ( kp /= k ) then
tmp = work( n+k )
work( n+k ) = work( n+kp )
work( n+kp ) = tmp
end if
k = k - 1_${ik}$
else
kp = -ipiv( k )
tmp = work( n+k )
work( n+k ) = work( n+kp )
work( n+kp ) = tmp
k = k - 2_${ik}$
endif
end do
end if
! compute the *inverse* of the max element growth factor. dividing
! by zero would imply the largest entry of the factor's column is
! zero. than can happen when either the column of a is zero or
! massive pivots made the factor underflow to zero. neither counts
! as growth in itself, so simply ignore terms with zero
! denominators.
if ( upper ) then
do i = ncols, n
umax = work( i )
amax = work( n+i )
if ( umax /= zero ) then
rpvgrw = min( amax / umax, rpvgrw )
end if
end do
else
do i = 1, ncols
umax = work( i )
amax = work( n+i )
if ( umax /= zero ) then
rpvgrw = min( amax / umax, rpvgrw )
end if
end do
end if
stdlib${ii}$_zla_syrpvgrw = rpvgrw
end function stdlib${ii}$_zla_syrpvgrw
#:for ck,ct,ci in CMPLX_KINDS_TYPES
#:if not ck in ["sp","dp"]
real(${ck}$) module function stdlib${ii}$_${ci}$la_syrpvgrw( uplo, n, info, a, lda, af,ldaf, ipiv, work )
!! ZLA_SYRPVGRW: computes the reciprocal pivot growth factor
!! norm(A)/norm(U). The "max absolute element" norm is used. If this is
!! much less than 1, the stability of the LU factorization of the
!! (equilibrated) matrix A could be poor. This also means that the
!! solution X, estimated condition numbers, and error bounds could be
!! unreliable.
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_${ck}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: uplo
integer(${ik}$), intent(in) :: n, info, lda, ldaf
! Array Arguments
complex(${ck}$), intent(in) :: a(lda,*), af(ldaf,*)
real(${ck}$), intent(out) :: work(*)
integer(${ik}$), intent(in) :: ipiv(*)
! =====================================================================
! Local Scalars
integer(${ik}$) :: ncols, i, j, k, kp
real(${ck}$) :: amax, umax, rpvgrw, tmp
logical(lk) :: upper
complex(${ck}$) :: zdum
! Intrinsic Functions
! Statement Functions
real(${ck}$) :: cabs1
! Statement Function Definitions
cabs1( zdum ) = abs( real( zdum,KIND=${ck}$) ) + abs( aimag ( zdum ) )
! Executable Statements
upper = stdlib_lsame( 'UPPER', uplo )
if ( info==0_${ik}$ ) then
if ( upper ) then
ncols = 1_${ik}$
else
ncols = n
end if
else
ncols = info
end if
rpvgrw = one
do i = 1, 2*n
work( i ) = zero
end do
! find the max magnitude entry of each column of a. compute the max
! for all n columns so we can apply the pivot permutation while
! looping below. assume a full factorization is the common case.
if ( upper ) then
do j = 1, n
do i = 1, j
work( n+i ) = max( cabs1( a( i, j ) ), work( n+i ) )
work( n+j ) = max( cabs1( a( i, j ) ), work( n+j ) )
end do
end do
else
do j = 1, n
do i = j, n
work( n+i ) = max( cabs1( a( i, j ) ), work( n+i ) )
work( n+j ) = max( cabs1( a( i, j ) ), work( n+j ) )
end do
end do
end if
! now find the max magnitude entry of each column of u or l. also
! permute the magnitudes of a above so they're in the same order as
! the factor.
! the iteration orders and permutations were copied from stdlib${ii}$_${ci}$sytrs.
! calls to stdlib${ii}$_dswap would be severe overkill.
if ( upper ) then
k = n
do while ( k < ncols .and. k>0 )
if ( ipiv( k )>0_${ik}$ ) then
! 1x1 pivot
kp = ipiv( k )
if ( kp /= k ) then
tmp = work( n+k )
work( n+k ) = work( n+kp )
work( n+kp ) = tmp
end if
do i = 1, k
work( k ) = max( cabs1( af( i, k ) ), work( k ) )
end do
k = k - 1_${ik}$
else
! 2x2 pivot
kp = -ipiv( k )
tmp = work( n+k-1 )
work( n+k-1 ) = work( n+kp )
work( n+kp ) = tmp
do i = 1, k-1
work( k ) = max( cabs1( af( i, k ) ), work( k ) )
work( k-1 ) =max( cabs1( af( i, k-1 ) ), work( k-1 ) )
end do
work( k ) = max( cabs1( af( k, k ) ), work( k ) )
k = k - 2_${ik}$
end if
end do
k = ncols
do while ( k <= n )
if ( ipiv( k )>0_${ik}$ ) then
kp = ipiv( k )
if ( kp /= k ) then
tmp = work( n+k )
work( n+k ) = work( n+kp )
work( n+kp ) = tmp
end if
k = k + 1_${ik}$
else
kp = -ipiv( k )
tmp = work( n+k )
work( n+k ) = work( n+kp )
work( n+kp ) = tmp
k = k + 2_${ik}$
end if
end do
else
k = 1_${ik}$
do while ( k <= ncols )
if ( ipiv( k )>0_${ik}$ ) then
! 1x1 pivot
kp = ipiv( k )
if ( kp /= k ) then
tmp = work( n+k )
work( n+k ) = work( n+kp )
work( n+kp ) = tmp
end if
do i = k, n
work( k ) = max( cabs1( af( i, k ) ), work( k ) )
end do
k = k + 1_${ik}$
else
! 2x2 pivot
kp = -ipiv( k )
tmp = work( n+k+1 )
work( n+k+1 ) = work( n+kp )
work( n+kp ) = tmp
do i = k+1, n
work( k ) = max( cabs1( af( i, k ) ), work( k ) )
work( k+1 ) =max( cabs1( af( i, k+1 ) ), work( k+1 ) )
end do
work( k ) = max( cabs1( af( k, k ) ), work( k ) )
k = k + 2_${ik}$
end if
end do
k = ncols
do while ( k >= 1 )
if ( ipiv( k )>0_${ik}$ ) then
kp = ipiv( k )
if ( kp /= k ) then
tmp = work( n+k )
work( n+k ) = work( n+kp )
work( n+kp ) = tmp
end if
k = k - 1_${ik}$
else
kp = -ipiv( k )
tmp = work( n+k )
work( n+k ) = work( n+kp )
work( n+kp ) = tmp
k = k - 2_${ik}$
endif
end do
end if
! compute the *inverse* of the max element growth factor. dividing
! by zero would imply the largest entry of the factor's column is
! zero. than can happen when either the column of a is zero or
! massive pivots made the factor underflow to zero. neither counts
! as growth in itself, so simply ignore terms with zero
! denominators.
if ( upper ) then
do i = ncols, n
umax = work( i )
amax = work( n+i )
if ( umax /= zero ) then
rpvgrw = min( amax / umax, rpvgrw )
end if
end do
else
do i = 1, ncols
umax = work( i )
amax = work( n+i )
if ( umax /= zero ) then
rpvgrw = min( amax / umax, rpvgrw )
end if
end do
end if
stdlib${ii}$_${ci}$la_syrpvgrw = rpvgrw
end function stdlib${ii}$_${ci}$la_syrpvgrw
#:endif
#:endfor
pure real(sp) module function stdlib${ii}$_sla_gerpvgrw( n, ncols, a, lda, af, ldaf )
!! SLA_GERPVGRW computes the reciprocal pivot growth factor
!! norm(A)/norm(U). The "max absolute element" norm is used. If this is
!! much less than 1, the stability of the LU factorization of the
!! (equilibrated) matrix A could be poor. This also means that the
!! solution X, estimated condition numbers, and error bounds could be
!! unreliable.
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
integer(${ik}$), intent(in) :: n, ncols, lda, ldaf
! Array Arguments
real(sp), intent(in) :: a(lda,*), af(ldaf,*)
! =====================================================================
! Local Scalars
integer(${ik}$) :: i, j
real(sp) :: amax, umax, rpvgrw
! Intrinsic Functions
! Executable Statements
rpvgrw = one
do j = 1, ncols
amax = zero
umax = zero
do i = 1, n
amax = max( abs( a( i, j ) ), amax )
end do
do i = 1, j
umax = max( abs( af( i, j ) ), umax )
end do
if ( umax /= 0.0_sp ) then
rpvgrw = min( amax / umax, rpvgrw )
end if
end do
stdlib${ii}$_sla_gerpvgrw = rpvgrw
end function stdlib${ii}$_sla_gerpvgrw
pure real(dp) module function stdlib${ii}$_dla_gerpvgrw( n, ncols, a, lda, af,ldaf )
!! DLA_GERPVGRW computes the reciprocal pivot growth factor
!! norm(A)/norm(U). The "max absolute element" norm is used. If this is
!! much less than 1, the stability of the LU factorization of the
!! (equilibrated) matrix A could be poor. This also means that the
!! solution X, estimated condition numbers, and error bounds could be
!! unreliable.
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
integer(${ik}$), intent(in) :: n, ncols, lda, ldaf
! Array Arguments
real(dp), intent(in) :: a(lda,*), af(ldaf,*)
! =====================================================================
! Local Scalars
integer(${ik}$) :: i, j
real(dp) :: amax, umax, rpvgrw
! Intrinsic Functions
! Executable Statements
rpvgrw = one
do j = 1, ncols
amax = zero
umax = zero
do i = 1, n
amax = max( abs( a( i, j ) ), amax )
end do
do i = 1, j
umax = max( abs( af( i, j ) ), umax )
end do
if ( umax /= zero ) then
rpvgrw = min( amax / umax, rpvgrw )
end if
end do
stdlib${ii}$_dla_gerpvgrw = rpvgrw
end function stdlib${ii}$_dla_gerpvgrw
#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
pure real(${rk}$) module function stdlib${ii}$_${ri}$la_gerpvgrw( n, ncols, a, lda, af,ldaf )
!! DLA_GERPVGRW: computes the reciprocal pivot growth factor
!! norm(A)/norm(U). The "max absolute element" norm is used. If this is
!! much less than 1, the stability of the LU factorization of the
!! (equilibrated) matrix A could be poor. This also means that the
!! solution X, estimated condition numbers, and error bounds could be
!! unreliable.
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_${rk}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
integer(${ik}$), intent(in) :: n, ncols, lda, ldaf
! Array Arguments
real(${rk}$), intent(in) :: a(lda,*), af(ldaf,*)
! =====================================================================
! Local Scalars
integer(${ik}$) :: i, j
real(${rk}$) :: amax, umax, rpvgrw
! Intrinsic Functions
! Executable Statements
rpvgrw = one
do j = 1, ncols
amax = zero
umax = zero
do i = 1, n
amax = max( abs( a( i, j ) ), amax )
end do
do i = 1, j
umax = max( abs( af( i, j ) ), umax )
end do
if ( umax /= zero ) then
rpvgrw = min( amax / umax, rpvgrw )
end if
end do
stdlib${ii}$_${ri}$la_gerpvgrw = rpvgrw
end function stdlib${ii}$_${ri}$la_gerpvgrw
#:endif
#:endfor
pure real(sp) module function stdlib${ii}$_cla_gerpvgrw( n, ncols, a, lda, af, ldaf )
!! CLA_GERPVGRW computes the reciprocal pivot growth factor
!! norm(A)/norm(U). The "max absolute element" norm is used. If this is
!! much less than 1, the stability of the LU factorization of the
!! (equilibrated) matrix A could be poor. This also means that the
!! solution X, estimated condition numbers, and error bounds could be
!! unreliable.
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
integer(${ik}$), intent(in) :: n, ncols, lda, ldaf
! Array Arguments
complex(sp), intent(in) :: a(lda,*), af(ldaf,*)
! =====================================================================
! Local Scalars
integer(${ik}$) :: i, j
real(sp) :: amax, umax, rpvgrw
complex(sp) :: zdum
! Intrinsic Functions
! Statement Functions
real(sp) :: cabs1
! Statement Function Definitions
cabs1( zdum ) = abs( real( zdum,KIND=sp) ) + abs( aimag( zdum ) )
! Executable Statements
rpvgrw = one
do j = 1, ncols
amax = zero
umax = zero
do i = 1, n
amax = max( cabs1( a( i, j ) ), amax )
end do
do i = 1, j
umax = max( cabs1( af( i, j ) ), umax )
end do
if ( umax /= 0.0_sp ) then
rpvgrw = min( amax / umax, rpvgrw )
end if
end do
stdlib${ii}$_cla_gerpvgrw = rpvgrw
end function stdlib${ii}$_cla_gerpvgrw
pure real(dp) module function stdlib${ii}$_zla_gerpvgrw( n, ncols, a, lda, af,ldaf )
!! ZLA_GERPVGRW computes the reciprocal pivot growth factor
!! norm(A)/norm(U). The "max absolute element" norm is used. If this is
!! much less than 1, the stability of the LU factorization of the
!! (equilibrated) matrix A could be poor. This also means that the
!! solution X, estimated condition numbers, and error bounds could be
!! unreliable.
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
integer(${ik}$), intent(in) :: n, ncols, lda, ldaf
! Array Arguments
complex(dp), intent(in) :: a(lda,*), af(ldaf,*)
! =====================================================================
! Local Scalars
integer(${ik}$) :: i, j
real(dp) :: amax, umax, rpvgrw
complex(dp) :: zdum
! Intrinsic Functions
! Statement Functions
real(dp) :: cabs1
! Statement Function Definitions
cabs1( zdum ) = abs( real( zdum,KIND=dp) ) + abs( aimag( zdum ) )
! Executable Statements
rpvgrw = one
do j = 1, ncols
amax = zero
umax = zero
do i = 1, n
amax = max( cabs1( a( i, j ) ), amax )
end do
do i = 1, j
umax = max( cabs1( af( i, j ) ), umax )
end do
if ( umax /= zero ) then
rpvgrw = min( amax / umax, rpvgrw )
end if
end do
stdlib${ii}$_zla_gerpvgrw = rpvgrw
end function stdlib${ii}$_zla_gerpvgrw
#:for ck,ct,ci in CMPLX_KINDS_TYPES
#:if not ck in ["sp","dp"]
pure real(${ck}$) module function stdlib${ii}$_${ci}$la_gerpvgrw( n, ncols, a, lda, af,ldaf )
!! ZLA_GERPVGRW: computes the reciprocal pivot growth factor
!! norm(A)/norm(U). The "max absolute element" norm is used. If this is
!! much less than 1, the stability of the LU factorization of the
!! (equilibrated) matrix A could be poor. This also means that the
!! solution X, estimated condition numbers, and error bounds could be
!! unreliable.
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_${ck}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
integer(${ik}$), intent(in) :: n, ncols, lda, ldaf
! Array Arguments
complex(${ck}$), intent(in) :: a(lda,*), af(ldaf,*)
! =====================================================================
! Local Scalars
integer(${ik}$) :: i, j
real(${ck}$) :: amax, umax, rpvgrw
complex(${ck}$) :: zdum
! Intrinsic Functions
! Statement Functions
real(${ck}$) :: cabs1
! Statement Function Definitions
cabs1( zdum ) = abs( real( zdum,KIND=${ck}$) ) + abs( aimag( zdum ) )
! Executable Statements
rpvgrw = one
do j = 1, ncols
amax = zero
umax = zero
do i = 1, n
amax = max( cabs1( a( i, j ) ), amax )
end do
do i = 1, j
umax = max( cabs1( af( i, j ) ), umax )
end do
if ( umax /= zero ) then
rpvgrw = min( amax / umax, rpvgrw )
end if
end do
stdlib${ii}$_${ci}$la_gerpvgrw = rpvgrw
end function stdlib${ii}$_${ci}$la_gerpvgrw
#:endif
#:endfor
real(sp) module function stdlib${ii}$_cla_gbrcond_c( trans, n, kl, ku, ab, ldab, afb,ldafb, ipiv, c, &
!! CLA_GBRCOND_C Computes the infinity norm condition number of
!! op(A) * inv(diag(C)) where C is a REAL vector.
capply, info, work,rwork )
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: trans
logical(lk), intent(in) :: capply
integer(${ik}$), intent(in) :: n, kl, ku, ldab, ldafb
integer(${ik}$) :: kd, ke
integer(${ik}$), intent(out) :: info
! Array Arguments
integer(${ik}$), intent(in) :: ipiv(*)
complex(sp), intent(in) :: ab(ldab,*), afb(ldafb,*)
complex(sp), intent(out) :: work(*)
real(sp), intent(in) :: c(*)
real(sp), intent(out) :: rwork(*)
! =====================================================================
! Local Scalars
logical(lk) :: notrans
integer(${ik}$) :: kase, i, j
real(sp) :: ainvnm, anorm, tmp
complex(sp) :: zdum
! Local Arrays
integer(${ik}$) :: isave(3_${ik}$)
! Intrinsic Functions
! Statement Functions
real(sp) :: cabs1
! Statement Function Definitions
cabs1( zdum ) = abs( real( zdum,KIND=sp) ) + abs( aimag( zdum ) )
! Executable Statements
stdlib${ii}$_cla_gbrcond_c = zero
info = 0_${ik}$
notrans = stdlib_lsame( trans, 'N' )
if ( .not. notrans .and. .not. stdlib_lsame( trans, 'T' ) .and. .not.stdlib_lsame( &
trans, 'C' ) ) then
info = -1_${ik}$
else if( n<0_${ik}$ ) then
info = -2_${ik}$
else if( kl<0_${ik}$ .or. kl>n-1 ) then
info = -3_${ik}$
else if( ku<0_${ik}$ .or. ku>n-1 ) then
info = -4_${ik}$
else if( ldab<kl+ku+1 ) then
info = -6_${ik}$
else if( ldafb<2_${ik}$*kl+ku+1 ) then
info = -8_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'CLA_GBRCOND_C', -info )
return
end if
! compute norm of op(a)*op2(c).
anorm = zero
kd = ku + 1_${ik}$
ke = kl + 1_${ik}$
if ( notrans ) then
do i = 1, n
tmp = zero
if ( capply ) then
do j = max( i-kl, 1 ), min( i+ku, n )
tmp = tmp + cabs1( ab( kd+i-j, j ) ) / c( j )
end do
else
do j = max( i-kl, 1 ), min( i+ku, n )
tmp = tmp + cabs1( ab( kd+i-j, j ) )
end do
end if
rwork( i ) = tmp
anorm = max( anorm, tmp )
end do
else
do i = 1, n
tmp = zero
if ( capply ) then
do j = max( i-kl, 1 ), min( i+ku, n )
tmp = tmp + cabs1( ab( ke-i+j, i ) ) / c( j )
end do
else
do j = max( i-kl, 1 ), min( i+ku, n )
tmp = tmp + cabs1( ab( ke-i+j, i ) )
end do
end if
rwork( i ) = tmp
anorm = max( anorm, tmp )
end do
end if
! quick return if possible.
if( n==0_${ik}$ ) then
stdlib${ii}$_cla_gbrcond_c = one
return
else if( anorm == zero ) then
return
end if
! estimate the norm of inv(op(a)).
ainvnm = zero
kase = 0_${ik}$
10 continue
call stdlib${ii}$_clacn2( n, work( n+1 ), work, ainvnm, kase, isave )
if( kase/=0_${ik}$ ) then
if( kase==2_${ik}$ ) then
! multiply by r.
do i = 1, n
work( i ) = work( i ) * rwork( i )
end do
if ( notrans ) then
call stdlib${ii}$_cgbtrs( 'NO TRANSPOSE', n, kl, ku, 1_${ik}$, afb, ldafb,ipiv, work, n, &
info )
else
call stdlib${ii}$_cgbtrs( 'CONJUGATE TRANSPOSE', n, kl, ku, 1_${ik}$, afb,ldafb, ipiv, &
work, n, info )
endif
! multiply by inv(c).
if ( capply ) then
do i = 1, n
work( i ) = work( i ) * c( i )
end do
end if
else
! multiply by inv(c**h).
if ( capply ) then
do i = 1, n
work( i ) = work( i ) * c( i )
end do
end if
if ( notrans ) then
call stdlib${ii}$_cgbtrs( 'CONJUGATE TRANSPOSE', n, kl, ku, 1_${ik}$, afb,ldafb, ipiv, &
work, n, info )
else
call stdlib${ii}$_cgbtrs( 'NO TRANSPOSE', n, kl, ku, 1_${ik}$, afb, ldafb,ipiv, work, n, &
info )
end if
! multiply by r.
do i = 1, n
work( i ) = work( i ) * rwork( i )
end do
end if
go to 10
end if
! compute the estimate of the reciprocal condition number.
if( ainvnm /= zero )stdlib${ii}$_cla_gbrcond_c = one / ainvnm
return
end function stdlib${ii}$_cla_gbrcond_c
real(dp) module function stdlib${ii}$_zla_gbrcond_c( trans, n, kl, ku, ab,ldab, afb, ldafb, ipiv,c, &
!! ZLA_GBRCOND_C Computes the infinity norm condition number of
!! op(A) * inv(diag(C)) where C is a DOUBLE PRECISION vector.
capply, info, work,rwork )
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: trans
logical(lk), intent(in) :: capply
integer(${ik}$), intent(in) :: n, kl, ku, ldab, ldafb
integer(${ik}$) :: kd, ke
integer(${ik}$), intent(out) :: info
! Array Arguments
integer(${ik}$), intent(in) :: ipiv(*)
complex(dp), intent(in) :: ab(ldab,*), afb(ldafb,*)
complex(dp), intent(out) :: work(*)
real(dp), intent(in) :: c(*)
real(dp), intent(out) :: rwork(*)
! =====================================================================
! Local Scalars
logical(lk) :: notrans
integer(${ik}$) :: kase, i, j
real(dp) :: ainvnm, anorm, tmp
complex(dp) :: zdum
! Local Arrays
integer(${ik}$) :: isave(3_${ik}$)
! Intrinsic Functions
! Statement Functions
real(dp) :: cabs1
! Statement Function Definitions
cabs1( zdum ) = abs( real( zdum,KIND=dp) ) + abs( aimag( zdum ) )
! Executable Statements
stdlib${ii}$_zla_gbrcond_c = zero
info = 0_${ik}$
notrans = stdlib_lsame( trans, 'N' )
if ( .not. notrans .and. .not. stdlib_lsame( trans, 'T' ) .and. .not.stdlib_lsame( &
trans, 'C' ) ) then
info = -1_${ik}$
else if( n<0_${ik}$ ) then
info = -2_${ik}$
else if( kl<0_${ik}$ .or. kl>n-1 ) then
info = -3_${ik}$
else if( ku<0_${ik}$ .or. ku>n-1 ) then
info = -4_${ik}$
else if( ldab<kl+ku+1 ) then
info = -6_${ik}$
else if( ldafb<2_${ik}$*kl+ku+1 ) then
info = -8_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'ZLA_GBRCOND_C', -info )
return
end if
! compute norm of op(a)*op2(c).
anorm = zero
kd = ku + 1_${ik}$
ke = kl + 1_${ik}$
if ( notrans ) then
do i = 1, n
tmp = zero
if ( capply ) then
do j = max( i-kl, 1 ), min( i+ku, n )
tmp = tmp + cabs1( ab( kd+i-j, j ) ) / c( j )
end do
else
do j = max( i-kl, 1 ), min( i+ku, n )
tmp = tmp + cabs1( ab( kd+i-j, j ) )
end do
end if
rwork( i ) = tmp
anorm = max( anorm, tmp )
end do
else
do i = 1, n
tmp = zero
if ( capply ) then
do j = max( i-kl, 1 ), min( i+ku, n )
tmp = tmp + cabs1( ab( ke-i+j, i ) ) / c( j )
end do
else
do j = max( i-kl, 1 ), min( i+ku, n )
tmp = tmp + cabs1( ab( ke-i+j, i ) )
end do
end if
rwork( i ) = tmp
anorm = max( anorm, tmp )
end do
end if
! quick return if possible.
if( n==0_${ik}$ ) then
stdlib${ii}$_zla_gbrcond_c = one
return
else if( anorm == zero ) then
return
end if
! estimate the norm of inv(op(a)).
ainvnm = zero
kase = 0_${ik}$
10 continue
call stdlib${ii}$_zlacn2( n, work( n+1 ), work, ainvnm, kase, isave )
if( kase/=0_${ik}$ ) then
if( kase==2_${ik}$ ) then
! multiply by r.
do i = 1, n
work( i ) = work( i ) * rwork( i )
end do
if ( notrans ) then
call stdlib${ii}$_zgbtrs( 'NO TRANSPOSE', n, kl, ku, 1_${ik}$, afb, ldafb,ipiv, work, n, &
info )
else
call stdlib${ii}$_zgbtrs( 'CONJUGATE TRANSPOSE', n, kl, ku, 1_${ik}$, afb,ldafb, ipiv, &
work, n, info )
endif
! multiply by inv(c).
if ( capply ) then
do i = 1, n
work( i ) = work( i ) * c( i )
end do
end if
else
! multiply by inv(c**h).
if ( capply ) then
do i = 1, n
work( i ) = work( i ) * c( i )
end do
end if
if ( notrans ) then
call stdlib${ii}$_zgbtrs( 'CONJUGATE TRANSPOSE', n, kl, ku, 1_${ik}$, afb,ldafb, ipiv, &
work, n, info )
else
call stdlib${ii}$_zgbtrs( 'NO TRANSPOSE', n, kl, ku, 1_${ik}$, afb, ldafb,ipiv, work, n, &
info )
end if
! multiply by r.
do i = 1, n
work( i ) = work( i ) * rwork( i )
end do
end if
go to 10
end if
! compute the estimate of the reciprocal condition number.
if( ainvnm /= zero )stdlib${ii}$_zla_gbrcond_c = one / ainvnm
return
end function stdlib${ii}$_zla_gbrcond_c
#:for ck,ct,ci in CMPLX_KINDS_TYPES
#:if not ck in ["sp","dp"]
real(${ck}$) module function stdlib${ii}$_${ci}$la_gbrcond_c( trans, n, kl, ku, ab,ldab, afb, ldafb, ipiv,c, &
!! ZLA_GBRCOND_C: Computes the infinity norm condition number of
!! op(A) * inv(diag(C)) where C is a DOUBLE PRECISION vector.
capply, info, work,rwork )
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_${ck}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: trans
logical(lk), intent(in) :: capply
integer(${ik}$), intent(in) :: n, kl, ku, ldab, ldafb
integer(${ik}$) :: kd, ke
integer(${ik}$), intent(out) :: info
! Array Arguments
integer(${ik}$), intent(in) :: ipiv(*)
complex(${ck}$), intent(in) :: ab(ldab,*), afb(ldafb,*)
complex(${ck}$), intent(out) :: work(*)
real(${ck}$), intent(in) :: c(*)
real(${ck}$), intent(out) :: rwork(*)
! =====================================================================
! Local Scalars
logical(lk) :: notrans
integer(${ik}$) :: kase, i, j
real(${ck}$) :: ainvnm, anorm, tmp
complex(${ck}$) :: zdum
! Local Arrays
integer(${ik}$) :: isave(3_${ik}$)
! Intrinsic Functions
! Statement Functions
real(${ck}$) :: cabs1
! Statement Function Definitions
cabs1( zdum ) = abs( real( zdum,KIND=${ck}$) ) + abs( aimag( zdum ) )
! Executable Statements
stdlib${ii}$_${ci}$la_gbrcond_c = zero
info = 0_${ik}$
notrans = stdlib_lsame( trans, 'N' )
if ( .not. notrans .and. .not. stdlib_lsame( trans, 'T' ) .and. .not.stdlib_lsame( &
trans, 'C' ) ) then
info = -1_${ik}$
else if( n<0_${ik}$ ) then
info = -2_${ik}$
else if( kl<0_${ik}$ .or. kl>n-1 ) then
info = -3_${ik}$
else if( ku<0_${ik}$ .or. ku>n-1 ) then
info = -4_${ik}$
else if( ldab<kl+ku+1 ) then
info = -6_${ik}$
else if( ldafb<2_${ik}$*kl+ku+1 ) then
info = -8_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'ZLA_GBRCOND_C', -info )
return
end if
! compute norm of op(a)*op2(c).
anorm = zero
kd = ku + 1_${ik}$
ke = kl + 1_${ik}$
if ( notrans ) then
do i = 1, n
tmp = zero
if ( capply ) then
do j = max( i-kl, 1 ), min( i+ku, n )
tmp = tmp + cabs1( ab( kd+i-j, j ) ) / c( j )
end do
else
do j = max( i-kl, 1 ), min( i+ku, n )
tmp = tmp + cabs1( ab( kd+i-j, j ) )
end do
end if
rwork( i ) = tmp
anorm = max( anorm, tmp )
end do
else
do i = 1, n
tmp = zero
if ( capply ) then
do j = max( i-kl, 1 ), min( i+ku, n )
tmp = tmp + cabs1( ab( ke-i+j, i ) ) / c( j )
end do
else
do j = max( i-kl, 1 ), min( i+ku, n )
tmp = tmp + cabs1( ab( ke-i+j, i ) )
end do
end if
rwork( i ) = tmp
anorm = max( anorm, tmp )
end do
end if
! quick return if possible.
if( n==0_${ik}$ ) then
stdlib${ii}$_${ci}$la_gbrcond_c = one
return
else if( anorm == zero ) then
return
end if
! estimate the norm of inv(op(a)).
ainvnm = zero
kase = 0_${ik}$
10 continue
call stdlib${ii}$_${ci}$lacn2( n, work( n+1 ), work, ainvnm, kase, isave )
if( kase/=0_${ik}$ ) then
if( kase==2_${ik}$ ) then
! multiply by r.
do i = 1, n
work( i ) = work( i ) * rwork( i )
end do
if ( notrans ) then
call stdlib${ii}$_${ci}$gbtrs( 'NO TRANSPOSE', n, kl, ku, 1_${ik}$, afb, ldafb,ipiv, work, n, &
info )
else
call stdlib${ii}$_${ci}$gbtrs( 'CONJUGATE TRANSPOSE', n, kl, ku, 1_${ik}$, afb,ldafb, ipiv, &
work, n, info )
endif
! multiply by inv(c).
if ( capply ) then
do i = 1, n
work( i ) = work( i ) * c( i )
end do
end if
else
! multiply by inv(c**h).
if ( capply ) then
do i = 1, n
work( i ) = work( i ) * c( i )
end do
end if
if ( notrans ) then
call stdlib${ii}$_${ci}$gbtrs( 'CONJUGATE TRANSPOSE', n, kl, ku, 1_${ik}$, afb,ldafb, ipiv, &
work, n, info )
else
call stdlib${ii}$_${ci}$gbtrs( 'NO TRANSPOSE', n, kl, ku, 1_${ik}$, afb, ldafb,ipiv, work, n, &
info )
end if
! multiply by r.
do i = 1, n
work( i ) = work( i ) * rwork( i )
end do
end if
go to 10
end if
! compute the estimate of the reciprocal condition number.
if( ainvnm /= zero )stdlib${ii}$_${ci}$la_gbrcond_c = one / ainvnm
return
end function stdlib${ii}$_${ci}$la_gbrcond_c
#:endif
#:endfor
real(sp) module function stdlib${ii}$_cla_gercond_c( trans, n, a, lda, af, ldaf, ipiv, c,capply, info, &
!! CLA_GERCOND_C computes the infinity norm condition number of
!! op(A) * inv(diag(C)) where C is a REAL vector.
work, rwork )
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: trans
logical(lk), intent(in) :: capply
integer(${ik}$), intent(in) :: n, lda, ldaf
integer(${ik}$), intent(out) :: info
! Array Arguments
integer(${ik}$), intent(in) :: ipiv(*)
complex(sp), intent(in) :: a(lda,*), af(ldaf,*)
complex(sp), intent(out) :: work(*)
real(sp), intent(in) :: c(*)
real(sp), intent(out) :: rwork(*)
! =====================================================================
! Local Scalars
logical(lk) :: notrans
integer(${ik}$) :: kase, i, j
real(sp) :: ainvnm, anorm, tmp
complex(sp) :: zdum
! Local Arrays
integer(${ik}$) :: isave(3_${ik}$)
! Intrinsic Functions
! Statement Functions
real(sp) :: cabs1
! Statement Function Definitions
cabs1( zdum ) = abs( real( zdum,KIND=sp) ) + abs( aimag( zdum ) )
! Executable Statements
stdlib${ii}$_cla_gercond_c = zero
info = 0_${ik}$
notrans = stdlib_lsame( trans, 'N' )
if ( .not. notrans .and. .not. stdlib_lsame( trans, 'T' ) .and. .not.stdlib_lsame( &
trans, 'C' ) ) then
info = -1_${ik}$
else if( n<0_${ik}$ ) then
info = -2_${ik}$
else if( lda<max( 1_${ik}$, n ) ) then
info = -4_${ik}$
else if( ldaf<max( 1_${ik}$, n ) ) then
info = -6_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'CLA_GERCOND_C', -info )
return
end if
! compute norm of op(a)*op2(c).
anorm = zero
if ( notrans ) then
do i = 1, n
tmp = zero
if ( capply ) then
do j = 1, n
tmp = tmp + cabs1( a( i, j ) ) / c( j )
end do
else
do j = 1, n
tmp = tmp + cabs1( a( i, j ) )
end do
end if
rwork( i ) = tmp
anorm = max( anorm, tmp )
end do
else
do i = 1, n
tmp = zero
if ( capply ) then
do j = 1, n
tmp = tmp + cabs1( a( j, i ) ) / c( j )
end do
else
do j = 1, n
tmp = tmp + cabs1( a( j, i ) )
end do
end if
rwork( i ) = tmp
anorm = max( anorm, tmp )
end do
end if
! quick return if possible.
if( n==0_${ik}$ ) then
stdlib${ii}$_cla_gercond_c = one
return
else if( anorm == zero ) then
return
end if
! estimate the norm of inv(op(a)).
ainvnm = zero
kase = 0_${ik}$
10 continue
call stdlib${ii}$_clacn2( n, work( n+1 ), work, ainvnm, kase, isave )
if( kase/=0_${ik}$ ) then
if( kase==2_${ik}$ ) then
! multiply by r.
do i = 1, n
work( i ) = work( i ) * rwork( i )
end do
if (notrans) then
call stdlib${ii}$_cgetrs( 'NO TRANSPOSE', n, 1_${ik}$, af, ldaf, ipiv,work, n, info )
else
call stdlib${ii}$_cgetrs( 'CONJUGATE TRANSPOSE', n, 1_${ik}$, af, ldaf, ipiv,work, n, info &
)
endif
! multiply by inv(c).
if ( capply ) then
do i = 1, n
work( i ) = work( i ) * c( i )
end do
end if
else
! multiply by inv(c**h).
if ( capply ) then
do i = 1, n
work( i ) = work( i ) * c( i )
end do
end if
if ( notrans ) then
call stdlib${ii}$_cgetrs( 'CONJUGATE TRANSPOSE', n, 1_${ik}$, af, ldaf, ipiv,work, n, info &
)
else
call stdlib${ii}$_cgetrs( 'NO TRANSPOSE', n, 1_${ik}$, af, ldaf, ipiv,work, n, info )
end if
! multiply by r.
do i = 1, n
work( i ) = work( i ) * rwork( i )
end do
end if
go to 10
end if
! compute the estimate of the reciprocal condition number.
if( ainvnm /= zero )stdlib${ii}$_cla_gercond_c = one / ainvnm
return
end function stdlib${ii}$_cla_gercond_c
real(dp) module function stdlib${ii}$_zla_gercond_c( trans, n, a, lda, af,ldaf, ipiv, c, capply,info, &
!! ZLA_GERCOND_C computes the infinity norm condition number of
!! op(A) * inv(diag(C)) where C is a DOUBLE PRECISION vector.
work, rwork )
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: trans
logical(lk), intent(in) :: capply
integer(${ik}$), intent(in) :: n, lda, ldaf
integer(${ik}$), intent(out) :: info
! Array Arguments
integer(${ik}$), intent(in) :: ipiv(*)
complex(dp), intent(in) :: a(lda,*), af(ldaf,*)
complex(dp), intent(out) :: work(*)
real(dp), intent(in) :: c(*)
real(dp), intent(out) :: rwork(*)
! =====================================================================
! Local Scalars
logical(lk) :: notrans
integer(${ik}$) :: kase, i, j
real(dp) :: ainvnm, anorm, tmp
complex(dp) :: zdum
! Local Arrays
integer(${ik}$) :: isave(3_${ik}$)
! Intrinsic Functions
! Statement Functions
real(dp) :: cabs1
! Statement Function Definitions
cabs1( zdum ) = abs( real( zdum,KIND=dp) ) + abs( aimag( zdum ) )
! Executable Statements
stdlib${ii}$_zla_gercond_c = zero
info = 0_${ik}$
notrans = stdlib_lsame( trans, 'N' )
if ( .not. notrans .and. .not. stdlib_lsame( trans, 'T' ) .and. .not.stdlib_lsame( &
trans, 'C' ) ) then
info = -1_${ik}$
else if( n<0_${ik}$ ) then
info = -2_${ik}$
else if( lda<max( 1_${ik}$, n ) ) then
info = -4_${ik}$
else if( ldaf<max( 1_${ik}$, n ) ) then
info = -6_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'ZLA_GERCOND_C', -info )
return
end if
! compute norm of op(a)*op2(c).
anorm = zero
if ( notrans ) then
do i = 1, n
tmp = zero
if ( capply ) then
do j = 1, n
tmp = tmp + cabs1( a( i, j ) ) / c( j )
end do
else
do j = 1, n
tmp = tmp + cabs1( a( i, j ) )
end do
end if
rwork( i ) = tmp
anorm = max( anorm, tmp )
end do
else
do i = 1, n
tmp = zero
if ( capply ) then
do j = 1, n
tmp = tmp + cabs1( a( j, i ) ) / c( j )
end do
else
do j = 1, n
tmp = tmp + cabs1( a( j, i ) )
end do
end if
rwork( i ) = tmp
anorm = max( anorm, tmp )
end do
end if
! quick return if possible.
if( n==0_${ik}$ ) then
stdlib${ii}$_zla_gercond_c = one
return
else if( anorm == zero ) then
return
end if
! estimate the norm of inv(op(a)).
ainvnm = zero
kase = 0_${ik}$
10 continue
call stdlib${ii}$_zlacn2( n, work( n+1 ), work, ainvnm, kase, isave )
if( kase/=0_${ik}$ ) then
if( kase==2_${ik}$ ) then
! multiply by r.
do i = 1, n
work( i ) = work( i ) * rwork( i )
end do
if (notrans) then
call stdlib${ii}$_zgetrs( 'NO TRANSPOSE', n, 1_${ik}$, af, ldaf, ipiv,work, n, info )
else
call stdlib${ii}$_zgetrs( 'CONJUGATE TRANSPOSE', n, 1_${ik}$, af, ldaf, ipiv,work, n, info &
)
endif
! multiply by inv(c).
if ( capply ) then
do i = 1, n
work( i ) = work( i ) * c( i )
end do
end if
else
! multiply by inv(c**h).
if ( capply ) then
do i = 1, n
work( i ) = work( i ) * c( i )
end do
end if
if ( notrans ) then
call stdlib${ii}$_zgetrs( 'CONJUGATE TRANSPOSE', n, 1_${ik}$, af, ldaf, ipiv,work, n, info &
)
else
call stdlib${ii}$_zgetrs( 'NO TRANSPOSE', n, 1_${ik}$, af, ldaf, ipiv,work, n, info )
end if
! multiply by r.
do i = 1, n
work( i ) = work( i ) * rwork( i )
end do
end if
go to 10
end if
! compute the estimate of the reciprocal condition number.
if( ainvnm /= zero )stdlib${ii}$_zla_gercond_c = one / ainvnm
return
end function stdlib${ii}$_zla_gercond_c
#:for ck,ct,ci in CMPLX_KINDS_TYPES
#:if not ck in ["sp","dp"]
real(${ck}$) module function stdlib${ii}$_${ci}$la_gercond_c( trans, n, a, lda, af,ldaf, ipiv, c, capply,info, &
!! ZLA_GERCOND_C: computes the infinity norm condition number of
!! op(A) * inv(diag(C)) where C is a DOUBLE PRECISION vector.
work, rwork )
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_${ck}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: trans
logical(lk), intent(in) :: capply
integer(${ik}$), intent(in) :: n, lda, ldaf
integer(${ik}$), intent(out) :: info
! Array Arguments
integer(${ik}$), intent(in) :: ipiv(*)
complex(${ck}$), intent(in) :: a(lda,*), af(ldaf,*)
complex(${ck}$), intent(out) :: work(*)
real(${ck}$), intent(in) :: c(*)
real(${ck}$), intent(out) :: rwork(*)
! =====================================================================
! Local Scalars
logical(lk) :: notrans
integer(${ik}$) :: kase, i, j
real(${ck}$) :: ainvnm, anorm, tmp
complex(${ck}$) :: zdum
! Local Arrays
integer(${ik}$) :: isave(3_${ik}$)
! Intrinsic Functions
! Statement Functions
real(${ck}$) :: cabs1
! Statement Function Definitions
cabs1( zdum ) = abs( real( zdum,KIND=${ck}$) ) + abs( aimag( zdum ) )
! Executable Statements
stdlib${ii}$_${ci}$la_gercond_c = zero
info = 0_${ik}$
notrans = stdlib_lsame( trans, 'N' )
if ( .not. notrans .and. .not. stdlib_lsame( trans, 'T' ) .and. .not.stdlib_lsame( &
trans, 'C' ) ) then
info = -1_${ik}$
else if( n<0_${ik}$ ) then
info = -2_${ik}$
else if( lda<max( 1_${ik}$, n ) ) then
info = -4_${ik}$
else if( ldaf<max( 1_${ik}$, n ) ) then
info = -6_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'ZLA_GERCOND_C', -info )
return
end if
! compute norm of op(a)*op2(c).
anorm = zero
if ( notrans ) then
do i = 1, n
tmp = zero
if ( capply ) then
do j = 1, n
tmp = tmp + cabs1( a( i, j ) ) / c( j )
end do
else
do j = 1, n
tmp = tmp + cabs1( a( i, j ) )
end do
end if
rwork( i ) = tmp
anorm = max( anorm, tmp )
end do
else
do i = 1, n
tmp = zero
if ( capply ) then
do j = 1, n
tmp = tmp + cabs1( a( j, i ) ) / c( j )
end do
else
do j = 1, n
tmp = tmp + cabs1( a( j, i ) )
end do
end if
rwork( i ) = tmp
anorm = max( anorm, tmp )
end do
end if
! quick return if possible.
if( n==0_${ik}$ ) then
stdlib${ii}$_${ci}$la_gercond_c = one
return
else if( anorm == zero ) then
return
end if
! estimate the norm of inv(op(a)).
ainvnm = zero
kase = 0_${ik}$
10 continue
call stdlib${ii}$_${ci}$lacn2( n, work( n+1 ), work, ainvnm, kase, isave )
if( kase/=0_${ik}$ ) then
if( kase==2_${ik}$ ) then
! multiply by r.
do i = 1, n
work( i ) = work( i ) * rwork( i )
end do
if (notrans) then
call stdlib${ii}$_${ci}$getrs( 'NO TRANSPOSE', n, 1_${ik}$, af, ldaf, ipiv,work, n, info )
else
call stdlib${ii}$_${ci}$getrs( 'CONJUGATE TRANSPOSE', n, 1_${ik}$, af, ldaf, ipiv,work, n, info &
)
endif
! multiply by inv(c).
if ( capply ) then
do i = 1, n
work( i ) = work( i ) * c( i )
end do
end if
else
! multiply by inv(c**h).
if ( capply ) then
do i = 1, n
work( i ) = work( i ) * c( i )
end do
end if
if ( notrans ) then
call stdlib${ii}$_${ci}$getrs( 'CONJUGATE TRANSPOSE', n, 1_${ik}$, af, ldaf, ipiv,work, n, info &
)
else
call stdlib${ii}$_${ci}$getrs( 'NO TRANSPOSE', n, 1_${ik}$, af, ldaf, ipiv,work, n, info )
end if
! multiply by r.
do i = 1, n
work( i ) = work( i ) * rwork( i )
end do
end if
go to 10
end if
! compute the estimate of the reciprocal condition number.
if( ainvnm /= zero )stdlib${ii}$_${ci}$la_gercond_c = one / ainvnm
return
end function stdlib${ii}$_${ci}$la_gercond_c
#:endif
#:endfor
real(sp) module function stdlib${ii}$_cla_hercond_c( uplo, n, a, lda, af, ldaf, ipiv, c,capply, info, &
!! CLA_HERCOND_C computes the infinity norm condition number of
!! op(A) * inv(diag(C)) where C is a REAL vector.
work, rwork )
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: uplo
logical(lk), intent(in) :: capply
integer(${ik}$), intent(in) :: n, lda, ldaf
integer(${ik}$), intent(out) :: info
! Array Arguments
integer(${ik}$), intent(in) :: ipiv(*)
complex(sp), intent(in) :: a(lda,*), af(ldaf,*)
complex(sp), intent(out) :: work(*)
real(sp), intent(in) :: c(*)
real(sp), intent(out) :: rwork(*)
! =====================================================================
! Local Scalars
integer(${ik}$) :: kase, i, j
real(sp) :: ainvnm, anorm, tmp
logical(lk) :: up, upper
complex(sp) :: zdum
! Local Arrays
integer(${ik}$) :: isave(3_${ik}$)
! Intrinsic Functions
! Statement Functions
real(sp) :: cabs1
! Statement Function Definitions
cabs1( zdum ) = abs( real( zdum,KIND=sp) ) + abs( aimag( zdum ) )
! Executable Statements
stdlib${ii}$_cla_hercond_c = zero
info = 0_${ik}$
upper = stdlib_lsame( uplo, 'U' )
if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then
info = -1_${ik}$
else if( n<0_${ik}$ ) then
info = -2_${ik}$
else if( lda<max( 1_${ik}$, n ) ) then
info = -4_${ik}$
else if( ldaf<max( 1_${ik}$, n ) ) then
info = -6_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'CLA_HERCOND_C', -info )
return
end if
up = .false.
if ( stdlib_lsame( uplo, 'U' ) ) up = .true.
! compute norm of op(a)*op2(c).
anorm = zero
if ( up ) then
do i = 1, n
tmp = zero
if ( capply ) then
do j = 1, i
tmp = tmp + cabs1( a( j, i ) ) / c( j )
end do
do j = i+1, n
tmp = tmp + cabs1( a( i, j ) ) / c( j )
end do
else
do j = 1, i
tmp = tmp + cabs1( a( j, i ) )
end do
do j = i+1, n
tmp = tmp + cabs1( a( i, j ) )
end do
end if
rwork( i ) = tmp
anorm = max( anorm, tmp )
end do
else
do i = 1, n
tmp = zero
if ( capply ) then
do j = 1, i
tmp = tmp + cabs1( a( i, j ) ) / c( j )
end do
do j = i+1, n
tmp = tmp + cabs1( a( j, i ) ) / c( j )
end do
else
do j = 1, i
tmp = tmp + cabs1( a( i, j ) )
end do
do j = i+1, n
tmp = tmp + cabs1( a( j, i ) )
end do
end if
rwork( i ) = tmp
anorm = max( anorm, tmp )
end do
end if
! quick return if possible.
if( n==0_${ik}$ ) then
stdlib${ii}$_cla_hercond_c = one
return
else if( anorm == zero ) then
return
end if
! estimate the norm of inv(op(a)).
ainvnm = zero
kase = 0_${ik}$
10 continue
call stdlib${ii}$_clacn2( n, work( n+1 ), work, ainvnm, kase, isave )
if( kase/=0_${ik}$ ) then
if( kase==2_${ik}$ ) then
! multiply by r.
do i = 1, n
work( i ) = work( i ) * rwork( i )
end do
if ( up ) then
call stdlib${ii}$_chetrs( 'U', n, 1_${ik}$, af, ldaf, ipiv,work, n, info )
else
call stdlib${ii}$_chetrs( 'L', n, 1_${ik}$, af, ldaf, ipiv,work, n, info )
endif
! multiply by inv(c).
if ( capply ) then
do i = 1, n
work( i ) = work( i ) * c( i )
end do
end if
else
! multiply by inv(c**h).
if ( capply ) then
do i = 1, n
work( i ) = work( i ) * c( i )
end do
end if
if ( up ) then
call stdlib${ii}$_chetrs( 'U', n, 1_${ik}$, af, ldaf, ipiv,work, n, info )
else
call stdlib${ii}$_chetrs( 'L', n, 1_${ik}$, af, ldaf, ipiv,work, n, info )
end if
! multiply by r.
do i = 1, n
work( i ) = work( i ) * rwork( i )
end do
end if
go to 10
end if
! compute the estimate of the reciprocal condition number.
if( ainvnm /= zero )stdlib${ii}$_cla_hercond_c = one / ainvnm
return
end function stdlib${ii}$_cla_hercond_c
real(dp) module function stdlib${ii}$_zla_hercond_c( uplo, n, a, lda, af,ldaf, ipiv, c, capply,info, work,&
!! ZLA_HERCOND_C computes the infinity norm condition number of
!! op(A) * inv(diag(C)) where C is a DOUBLE PRECISION vector.
rwork )
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: uplo
logical(lk), intent(in) :: capply
integer(${ik}$), intent(in) :: n, lda, ldaf
integer(${ik}$), intent(out) :: info
! Array Arguments
integer(${ik}$), intent(in) :: ipiv(*)
complex(dp), intent(in) :: a(lda,*), af(ldaf,*)
complex(dp), intent(out) :: work(*)
real(dp), intent(in) :: c(*)
real(dp), intent(out) :: rwork(*)
! =====================================================================
! Local Scalars
integer(${ik}$) :: kase, i, j
real(dp) :: ainvnm, anorm, tmp
logical(lk) :: up, upper
complex(dp) :: zdum
! Local Arrays
integer(${ik}$) :: isave(3_${ik}$)
! Intrinsic Functions
! Statement Functions
real(dp) :: cabs1
! Statement Function Definitions
cabs1( zdum ) = abs( real( zdum,KIND=dp) ) + abs( aimag( zdum ) )
! Executable Statements
stdlib${ii}$_zla_hercond_c = zero
info = 0_${ik}$
upper = stdlib_lsame( uplo, 'U' )
if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then
info = -1_${ik}$
else if( n<0_${ik}$ ) then
info = -2_${ik}$
else if( lda<max( 1_${ik}$, n ) ) then
info = -4_${ik}$
else if( ldaf<max( 1_${ik}$, n ) ) then
info = -6_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'ZLA_HERCOND_C', -info )
return
end if
up = .false.
if ( stdlib_lsame( uplo, 'U' ) ) up = .true.
! compute norm of op(a)*op2(c).
anorm = zero
if ( up ) then
do i = 1, n
tmp = zero
if ( capply ) then
do j = 1, i
tmp = tmp + cabs1( a( j, i ) ) / c( j )
end do
do j = i+1, n
tmp = tmp + cabs1( a( i, j ) ) / c( j )
end do
else
do j = 1, i
tmp = tmp + cabs1( a( j, i ) )
end do
do j = i+1, n
tmp = tmp + cabs1( a( i, j ) )
end do
end if
rwork( i ) = tmp
anorm = max( anorm, tmp )
end do
else
do i = 1, n
tmp = zero
if ( capply ) then
do j = 1, i
tmp = tmp + cabs1( a( i, j ) ) / c( j )
end do
do j = i+1, n
tmp = tmp + cabs1( a( j, i ) ) / c( j )
end do
else
do j = 1, i
tmp = tmp + cabs1( a( i, j ) )
end do
do j = i+1, n
tmp = tmp + cabs1( a( j, i ) )
end do
end if
rwork( i ) = tmp
anorm = max( anorm, tmp )
end do
end if
! quick return if possible.
if( n==0_${ik}$ ) then
stdlib${ii}$_zla_hercond_c = one
return
else if( anorm == zero ) then
return
end if
! estimate the norm of inv(op(a)).
ainvnm = zero
kase = 0_${ik}$
10 continue
call stdlib${ii}$_zlacn2( n, work( n+1 ), work, ainvnm, kase, isave )
if( kase/=0_${ik}$ ) then
if( kase==2_${ik}$ ) then
! multiply by r.
do i = 1, n
work( i ) = work( i ) * rwork( i )
end do
if ( up ) then
call stdlib${ii}$_zhetrs( 'U', n, 1_${ik}$, af, ldaf, ipiv,work, n, info )
else
call stdlib${ii}$_zhetrs( 'L', n, 1_${ik}$, af, ldaf, ipiv,work, n, info )
endif
! multiply by inv(c).
if ( capply ) then
do i = 1, n
work( i ) = work( i ) * c( i )
end do
end if
else
! multiply by inv(c**h).
if ( capply ) then
do i = 1, n
work( i ) = work( i ) * c( i )
end do
end if
if ( up ) then
call stdlib${ii}$_zhetrs( 'U', n, 1_${ik}$, af, ldaf, ipiv,work, n, info )
else
call stdlib${ii}$_zhetrs( 'L', n, 1_${ik}$, af, ldaf, ipiv,work, n, info )
end if
! multiply by r.
do i = 1, n
work( i ) = work( i ) * rwork( i )
end do
end if
go to 10
end if
! compute the estimate of the reciprocal condition number.
if( ainvnm /= zero )stdlib${ii}$_zla_hercond_c = one / ainvnm
return
end function stdlib${ii}$_zla_hercond_c
#:for ck,ct,ci in CMPLX_KINDS_TYPES
#:if not ck in ["sp","dp"]
real(${ck}$) module function stdlib${ii}$_${ci}$la_hercond_c( uplo, n, a, lda, af,ldaf, ipiv, c, capply,info, work,&
!! ZLA_HERCOND_C: computes the infinity norm condition number of
!! op(A) * inv(diag(C)) where C is a DOUBLE PRECISION vector.
rwork )
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_${ck}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: uplo
logical(lk), intent(in) :: capply
integer(${ik}$), intent(in) :: n, lda, ldaf
integer(${ik}$), intent(out) :: info
! Array Arguments
integer(${ik}$), intent(in) :: ipiv(*)
complex(${ck}$), intent(in) :: a(lda,*), af(ldaf,*)
complex(${ck}$), intent(out) :: work(*)
real(${ck}$), intent(in) :: c(*)
real(${ck}$), intent(out) :: rwork(*)
! =====================================================================
! Local Scalars
integer(${ik}$) :: kase, i, j
real(${ck}$) :: ainvnm, anorm, tmp
logical(lk) :: up, upper
complex(${ck}$) :: zdum
! Local Arrays
integer(${ik}$) :: isave(3_${ik}$)
! Intrinsic Functions
! Statement Functions
real(${ck}$) :: cabs1
! Statement Function Definitions
cabs1( zdum ) = abs( real( zdum,KIND=${ck}$) ) + abs( aimag( zdum ) )
! Executable Statements
stdlib${ii}$_${ci}$la_hercond_c = zero
info = 0_${ik}$
upper = stdlib_lsame( uplo, 'U' )
if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then
info = -1_${ik}$
else if( n<0_${ik}$ ) then
info = -2_${ik}$
else if( lda<max( 1_${ik}$, n ) ) then
info = -4_${ik}$
else if( ldaf<max( 1_${ik}$, n ) ) then
info = -6_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'ZLA_HERCOND_C', -info )
return
end if
up = .false.
if ( stdlib_lsame( uplo, 'U' ) ) up = .true.
! compute norm of op(a)*op2(c).
anorm = zero
if ( up ) then
do i = 1, n
tmp = zero
if ( capply ) then
do j = 1, i
tmp = tmp + cabs1( a( j, i ) ) / c( j )
end do
do j = i+1, n
tmp = tmp + cabs1( a( i, j ) ) / c( j )
end do
else
do j = 1, i
tmp = tmp + cabs1( a( j, i ) )
end do
do j = i+1, n
tmp = tmp + cabs1( a( i, j ) )
end do
end if
rwork( i ) = tmp
anorm = max( anorm, tmp )
end do
else
do i = 1, n
tmp = zero
if ( capply ) then
do j = 1, i
tmp = tmp + cabs1( a( i, j ) ) / c( j )
end do
do j = i+1, n
tmp = tmp + cabs1( a( j, i ) ) / c( j )
end do
else
do j = 1, i
tmp = tmp + cabs1( a( i, j ) )
end do
do j = i+1, n
tmp = tmp + cabs1( a( j, i ) )
end do
end if
rwork( i ) = tmp
anorm = max( anorm, tmp )
end do
end if
! quick return if possible.
if( n==0_${ik}$ ) then
stdlib${ii}$_${ci}$la_hercond_c = one
return
else if( anorm == zero ) then
return
end if
! estimate the norm of inv(op(a)).
ainvnm = zero
kase = 0_${ik}$
10 continue
call stdlib${ii}$_${ci}$lacn2( n, work( n+1 ), work, ainvnm, kase, isave )
if( kase/=0_${ik}$ ) then
if( kase==2_${ik}$ ) then
! multiply by r.
do i = 1, n
work( i ) = work( i ) * rwork( i )
end do
if ( up ) then
call stdlib${ii}$_${ci}$hetrs( 'U', n, 1_${ik}$, af, ldaf, ipiv,work, n, info )
else
call stdlib${ii}$_${ci}$hetrs( 'L', n, 1_${ik}$, af, ldaf, ipiv,work, n, info )
endif
! multiply by inv(c).
if ( capply ) then
do i = 1, n
work( i ) = work( i ) * c( i )
end do
end if
else
! multiply by inv(c**h).
if ( capply ) then
do i = 1, n
work( i ) = work( i ) * c( i )
end do
end if
if ( up ) then
call stdlib${ii}$_${ci}$hetrs( 'U', n, 1_${ik}$, af, ldaf, ipiv,work, n, info )
else
call stdlib${ii}$_${ci}$hetrs( 'L', n, 1_${ik}$, af, ldaf, ipiv,work, n, info )
end if
! multiply by r.
do i = 1, n
work( i ) = work( i ) * rwork( i )
end do
end if
go to 10
end if
! compute the estimate of the reciprocal condition number.
if( ainvnm /= zero )stdlib${ii}$_${ci}$la_hercond_c = one / ainvnm
return
end function stdlib${ii}$_${ci}$la_hercond_c
#:endif
#:endfor
module subroutine stdlib${ii}$_sla_syamv( uplo, n, alpha, a, lda, x, incx, beta, y,incy )
!! SLA_SYAMV performs the matrix-vector operation
!! y := alpha*abs(A)*abs(x) + beta*abs(y),
!! where alpha and beta are scalars, x and y are vectors and A is an
!! n by n symmetric matrix.
!! This function is primarily used in calculating error bounds.
!! To protect against underflow during evaluation, components in
!! the resulting vector are perturbed away from zero by (N+1)
!! times the underflow threshold. To prevent unnecessarily large
!! errors for block-structure embedded in general matrices,
!! "symbolically" zero components are not perturbed. A zero
!! entry is considered "symbolic" if all multiplications involved
!! in computing that entry have at least one zero multiplicand.
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
real(sp), intent(in) :: alpha, beta
integer(${ik}$), intent(in) :: incx, incy, lda, n, uplo
! Array Arguments
real(sp), intent(in) :: a(lda,*), x(*)
real(sp), intent(inout) :: y(*)
! =====================================================================
! Local Scalars
logical(lk) :: symb_zero
real(sp) :: temp, safe1
integer(${ik}$) :: i, info, iy, j, jx, kx, ky
! Intrinsic Functions
! Executable Statements
! test the input parameters.
info = 0_${ik}$
if ( uplo/=stdlib${ii}$_ilauplo( 'U' ) .and.uplo/=stdlib${ii}$_ilauplo( 'L' ) ) then
info = 1_${ik}$
else if( n<0_${ik}$ )then
info = 2_${ik}$
else if( lda<max( 1_${ik}$, n ) )then
info = 5_${ik}$
else if( incx==0_${ik}$ )then
info = 7_${ik}$
else if( incy==0_${ik}$ )then
info = 10_${ik}$
end if
if( info/=0_${ik}$ )then
call stdlib${ii}$_xerbla( 'SLA_SYAMV', info )
return
end if
! quick return if possible.
if( ( n==0 ).or.( ( alpha==zero ).and.( beta==one ) ) )return
! set up the start points in x and y.
if( incx>0_${ik}$ )then
kx = 1_${ik}$
else
kx = 1_${ik}$ - ( n - 1_${ik}$ )*incx
end if
if( incy>0_${ik}$ )then
ky = 1_${ik}$
else
ky = 1_${ik}$ - ( n - 1_${ik}$ )*incy
end if
! set safe1 essentially to be the underflow threshold times the
! number of additions in each row.
safe1 = stdlib${ii}$_slamch( 'SAFE MINIMUM' )
safe1 = (n+1)*safe1
! form y := alpha*abs(a)*abs(x) + beta*abs(y).
! the o(n^2) symb_zero tests could be replaced by o(n) queries to
! the inexact flag. still doesn't help change the iteration order
! to per-column.
iy = ky
if ( incx==1_${ik}$ ) then
if ( uplo == stdlib${ii}$_ilauplo( 'U' ) ) then
do i = 1, n
if ( beta == zero ) then
symb_zero = .true.
y( iy ) = zero
else if ( y( iy ) == zero ) then
symb_zero = .true.
else
symb_zero = .false.
y( iy ) = beta * abs( y( iy ) )
end if
if ( alpha /= zero ) then
do j = 1, i
temp = abs( a( j, i ) )
symb_zero = symb_zero .and.( x( j ) == zero .or. temp == zero )
y( iy ) = y( iy ) + alpha*abs( x( j ) )*temp
end do
do j = i+1, n
temp = abs( a( i, j ) )
symb_zero = symb_zero .and.( x( j ) == zero .or. temp == zero )
y( iy ) = y( iy ) + alpha*abs( x( j ) )*temp
end do
end if
if ( .not.symb_zero )y( iy ) = y( iy ) + sign( safe1, y( iy ) )
iy = iy + incy
end do
else
do i = 1, n
if ( beta == zero ) then
symb_zero = .true.
y( iy ) = zero
else if ( y( iy ) == zero ) then
symb_zero = .true.
else
symb_zero = .false.
y( iy ) = beta * abs( y( iy ) )
end if
if ( alpha /= zero ) then
do j = 1, i
temp = abs( a( i, j ) )
symb_zero = symb_zero .and.( x( j ) == zero .or. temp == zero )
y( iy ) = y( iy ) + alpha*abs( x( j ) )*temp
end do
do j = i+1, n
temp = abs( a( j, i ) )
symb_zero = symb_zero .and.( x( j ) == zero .or. temp == zero )
y( iy ) = y( iy ) + alpha*abs( x( j ) )*temp
end do
end if
if ( .not.symb_zero )y( iy ) = y( iy ) + sign( safe1, y( iy ) )
iy = iy + incy
end do
end if
else
if ( uplo == stdlib${ii}$_ilauplo( 'U' ) ) then
do i = 1, n
if ( beta == zero ) then
symb_zero = .true.
y( iy ) = zero
else if ( y( iy ) == zero ) then
symb_zero = .true.
else
symb_zero = .false.
y( iy ) = beta * abs( y( iy ) )
end if
jx = kx
if ( alpha /= zero ) then
do j = 1, i
temp = abs( a( j, i ) )
symb_zero = symb_zero .and.( x( j ) == zero .or. temp == zero )
y( iy ) = y( iy ) + alpha*abs( x( jx ) )*temp
jx = jx + incx
end do
do j = i+1, n
temp = abs( a( i, j ) )
symb_zero = symb_zero .and.( x( j ) == zero .or. temp == zero )
y( iy ) = y( iy ) + alpha*abs( x( jx ) )*temp
jx = jx + incx
end do
end if
if ( .not.symb_zero )y( iy ) = y( iy ) + sign( safe1, y( iy ) )
iy = iy + incy
end do
else
do i = 1, n
if ( beta == zero ) then
symb_zero = .true.
y( iy ) = zero
else if ( y( iy ) == zero ) then
symb_zero = .true.
else
symb_zero = .false.
y( iy ) = beta * abs( y( iy ) )
end if
jx = kx
if ( alpha /= zero ) then
do j = 1, i
temp = abs( a( i, j ) )
symb_zero = symb_zero .and.( x( j ) == zero .or. temp == zero )
y( iy ) = y( iy ) + alpha*abs( x( jx ) )*temp
jx = jx + incx
end do
do j = i+1, n
temp = abs( a( j, i ) )
symb_zero = symb_zero .and.( x( j ) == zero .or. temp == zero )
y( iy ) = y( iy ) + alpha*abs( x( jx ) )*temp
jx = jx + incx
end do
end if
if ( .not.symb_zero )y( iy ) = y( iy ) + sign( safe1, y( iy ) )
iy = iy + incy
end do
end if
end if
return
end subroutine stdlib${ii}$_sla_syamv
module subroutine stdlib${ii}$_dla_syamv( uplo, n, alpha, a, lda, x, incx, beta, y,incy )
!! DLA_SYAMV performs the matrix-vector operation
!! y := alpha*abs(A)*abs(x) + beta*abs(y),
!! where alpha and beta are scalars, x and y are vectors and A is an
!! n by n symmetric matrix.
!! This function is primarily used in calculating error bounds.
!! To protect against underflow during evaluation, components in
!! the resulting vector are perturbed away from zero by (N+1)
!! times the underflow threshold. To prevent unnecessarily large
!! errors for block-structure embedded in general matrices,
!! "symbolically" zero components are not perturbed. A zero
!! entry is considered "symbolic" if all multiplications involved
!! in computing that entry have at least one zero multiplicand.
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
real(dp), intent(in) :: alpha, beta
integer(${ik}$), intent(in) :: incx, incy, lda, n, uplo
! Array Arguments
real(dp), intent(in) :: a(lda,*), x(*)
real(dp), intent(inout) :: y(*)
! =====================================================================
! Local Scalars
logical(lk) :: symb_zero
real(dp) :: temp, safe1
integer(${ik}$) :: i, info, iy, j, jx, kx, ky
! Intrinsic Functions
! Executable Statements
! test the input parameters.
info = 0_${ik}$
if ( uplo/=stdlib${ii}$_ilauplo( 'U' ) .and.uplo/=stdlib${ii}$_ilauplo( 'L' ) ) then
info = 1_${ik}$
else if( n<0_${ik}$ )then
info = 2_${ik}$
else if( lda<max( 1_${ik}$, n ) )then
info = 5_${ik}$
else if( incx==0_${ik}$ )then
info = 7_${ik}$
else if( incy==0_${ik}$ )then
info = 10_${ik}$
end if
if( info/=0_${ik}$ )then
call stdlib${ii}$_xerbla( 'DLA_SYAMV', info )
return
end if
! quick return if possible.
if( ( n==0 ).or.( ( alpha==zero ).and.( beta==one ) ) )return
! set up the start points in x and y.
if( incx>0_${ik}$ )then
kx = 1_${ik}$
else
kx = 1_${ik}$ - ( n - 1_${ik}$ )*incx
end if
if( incy>0_${ik}$ )then
ky = 1_${ik}$
else
ky = 1_${ik}$ - ( n - 1_${ik}$ )*incy
end if
! set safe1 essentially to be the underflow threshold times the
! number of additions in each row.
safe1 = stdlib${ii}$_dlamch( 'SAFE MINIMUM' )
safe1 = (n+1)*safe1
! form y := alpha*abs(a)*abs(x) + beta*abs(y).
! the o(n^2) symb_zero tests could be replaced by o(n) queries to
! the inexact flag. still doesn't help change the iteration order
! to per-column.
iy = ky
if ( incx==1_${ik}$ ) then
if ( uplo == stdlib${ii}$_ilauplo( 'U' ) ) then
do i = 1, n
if ( beta == zero ) then
symb_zero = .true.
y( iy ) = zero
else if ( y( iy ) == zero ) then
symb_zero = .true.
else
symb_zero = .false.
y( iy ) = beta * abs( y( iy ) )
end if
if ( alpha /= zero ) then
do j = 1, i
temp = abs( a( j, i ) )
symb_zero = symb_zero .and.( x( j ) == zero .or. temp == zero )
y( iy ) = y( iy ) + alpha*abs( x( j ) )*temp
end do
do j = i+1, n
temp = abs( a( i, j ) )
symb_zero = symb_zero .and.( x( j ) == zero .or. temp == zero )
y( iy ) = y( iy ) + alpha*abs( x( j ) )*temp
end do
end if
if ( .not.symb_zero )y( iy ) = y( iy ) + sign( safe1, y( iy ) )
iy = iy + incy
end do
else
do i = 1, n
if ( beta == zero ) then
symb_zero = .true.
y( iy ) = zero
else if ( y( iy ) == zero ) then
symb_zero = .true.
else
symb_zero = .false.
y( iy ) = beta * abs( y( iy ) )
end if
if ( alpha /= zero ) then
do j = 1, i
temp = abs( a( i, j ) )
symb_zero = symb_zero .and.( x( j ) == zero .or. temp == zero )
y( iy ) = y( iy ) + alpha*abs( x( j ) )*temp
end do
do j = i+1, n
temp = abs( a( j, i ) )
symb_zero = symb_zero .and.( x( j ) == zero .or. temp == zero )
y( iy ) = y( iy ) + alpha*abs( x( j ) )*temp
end do
end if
if ( .not.symb_zero )y( iy ) = y( iy ) + sign( safe1, y( iy ) )
iy = iy + incy
end do
end if
else
if ( uplo == stdlib${ii}$_ilauplo( 'U' ) ) then
do i = 1, n
if ( beta == zero ) then
symb_zero = .true.
y( iy ) = zero
else if ( y( iy ) == zero ) then
symb_zero = .true.
else
symb_zero = .false.
y( iy ) = beta * abs( y( iy ) )
end if
jx = kx
if ( alpha /= zero ) then
do j = 1, i
temp = abs( a( j, i ) )
symb_zero = symb_zero .and.( x( j ) == zero .or. temp == zero )
y( iy ) = y( iy ) + alpha*abs( x( jx ) )*temp
jx = jx + incx
end do
do j = i+1, n
temp = abs( a( i, j ) )
symb_zero = symb_zero .and.( x( j ) == zero .or. temp == zero )
y( iy ) = y( iy ) + alpha*abs( x( jx ) )*temp
jx = jx + incx
end do
end if
if ( .not.symb_zero )y( iy ) = y( iy ) + sign( safe1, y( iy ) )
iy = iy + incy
end do
else
do i = 1, n
if ( beta == zero ) then
symb_zero = .true.
y( iy ) = zero
else if ( y( iy ) == zero ) then
symb_zero = .true.
else
symb_zero = .false.
y( iy ) = beta * abs( y( iy ) )
end if
jx = kx
if ( alpha /= zero ) then
do j = 1, i
temp = abs( a( i, j ) )
symb_zero = symb_zero .and.( x( j ) == zero .or. temp == zero )
y( iy ) = y( iy ) + alpha*abs( x( jx ) )*temp
jx = jx + incx
end do
do j = i+1, n
temp = abs( a( j, i ) )
symb_zero = symb_zero .and.( x( j ) == zero .or. temp == zero )
y( iy ) = y( iy ) + alpha*abs( x( jx ) )*temp
jx = jx + incx
end do
end if
if ( .not.symb_zero )y( iy ) = y( iy ) + sign( safe1, y( iy ) )
iy = iy + incy
end do
end if
end if
return
end subroutine stdlib${ii}$_dla_syamv
#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
module subroutine stdlib${ii}$_${ri}$la_syamv( uplo, n, alpha, a, lda, x, incx, beta, y,incy )
!! DLA_SYAMV: performs the matrix-vector operation
!! y := alpha*abs(A)*abs(x) + beta*abs(y),
!! where alpha and beta are scalars, x and y are vectors and A is an
!! n by n symmetric matrix.
!! This function is primarily used in calculating error bounds.
!! To protect against underflow during evaluation, components in
!! the resulting vector are perturbed away from zero by (N+1)
!! times the underflow threshold. To prevent unnecessarily large
!! errors for block-structure embedded in general matrices,
!! "symbolically" zero components are not perturbed. A zero
!! entry is considered "symbolic" if all multiplications involved
!! in computing that entry have at least one zero multiplicand.
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_${rk}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
real(${rk}$), intent(in) :: alpha, beta
integer(${ik}$), intent(in) :: incx, incy, lda, n, uplo
! Array Arguments
real(${rk}$), intent(in) :: a(lda,*), x(*)
real(${rk}$), intent(inout) :: y(*)
! =====================================================================
! Local Scalars
logical(lk) :: symb_wero
real(${rk}$) :: temp, safe1
integer(${ik}$) :: i, info, iy, j, jx, kx, ky
! Intrinsic Functions
! Executable Statements
! test the input parameters.
info = 0_${ik}$
if ( uplo/=stdlib${ii}$_ilauplo( 'U' ) .and.uplo/=stdlib${ii}$_ilauplo( 'L' ) ) then
info = 1_${ik}$
else if( n<0_${ik}$ )then
info = 2_${ik}$
else if( lda<max( 1_${ik}$, n ) )then
info = 5_${ik}$
else if( incx==0_${ik}$ )then
info = 7_${ik}$
else if( incy==0_${ik}$ )then
info = 10_${ik}$
end if
if( info/=0_${ik}$ )then
call stdlib${ii}$_xerbla( 'DLA_SYAMV', info )
return
end if
! quick return if possible.
if( ( n==0 ).or.( ( alpha==zero ).and.( beta==one ) ) )return
! set up the start points in x and y.
if( incx>0_${ik}$ )then
kx = 1_${ik}$
else
kx = 1_${ik}$ - ( n - 1_${ik}$ )*incx
end if
if( incy>0_${ik}$ )then
ky = 1_${ik}$
else
ky = 1_${ik}$ - ( n - 1_${ik}$ )*incy
end if
! set safe1 essentially to be the underflow threshold times the
! number of additions in each row.
safe1 = stdlib${ii}$_${ri}$lamch( 'SAFE MINIMUM' )
safe1 = (n+1)*safe1
! form y := alpha*abs(a)*abs(x) + beta*abs(y).
! the o(n^2) symb_wero tests could be replaced by o(n) queries to
! the inexact flag. still doesn't help change the iteration order
! to per-column.
iy = ky
if ( incx==1_${ik}$ ) then
if ( uplo == stdlib${ii}$_ilauplo( 'U' ) ) then
do i = 1, n
if ( beta == zero ) then
symb_wero = .true.
y( iy ) = zero
else if ( y( iy ) == zero ) then
symb_wero = .true.
else
symb_wero = .false.
y( iy ) = beta * abs( y( iy ) )
end if
if ( alpha /= zero ) then
do j = 1, i
temp = abs( a( j, i ) )
symb_wero = symb_wero .and.( x( j ) == zero .or. temp == zero )
y( iy ) = y( iy ) + alpha*abs( x( j ) )*temp
end do
do j = i+1, n
temp = abs( a( i, j ) )
symb_wero = symb_wero .and.( x( j ) == zero .or. temp == zero )
y( iy ) = y( iy ) + alpha*abs( x( j ) )*temp
end do
end if
if ( .not.symb_wero )y( iy ) = y( iy ) + sign( safe1, y( iy ) )
iy = iy + incy
end do
else
do i = 1, n
if ( beta == zero ) then
symb_wero = .true.
y( iy ) = zero
else if ( y( iy ) == zero ) then
symb_wero = .true.
else
symb_wero = .false.
y( iy ) = beta * abs( y( iy ) )
end if
if ( alpha /= zero ) then
do j = 1, i
temp = abs( a( i, j ) )
symb_wero = symb_wero .and.( x( j ) == zero .or. temp == zero )
y( iy ) = y( iy ) + alpha*abs( x( j ) )*temp
end do
do j = i+1, n
temp = abs( a( j, i ) )
symb_wero = symb_wero .and.( x( j ) == zero .or. temp == zero )
y( iy ) = y( iy ) + alpha*abs( x( j ) )*temp
end do
end if
if ( .not.symb_wero )y( iy ) = y( iy ) + sign( safe1, y( iy ) )
iy = iy + incy
end do
end if
else
if ( uplo == stdlib${ii}$_ilauplo( 'U' ) ) then
do i = 1, n
if ( beta == zero ) then
symb_wero = .true.
y( iy ) = zero
else if ( y( iy ) == zero ) then
symb_wero = .true.
else
symb_wero = .false.
y( iy ) = beta * abs( y( iy ) )
end if
jx = kx
if ( alpha /= zero ) then
do j = 1, i
temp = abs( a( j, i ) )
symb_wero = symb_wero .and.( x( j ) == zero .or. temp == zero )
y( iy ) = y( iy ) + alpha*abs( x( jx ) )*temp
jx = jx + incx
end do
do j = i+1, n
temp = abs( a( i, j ) )
symb_wero = symb_wero .and.( x( j ) == zero .or. temp == zero )
y( iy ) = y( iy ) + alpha*abs( x( jx ) )*temp
jx = jx + incx
end do
end if
if ( .not.symb_wero )y( iy ) = y( iy ) + sign( safe1, y( iy ) )
iy = iy + incy
end do
else
do i = 1, n
if ( beta == zero ) then
symb_wero = .true.
y( iy ) = zero
else if ( y( iy ) == zero ) then
symb_wero = .true.
else
symb_wero = .false.
y( iy ) = beta * abs( y( iy ) )
end if
jx = kx
if ( alpha /= zero ) then
do j = 1, i
temp = abs( a( i, j ) )
symb_wero = symb_wero .and.( x( j ) == zero .or. temp == zero )
y( iy ) = y( iy ) + alpha*abs( x( jx ) )*temp
jx = jx + incx
end do
do j = i+1, n
temp = abs( a( j, i ) )
symb_wero = symb_wero .and.( x( j ) == zero .or. temp == zero )
y( iy ) = y( iy ) + alpha*abs( x( jx ) )*temp
jx = jx + incx
end do
end if
if ( .not.symb_wero )y( iy ) = y( iy ) + sign( safe1, y( iy ) )
iy = iy + incy
end do
end if
end if
return
end subroutine stdlib${ii}$_${ri}$la_syamv
#:endif
#:endfor
module subroutine stdlib${ii}$_cla_syamv( uplo, n, alpha, a, lda, x, incx, beta, y,incy )
!! CLA_SYAMV performs the matrix-vector operation
!! y := alpha*abs(A)*abs(x) + beta*abs(y),
!! where alpha and beta are scalars, x and y are vectors and A is an
!! n by n symmetric matrix.
!! This function is primarily used in calculating error bounds.
!! To protect against underflow during evaluation, components in
!! the resulting vector are perturbed away from zero by (N+1)
!! times the underflow threshold. To prevent unnecessarily large
!! errors for block-structure embedded in general matrices,
!! "symbolically" zero components are not perturbed. A zero
!! entry is considered "symbolic" if all multiplications involved
!! in computing that entry have at least one zero multiplicand.
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
real(sp), intent(in) :: alpha, beta
integer(${ik}$), intent(in) :: incx, incy, lda, n
integer(${ik}$), intent(in) :: uplo
! Array Arguments
complex(sp), intent(in) :: a(lda,*), x(*)
real(sp), intent(inout) :: y(*)
! =====================================================================
! Local Scalars
logical(lk) :: symb_zero
real(sp) :: temp, safe1
integer(${ik}$) :: i, info, iy, j, jx, kx, ky
complex(sp) :: zdum
! Intrinsic Functions
! Statement Functions
real(sp) :: cabs1
! Statement Function Definitions
cabs1( zdum ) = abs( real( zdum,KIND=sp) ) + abs( aimag ( zdum ) )
! Executable Statements
! test the input parameters.
info = 0_${ik}$
if ( uplo/=stdlib${ii}$_ilauplo( 'U' ) .and.uplo/=stdlib${ii}$_ilauplo( 'L' ) )then
info = 1_${ik}$
else if( n<0_${ik}$ )then
info = 2_${ik}$
else if( lda<max( 1_${ik}$, n ) )then
info = 5_${ik}$
else if( incx==0_${ik}$ )then
info = 7_${ik}$
else if( incy==0_${ik}$ )then
info = 10_${ik}$
end if
if( info/=0_${ik}$ )then
call stdlib${ii}$_xerbla( 'CLA_SYAMV', info )
return
end if
! quick return if possible.
if( ( n==0 ).or.( ( alpha==zero ).and.( beta==one ) ) )return
! set up the start points in x and y.
if( incx>0_${ik}$ )then
kx = 1_${ik}$
else
kx = 1_${ik}$ - ( n - 1_${ik}$ )*incx
end if
if( incy>0_${ik}$ )then
ky = 1_${ik}$
else
ky = 1_${ik}$ - ( n - 1_${ik}$ )*incy
end if
! set safe1 essentially to be the underflow threshold times the
! number of additions in each row.
safe1 = stdlib${ii}$_slamch( 'SAFE MINIMUM' )
safe1 = (n+1)*safe1
! form y := alpha*abs(a)*abs(x) + beta*abs(y).
! the o(n^2) symb_zero tests could be replaced by o(n) queries to
! the inexact flag. still doesn't help change the iteration order
! to per-column.
iy = ky
if ( incx==1_${ik}$ ) then
if ( uplo == stdlib${ii}$_ilauplo( 'U' ) ) then
do i = 1, n
if ( beta == zero ) then
symb_zero = .true.
y( iy ) = zero
else if ( y( iy ) == zero ) then
symb_zero = .true.
else
symb_zero = .false.
y( iy ) = beta * abs( y( iy ) )
end if
if ( alpha /= zero ) then
do j = 1, i
temp = cabs1( a( j, i ) )
symb_zero = symb_zero .and.( x( j ) == zero .or. temp == zero )
y( iy ) = y( iy ) + alpha*cabs1( x( j ) )*temp
end do
do j = i+1, n
temp = cabs1( a( i, j ) )
symb_zero = symb_zero .and.( x( j ) == zero .or. temp == zero )
y( iy ) = y( iy ) + alpha*cabs1( x( j ) )*temp
end do
end if
if ( .not.symb_zero )y( iy ) = y( iy ) + sign( safe1, y( iy ) )
iy = iy + incy
end do
else
do i = 1, n
if ( beta == zero ) then
symb_zero = .true.
y( iy ) = zero
else if ( y( iy ) == zero ) then
symb_zero = .true.
else
symb_zero = .false.
y( iy ) = beta * abs( y( iy ) )
end if
if ( alpha /= zero ) then
do j = 1, i
temp = cabs1( a( i, j ) )
symb_zero = symb_zero .and.( x( j ) == zero .or. temp == zero )
y( iy ) = y( iy ) + alpha*cabs1( x( j ) )*temp
end do
do j = i+1, n
temp = cabs1( a( j, i ) )
symb_zero = symb_zero .and.( x( j ) == zero .or. temp == zero )
y( iy ) = y( iy ) + alpha*cabs1( x( j ) )*temp
end do
end if
if ( .not.symb_zero )y( iy ) = y( iy ) + sign( safe1, y( iy ) )
iy = iy + incy
end do
end if
else
if ( uplo == stdlib${ii}$_ilauplo( 'U' ) ) then
do i = 1, n
if ( beta == zero ) then
symb_zero = .true.
y( iy ) = zero
else if ( y( iy ) == zero ) then
symb_zero = .true.
else
symb_zero = .false.
y( iy ) = beta * abs( y( iy ) )
end if
jx = kx
if ( alpha /= zero ) then
do j = 1, i
temp = cabs1( a( j, i ) )
symb_zero = symb_zero .and.( x( j ) == zero .or. temp == zero )
y( iy ) = y( iy ) + alpha*cabs1( x( jx ) )*temp
jx = jx + incx
end do
do j = i+1, n
temp = cabs1( a( i, j ) )
symb_zero = symb_zero .and.( x( j ) == zero .or. temp == zero )
y( iy ) = y( iy ) + alpha*cabs1( x( jx ) )*temp
jx = jx + incx
end do
end if
if ( .not.symb_zero )y( iy ) = y( iy ) + sign( safe1, y( iy ) )
iy = iy + incy
end do
else
do i = 1, n
if ( beta == zero ) then
symb_zero = .true.
y( iy ) = zero
else if ( y( iy ) == zero ) then
symb_zero = .true.
else
symb_zero = .false.
y( iy ) = beta * abs( y( iy ) )
end if
jx = kx
if ( alpha /= zero ) then
do j = 1, i
temp = cabs1( a( i, j ) )
symb_zero = symb_zero .and.( x( j ) == zero .or. temp == zero )
y( iy ) = y( iy ) + alpha*cabs1( x( jx ) )*temp
jx = jx + incx
end do
do j = i+1, n
temp = cabs1( a( j, i ) )
symb_zero = symb_zero .and.( x( j ) == zero .or. temp == zero )
y( iy ) = y( iy ) + alpha*cabs1( x( jx ) )*temp
jx = jx + incx
end do
end if
if ( .not.symb_zero )y( iy ) = y( iy ) + sign( safe1, y( iy ) )
iy = iy + incy
end do
end if
end if
return
end subroutine stdlib${ii}$_cla_syamv
module subroutine stdlib${ii}$_zla_syamv( uplo, n, alpha, a, lda, x, incx, beta, y,incy )
!! ZLA_SYAMV performs the matrix-vector operation
!! y := alpha*abs(A)*abs(x) + beta*abs(y),
!! where alpha and beta are scalars, x and y are vectors and A is an
!! n by n symmetric matrix.
!! This function is primarily used in calculating error bounds.
!! To protect against underflow during evaluation, components in
!! the resulting vector are perturbed away from zero by (N+1)
!! times the underflow threshold. To prevent unnecessarily large
!! errors for block-structure embedded in general matrices,
!! "symbolically" zero components are not perturbed. A zero
!! entry is considered "symbolic" if all multiplications involved
!! in computing that entry have at least one zero multiplicand.
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
real(dp), intent(in) :: alpha, beta
integer(${ik}$), intent(in) :: incx, incy, lda, n
integer(${ik}$), intent(in) :: uplo
! Array Arguments
complex(dp), intent(in) :: a(lda,*), x(*)
real(dp), intent(inout) :: y(*)
! =====================================================================
! Local Scalars
logical(lk) :: symb_zero
real(dp) :: temp, safe1
integer(${ik}$) :: i, info, iy, j, jx, kx, ky
complex(dp) :: zdum
! Intrinsic Functions
! Statement Functions
real(dp) :: cabs1
! Statement Function Definitions
cabs1( zdum ) = abs( real( zdum,KIND=dp) ) + abs( aimag ( zdum ) )
! Executable Statements
! test the input parameters.
info = 0_${ik}$
if ( uplo/=stdlib${ii}$_ilauplo( 'U' ) .and.uplo/=stdlib${ii}$_ilauplo( 'L' ) )then
info = 1_${ik}$
else if( n<0_${ik}$ )then
info = 2_${ik}$
else if( lda<max( 1_${ik}$, n ) )then
info = 5_${ik}$
else if( incx==0_${ik}$ )then
info = 7_${ik}$
else if( incy==0_${ik}$ )then
info = 10_${ik}$
end if
if( info/=0_${ik}$ )then
call stdlib${ii}$_xerbla( 'ZLA_SYAMV', info )
return
end if
! quick return if possible.
if( ( n==0 ).or.( ( alpha==zero ).and.( beta==one ) ) )return
! set up the start points in x and y.
if( incx>0_${ik}$ )then
kx = 1_${ik}$
else
kx = 1_${ik}$ - ( n - 1_${ik}$ )*incx
end if
if( incy>0_${ik}$ )then
ky = 1_${ik}$
else
ky = 1_${ik}$ - ( n - 1_${ik}$ )*incy
end if
! set safe1 essentially to be the underflow threshold times the
! number of additions in each row.
safe1 = stdlib${ii}$_dlamch( 'SAFE MINIMUM' )
safe1 = (n+1)*safe1
! form y := alpha*abs(a)*abs(x) + beta*abs(y).
! the o(n^2) symb_zero tests could be replaced by o(n) queries to
! the inexact flag. still doesn't help change the iteration order
! to per-column.
iy = ky
if ( incx==1_${ik}$ ) then
if ( uplo == stdlib${ii}$_ilauplo( 'U' ) ) then
do i = 1, n
if ( beta == zero ) then
symb_zero = .true.
y( iy ) = zero
else if ( y( iy ) == zero ) then
symb_zero = .true.
else
symb_zero = .false.
y( iy ) = beta * abs( y( iy ) )
end if
if ( alpha /= zero ) then
do j = 1, i
temp = cabs1( a( j, i ) )
symb_zero = symb_zero .and.( x( j ) == zero .or. temp == zero )
y( iy ) = y( iy ) + alpha*cabs1( x( j ) )*temp
end do
do j = i+1, n
temp = cabs1( a( i, j ) )
symb_zero = symb_zero .and.( x( j ) == zero .or. temp == zero )
y( iy ) = y( iy ) + alpha*cabs1( x( j ) )*temp
end do
end if
if ( .not.symb_zero )y( iy ) = y( iy ) + sign( safe1, y( iy ) )
iy = iy + incy
end do
else
do i = 1, n
if ( beta == zero ) then
symb_zero = .true.
y( iy ) = zero
else if ( y( iy ) == zero ) then
symb_zero = .true.
else
symb_zero = .false.
y( iy ) = beta * abs( y( iy ) )
end if
if ( alpha /= zero ) then
do j = 1, i
temp = cabs1( a( i, j ) )
symb_zero = symb_zero .and.( x( j ) == zero .or. temp == zero )
y( iy ) = y( iy ) + alpha*cabs1( x( j ) )*temp
end do
do j = i+1, n
temp = cabs1( a( j, i ) )
symb_zero = symb_zero .and.( x( j ) == zero .or. temp == zero )
y( iy ) = y( iy ) + alpha*cabs1( x( j ) )*temp
end do
end if
if ( .not.symb_zero )y( iy ) = y( iy ) + sign( safe1, y( iy ) )
iy = iy + incy
end do
end if
else
if ( uplo == stdlib${ii}$_ilauplo( 'U' ) ) then
do i = 1, n
if ( beta == zero ) then
symb_zero = .true.
y( iy ) = zero
else if ( y( iy ) == zero ) then
symb_zero = .true.
else
symb_zero = .false.
y( iy ) = beta * abs( y( iy ) )
end if
jx = kx
if ( alpha /= zero ) then
do j = 1, i
temp = cabs1( a( j, i ) )
symb_zero = symb_zero .and.( x( j ) == zero .or. temp == zero )
y( iy ) = y( iy ) + alpha*cabs1( x( jx ) )*temp
jx = jx + incx
end do
do j = i+1, n
temp = cabs1( a( i, j ) )
symb_zero = symb_zero .and.( x( j ) == zero .or. temp == zero )
y( iy ) = y( iy ) + alpha*cabs1( x( jx ) )*temp
jx = jx + incx
end do
end if
if ( .not.symb_zero )y( iy ) = y( iy ) + sign( safe1, y( iy ) )
iy = iy + incy
end do
else
do i = 1, n
if ( beta == zero ) then
symb_zero = .true.
y( iy ) = zero
else if ( y( iy ) == zero ) then
symb_zero = .true.
else
symb_zero = .false.
y( iy ) = beta * abs( y( iy ) )
end if
jx = kx
if ( alpha /= zero ) then
do j = 1, i
temp = cabs1( a( i, j ) )
symb_zero = symb_zero .and.( x( j ) == zero .or. temp == zero )
y( iy ) = y( iy ) + alpha*cabs1( x( jx ) )*temp
jx = jx + incx
end do
do j = i+1, n
temp = cabs1( a( j, i ) )
symb_zero = symb_zero .and.( x( j ) == zero .or. temp == zero )
y( iy ) = y( iy ) + alpha*cabs1( x( jx ) )*temp
jx = jx + incx
end do
end if
if ( .not.symb_zero )y( iy ) = y( iy ) + sign( safe1, y( iy ) )
iy = iy + incy
end do
end if
end if
return
end subroutine stdlib${ii}$_zla_syamv
#:for ck,ct,ci in CMPLX_KINDS_TYPES
#:if not ck in ["sp","dp"]
module subroutine stdlib${ii}$_${ci}$la_syamv( uplo, n, alpha, a, lda, x, incx, beta, y,incy )
!! ZLA_SYAMV: performs the matrix-vector operation
!! y := alpha*abs(A)*abs(x) + beta*abs(y),
!! where alpha and beta are scalars, x and y are vectors and A is an
!! n by n symmetric matrix.
!! This function is primarily used in calculating error bounds.
!! To protect against underflow during evaluation, components in
!! the resulting vector are perturbed away from zero by (N+1)
!! times the underflow threshold. To prevent unnecessarily large
!! errors for block-structure embedded in general matrices,
!! "symbolically" zero components are not perturbed. A zero
!! entry is considered "symbolic" if all multiplications involved
!! in computing that entry have at least one zero multiplicand.
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_${ck}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
real(${ck}$), intent(in) :: alpha, beta
integer(${ik}$), intent(in) :: incx, incy, lda, n
integer(${ik}$), intent(in) :: uplo
! Array Arguments
complex(${ck}$), intent(in) :: a(lda,*), x(*)
real(${ck}$), intent(inout) :: y(*)
! =====================================================================
! Local Scalars
logical(lk) :: symb_wero
real(${ck}$) :: temp, safe1
integer(${ik}$) :: i, info, iy, j, jx, kx, ky
complex(${ck}$) :: zdum
! Intrinsic Functions
! Statement Functions
real(${ck}$) :: cabs1
! Statement Function Definitions
cabs1( zdum ) = abs( real( zdum,KIND=${ck}$) ) + abs( aimag ( zdum ) )
! Executable Statements
! test the input parameters.
info = 0_${ik}$
if ( uplo/=stdlib${ii}$_ilauplo( 'U' ) .and.uplo/=stdlib${ii}$_ilauplo( 'L' ) )then
info = 1_${ik}$
else if( n<0_${ik}$ )then
info = 2_${ik}$
else if( lda<max( 1_${ik}$, n ) )then
info = 5_${ik}$
else if( incx==0_${ik}$ )then
info = 7_${ik}$
else if( incy==0_${ik}$ )then
info = 10_${ik}$
end if
if( info/=0_${ik}$ )then
call stdlib${ii}$_xerbla( 'ZLA_SYAMV', info )
return
end if
! quick return if possible.
if( ( n==0 ).or.( ( alpha==zero ).and.( beta==one ) ) )return
! set up the start points in x and y.
if( incx>0_${ik}$ )then
kx = 1_${ik}$
else
kx = 1_${ik}$ - ( n - 1_${ik}$ )*incx
end if
if( incy>0_${ik}$ )then
ky = 1_${ik}$
else
ky = 1_${ik}$ - ( n - 1_${ik}$ )*incy
end if
! set safe1 essentially to be the underflow threshold times the
! number of additions in each row.
safe1 = stdlib${ii}$_${c2ri(ci)}$lamch( 'SAFE MINIMUM' )
safe1 = (n+1)*safe1
! form y := alpha*abs(a)*abs(x) + beta*abs(y).
! the o(n^2) symb_wero tests could be replaced by o(n) queries to
! the inexact flag. still doesn't help change the iteration order
! to per-column.
iy = ky
if ( incx==1_${ik}$ ) then
if ( uplo == stdlib${ii}$_ilauplo( 'U' ) ) then
do i = 1, n
if ( beta == zero ) then
symb_wero = .true.
y( iy ) = zero
else if ( y( iy ) == zero ) then
symb_wero = .true.
else
symb_wero = .false.
y( iy ) = beta * abs( y( iy ) )
end if
if ( alpha /= zero ) then
do j = 1, i
temp = cabs1( a( j, i ) )
symb_wero = symb_wero .and.( x( j ) == zero .or. temp == zero )
y( iy ) = y( iy ) + alpha*cabs1( x( j ) )*temp
end do
do j = i+1, n
temp = cabs1( a( i, j ) )
symb_wero = symb_wero .and.( x( j ) == zero .or. temp == zero )
y( iy ) = y( iy ) + alpha*cabs1( x( j ) )*temp
end do
end if
if ( .not.symb_wero )y( iy ) = y( iy ) + sign( safe1, y( iy ) )
iy = iy + incy
end do
else
do i = 1, n
if ( beta == zero ) then
symb_wero = .true.
y( iy ) = zero
else if ( y( iy ) == zero ) then
symb_wero = .true.
else
symb_wero = .false.
y( iy ) = beta * abs( y( iy ) )
end if
if ( alpha /= zero ) then
do j = 1, i
temp = cabs1( a( i, j ) )
symb_wero = symb_wero .and.( x( j ) == zero .or. temp == zero )
y( iy ) = y( iy ) + alpha*cabs1( x( j ) )*temp
end do
do j = i+1, n
temp = cabs1( a( j, i ) )
symb_wero = symb_wero .and.( x( j ) == zero .or. temp == zero )
y( iy ) = y( iy ) + alpha*cabs1( x( j ) )*temp
end do
end if
if ( .not.symb_wero )y( iy ) = y( iy ) + sign( safe1, y( iy ) )
iy = iy + incy
end do
end if
else
if ( uplo == stdlib${ii}$_ilauplo( 'U' ) ) then
do i = 1, n
if ( beta == zero ) then
symb_wero = .true.
y( iy ) = zero
else if ( y( iy ) == zero ) then
symb_wero = .true.
else
symb_wero = .false.
y( iy ) = beta * abs( y( iy ) )
end if
jx = kx
if ( alpha /= zero ) then
do j = 1, i
temp = cabs1( a( j, i ) )
symb_wero = symb_wero .and.( x( j ) == zero .or. temp == zero )
y( iy ) = y( iy ) + alpha*cabs1( x( jx ) )*temp
jx = jx + incx
end do
do j = i+1, n
temp = cabs1( a( i, j ) )
symb_wero = symb_wero .and.( x( j ) == zero .or. temp == zero )
y( iy ) = y( iy ) + alpha*cabs1( x( jx ) )*temp
jx = jx + incx
end do
end if
if ( .not.symb_wero )y( iy ) = y( iy ) + sign( safe1, y( iy ) )
iy = iy + incy
end do
else
do i = 1, n
if ( beta == zero ) then
symb_wero = .true.
y( iy ) = zero
else if ( y( iy ) == zero ) then
symb_wero = .true.
else
symb_wero = .false.
y( iy ) = beta * abs( y( iy ) )
end if
jx = kx
if ( alpha /= zero ) then
do j = 1, i
temp = cabs1( a( i, j ) )
symb_wero = symb_wero .and.( x( j ) == zero .or. temp == zero )
y( iy ) = y( iy ) + alpha*cabs1( x( jx ) )*temp
jx = jx + incx
end do
do j = i+1, n
temp = cabs1( a( j, i ) )
symb_wero = symb_wero .and.( x( j ) == zero .or. temp == zero )
y( iy ) = y( iy ) + alpha*cabs1( x( jx ) )*temp
jx = jx + incx
end do
end if
if ( .not.symb_wero )y( iy ) = y( iy ) + sign( safe1, y( iy ) )
iy = iy + incy
end do
end if
end if
return
end subroutine stdlib${ii}$_${ci}$la_syamv
#:endif
#:endfor
real(sp) module function stdlib${ii}$_sla_syrcond( uplo, n, a, lda, af, ldaf, ipiv, cmode,c, info, work, &
!! SLA_SYRCOND estimates the Skeel condition number of op(A) * op2(C)
!! where op2 is determined by CMODE as follows
!! CMODE = 1 op2(C) = C
!! CMODE = 0 op2(C) = I
!! CMODE = -1 op2(C) = inv(C)
!! The Skeel condition number cond(A) = norminf( |inv(A)||A| )
!! is computed by computing scaling factors R such that
!! diag(R)*A*op2(C) is row equilibrated and computing the standard
!! infinity-norm condition number.
iwork )
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: uplo
integer(${ik}$), intent(in) :: n, lda, ldaf, cmode
integer(${ik}$), intent(out) :: info
! Array Arguments
integer(${ik}$), intent(out) :: iwork(*)
integer(${ik}$), intent(in) :: ipiv(*)
real(sp), intent(in) :: a(lda,*), af(ldaf,*), c(*)
real(sp), intent(out) :: work(*)
! =====================================================================
! Local Scalars
character :: normin
integer(${ik}$) :: kase, i, j
real(sp) :: ainvnm, smlnum, tmp
logical(lk) :: up
! Local Arrays
integer(${ik}$) :: isave(3_${ik}$)
! Intrinsic Functions
! Executable Statements
stdlib${ii}$_sla_syrcond = zero
info = 0_${ik}$
if( n<0_${ik}$ ) then
info = -2_${ik}$
else if( lda<max( 1_${ik}$, n ) ) then
info = -4_${ik}$
else if( ldaf<max( 1_${ik}$, n ) ) then
info = -6_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'SLA_SYRCOND', -info )
return
end if
if( n==0_${ik}$ ) then
stdlib${ii}$_sla_syrcond = one
return
end if
up = .false.
if ( stdlib_lsame( uplo, 'U' ) ) up = .true.
! compute the equilibration matrix r such that
! inv(r)*a*c has unit 1-norm.
if ( up ) then
do i = 1, n
tmp = zero
if ( cmode == 1_${ik}$ ) then
do j = 1, i
tmp = tmp + abs( a( j, i ) * c( j ) )
end do
do j = i+1, n
tmp = tmp + abs( a( i, j ) * c( j ) )
end do
else if ( cmode == 0_${ik}$ ) then
do j = 1, i
tmp = tmp + abs( a( j, i ) )
end do
do j = i+1, n
tmp = tmp + abs( a( i, j ) )
end do
else
do j = 1, i
tmp = tmp + abs( a( j, i ) / c( j ) )
end do
do j = i+1, n
tmp = tmp + abs( a( i, j ) / c( j ) )
end do
end if
work( 2_${ik}$*n+i ) = tmp
end do
else
do i = 1, n
tmp = zero
if ( cmode == 1_${ik}$ ) then
do j = 1, i
tmp = tmp + abs( a( i, j ) * c( j ) )
end do
do j = i+1, n
tmp = tmp + abs( a( j, i ) * c( j ) )
end do
else if ( cmode == 0_${ik}$ ) then
do j = 1, i
tmp = tmp + abs( a( i, j ) )
end do
do j = i+1, n
tmp = tmp + abs( a( j, i ) )
end do
else
do j = 1, i
tmp = tmp + abs( a( i, j) / c( j ) )
end do
do j = i+1, n
tmp = tmp + abs( a( j, i) / c( j ) )
end do
end if
work( 2_${ik}$*n+i ) = tmp
end do
endif
! estimate the norm of inv(op(a)).
smlnum = stdlib${ii}$_slamch( 'SAFE MINIMUM' )
ainvnm = zero
normin = 'N'
kase = 0_${ik}$
10 continue
call stdlib${ii}$_slacn2( n, work( n+1 ), work, iwork, ainvnm, kase, isave )
if( kase/=0_${ik}$ ) then
if( kase==2_${ik}$ ) then
! multiply by r.
do i = 1, n
work( i ) = work( i ) * work( 2_${ik}$*n+i )
end do
if ( up ) then
call stdlib${ii}$_ssytrs( 'U', n, 1_${ik}$, af, ldaf, ipiv, work, n, info )
else
call stdlib${ii}$_ssytrs( 'L', n, 1_${ik}$, af, ldaf, ipiv, work, n, info )
endif
! multiply by inv(c).
if ( cmode == 1_${ik}$ ) then
do i = 1, n
work( i ) = work( i ) / c( i )
end do
else if ( cmode == -1_${ik}$ ) then
do i = 1, n
work( i ) = work( i ) * c( i )
end do
end if
else
! multiply by inv(c**t).
if ( cmode == 1_${ik}$ ) then
do i = 1, n
work( i ) = work( i ) / c( i )
end do
else if ( cmode == -1_${ik}$ ) then
do i = 1, n
work( i ) = work( i ) * c( i )
end do
end if
if ( up ) then
call stdlib${ii}$_ssytrs( 'U', n, 1_${ik}$, af, ldaf, ipiv, work, n, info )
else
call stdlib${ii}$_ssytrs( 'L', n, 1_${ik}$, af, ldaf, ipiv, work, n, info )
endif
! multiply by r.
do i = 1, n
work( i ) = work( i ) * work( 2_${ik}$*n+i )
end do
end if
go to 10
end if
! compute the estimate of the reciprocal condition number.
if( ainvnm /= 0.0_sp )stdlib${ii}$_sla_syrcond = ( one / ainvnm )
return
end function stdlib${ii}$_sla_syrcond
real(dp) module function stdlib${ii}$_dla_syrcond( uplo, n, a, lda, af, ldaf,ipiv, cmode, c, info, work,&
!! DLA_SYRCOND estimates the Skeel condition number of op(A) * op2(C)
!! where op2 is determined by CMODE as follows
!! CMODE = 1 op2(C) = C
!! CMODE = 0 op2(C) = I
!! CMODE = -1 op2(C) = inv(C)
!! The Skeel condition number cond(A) = norminf( |inv(A)||A| )
!! is computed by computing scaling factors R such that
!! diag(R)*A*op2(C) is row equilibrated and computing the standard
!! infinity-norm condition number.
iwork )
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: uplo
integer(${ik}$), intent(in) :: n, lda, ldaf, cmode
integer(${ik}$), intent(out) :: info
! Array Arguments
integer(${ik}$), intent(out) :: iwork(*)
integer(${ik}$), intent(in) :: ipiv(*)
real(dp), intent(in) :: a(lda,*), af(ldaf,*), c(*)
real(dp), intent(out) :: work(*)
! =====================================================================
! Local Scalars
character :: normin
integer(${ik}$) :: kase, i, j
real(dp) :: ainvnm, smlnum, tmp
logical(lk) :: up
! Local Arrays
integer(${ik}$) :: isave(3_${ik}$)
! Intrinsic Functions
! Executable Statements
stdlib${ii}$_dla_syrcond = zero
info = 0_${ik}$
if( n<0_${ik}$ ) then
info = -2_${ik}$
else if( lda<max( 1_${ik}$, n ) ) then
info = -4_${ik}$
else if( ldaf<max( 1_${ik}$, n ) ) then
info = -6_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'DLA_SYRCOND', -info )
return
end if
if( n==0_${ik}$ ) then
stdlib${ii}$_dla_syrcond = one
return
end if
up = .false.
if ( stdlib_lsame( uplo, 'U' ) ) up = .true.
! compute the equilibration matrix r such that
! inv(r)*a*c has unit 1-norm.
if ( up ) then
do i = 1, n
tmp = zero
if ( cmode == 1_${ik}$ ) then
do j = 1, i
tmp = tmp + abs( a( j, i ) * c( j ) )
end do
do j = i+1, n
tmp = tmp + abs( a( i, j ) * c( j ) )
end do
else if ( cmode == 0_${ik}$ ) then
do j = 1, i
tmp = tmp + abs( a( j, i ) )
end do
do j = i+1, n
tmp = tmp + abs( a( i, j ) )
end do
else
do j = 1, i
tmp = tmp + abs( a( j, i ) / c( j ) )
end do
do j = i+1, n
tmp = tmp + abs( a( i, j ) / c( j ) )
end do
end if
work( 2_${ik}$*n+i ) = tmp
end do
else
do i = 1, n
tmp = zero
if ( cmode == 1_${ik}$ ) then
do j = 1, i
tmp = tmp + abs( a( i, j ) * c( j ) )
end do
do j = i+1, n
tmp = tmp + abs( a( j, i ) * c( j ) )
end do
else if ( cmode == 0_${ik}$ ) then
do j = 1, i
tmp = tmp + abs( a( i, j ) )
end do
do j = i+1, n
tmp = tmp + abs( a( j, i ) )
end do
else
do j = 1, i
tmp = tmp + abs( a( i, j) / c( j ) )
end do
do j = i+1, n
tmp = tmp + abs( a( j, i) / c( j ) )
end do
end if
work( 2_${ik}$*n+i ) = tmp
end do
endif
! estimate the norm of inv(op(a)).
smlnum = stdlib${ii}$_dlamch( 'SAFE MINIMUM' )
ainvnm = zero
normin = 'N'
kase = 0_${ik}$
10 continue
call stdlib${ii}$_dlacn2( n, work( n+1 ), work, iwork, ainvnm, kase, isave )
if( kase/=0_${ik}$ ) then
if( kase==2_${ik}$ ) then
! multiply by r.
do i = 1, n
work( i ) = work( i ) * work( 2_${ik}$*n+i )
end do
if ( up ) then
call stdlib${ii}$_dsytrs( 'U', n, 1_${ik}$, af, ldaf, ipiv, work, n, info )
else
call stdlib${ii}$_dsytrs( 'L', n, 1_${ik}$, af, ldaf, ipiv, work, n, info )
endif
! multiply by inv(c).
if ( cmode == 1_${ik}$ ) then
do i = 1, n
work( i ) = work( i ) / c( i )
end do
else if ( cmode == -1_${ik}$ ) then
do i = 1, n
work( i ) = work( i ) * c( i )
end do
end if
else
! multiply by inv(c**t).
if ( cmode == 1_${ik}$ ) then
do i = 1, n
work( i ) = work( i ) / c( i )
end do
else if ( cmode == -1_${ik}$ ) then
do i = 1, n
work( i ) = work( i ) * c( i )
end do
end if
if ( up ) then
call stdlib${ii}$_dsytrs( 'U', n, 1_${ik}$, af, ldaf, ipiv, work, n, info )
else
call stdlib${ii}$_dsytrs( 'L', n, 1_${ik}$, af, ldaf, ipiv, work, n, info )
endif
! multiply by r.
do i = 1, n
work( i ) = work( i ) * work( 2_${ik}$*n+i )
end do
end if
go to 10
end if
! compute the estimate of the reciprocal condition number.
if( ainvnm /= zero )stdlib${ii}$_dla_syrcond = ( one / ainvnm )
return
end function stdlib${ii}$_dla_syrcond
#:for rk,rt,ri in REAL_KINDS_TYPES
#:if not rk in ["sp","dp"]
real(${rk}$) module function stdlib${ii}$_${ri}$la_syrcond( uplo, n, a, lda, af, ldaf,ipiv, cmode, c, info, work,&
!! DLA_SYRCOND: estimates the Skeel condition number of op(A) * op2(C)
!! where op2 is determined by CMODE as follows
!! CMODE = 1 op2(C) = C
!! CMODE = 0 op2(C) = I
!! CMODE = -1 op2(C) = inv(C)
!! The Skeel condition number cond(A) = norminf( |inv(A)||A| )
!! is computed by computing scaling factors R such that
!! diag(R)*A*op2(C) is row equilibrated and computing the standard
!! infinity-norm condition number.
iwork )
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_${rk}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: uplo
integer(${ik}$), intent(in) :: n, lda, ldaf, cmode
integer(${ik}$), intent(out) :: info
! Array Arguments
integer(${ik}$), intent(out) :: iwork(*)
integer(${ik}$), intent(in) :: ipiv(*)
real(${rk}$), intent(in) :: a(lda,*), af(ldaf,*), c(*)
real(${rk}$), intent(out) :: work(*)
! =====================================================================
! Local Scalars
character :: normin
integer(${ik}$) :: kase, i, j
real(${rk}$) :: ainvnm, smlnum, tmp
logical(lk) :: up
! Local Arrays
integer(${ik}$) :: isave(3_${ik}$)
! Intrinsic Functions
! Executable Statements
stdlib${ii}$_${ri}$la_syrcond = zero
info = 0_${ik}$
if( n<0_${ik}$ ) then
info = -2_${ik}$
else if( lda<max( 1_${ik}$, n ) ) then
info = -4_${ik}$
else if( ldaf<max( 1_${ik}$, n ) ) then
info = -6_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'DLA_SYRCOND', -info )
return
end if
if( n==0_${ik}$ ) then
stdlib${ii}$_${ri}$la_syrcond = one
return
end if
up = .false.
if ( stdlib_lsame( uplo, 'U' ) ) up = .true.
! compute the equilibration matrix r such that
! inv(r)*a*c has unit 1-norm.
if ( up ) then
do i = 1, n
tmp = zero
if ( cmode == 1_${ik}$ ) then
do j = 1, i
tmp = tmp + abs( a( j, i ) * c( j ) )
end do
do j = i+1, n
tmp = tmp + abs( a( i, j ) * c( j ) )
end do
else if ( cmode == 0_${ik}$ ) then
do j = 1, i
tmp = tmp + abs( a( j, i ) )
end do
do j = i+1, n
tmp = tmp + abs( a( i, j ) )
end do
else
do j = 1, i
tmp = tmp + abs( a( j, i ) / c( j ) )
end do
do j = i+1, n
tmp = tmp + abs( a( i, j ) / c( j ) )
end do
end if
work( 2_${ik}$*n+i ) = tmp
end do
else
do i = 1, n
tmp = zero
if ( cmode == 1_${ik}$ ) then
do j = 1, i
tmp = tmp + abs( a( i, j ) * c( j ) )
end do
do j = i+1, n
tmp = tmp + abs( a( j, i ) * c( j ) )
end do
else if ( cmode == 0_${ik}$ ) then
do j = 1, i
tmp = tmp + abs( a( i, j ) )
end do
do j = i+1, n
tmp = tmp + abs( a( j, i ) )
end do
else
do j = 1, i
tmp = tmp + abs( a( i, j) / c( j ) )
end do
do j = i+1, n
tmp = tmp + abs( a( j, i) / c( j ) )
end do
end if
work( 2_${ik}$*n+i ) = tmp
end do
endif
! estimate the norm of inv(op(a)).
smlnum = stdlib${ii}$_${ri}$lamch( 'SAFE MINIMUM' )
ainvnm = zero
normin = 'N'
kase = 0_${ik}$
10 continue
call stdlib${ii}$_${ri}$lacn2( n, work( n+1 ), work, iwork, ainvnm, kase, isave )
if( kase/=0_${ik}$ ) then
if( kase==2_${ik}$ ) then
! multiply by r.
do i = 1, n
work( i ) = work( i ) * work( 2_${ik}$*n+i )
end do
if ( up ) then
call stdlib${ii}$_${ri}$sytrs( 'U', n, 1_${ik}$, af, ldaf, ipiv, work, n, info )
else
call stdlib${ii}$_${ri}$sytrs( 'L', n, 1_${ik}$, af, ldaf, ipiv, work, n, info )
endif
! multiply by inv(c).
if ( cmode == 1_${ik}$ ) then
do i = 1, n
work( i ) = work( i ) / c( i )
end do
else if ( cmode == -1_${ik}$ ) then
do i = 1, n
work( i ) = work( i ) * c( i )
end do
end if
else
! multiply by inv(c**t).
if ( cmode == 1_${ik}$ ) then
do i = 1, n
work( i ) = work( i ) / c( i )
end do
else if ( cmode == -1_${ik}$ ) then
do i = 1, n
work( i ) = work( i ) * c( i )
end do
end if
if ( up ) then
call stdlib${ii}$_${ri}$sytrs( 'U', n, 1_${ik}$, af, ldaf, ipiv, work, n, info )
else
call stdlib${ii}$_${ri}$sytrs( 'L', n, 1_${ik}$, af, ldaf, ipiv, work, n, info )
endif
! multiply by r.
do i = 1, n
work( i ) = work( i ) * work( 2_${ik}$*n+i )
end do
end if
go to 10
end if
! compute the estimate of the reciprocal condition number.
if( ainvnm /= zero )stdlib${ii}$_${ri}$la_syrcond = ( one / ainvnm )
return
end function stdlib${ii}$_${ri}$la_syrcond
#:endif
#:endfor
real(sp) module function stdlib${ii}$_cla_syrcond_c( uplo, n, a, lda, af, ldaf, ipiv, c,capply, info, &
!! CLA_SYRCOND_C Computes the infinity norm condition number of
!! op(A) * inv(diag(C)) where C is a REAL vector.
work, rwork )
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: uplo
logical(lk), intent(in) :: capply
integer(${ik}$), intent(in) :: n, lda, ldaf
integer(${ik}$), intent(out) :: info
! Array Arguments
integer(${ik}$), intent(in) :: ipiv(*)
complex(sp), intent(in) :: a(lda,*), af(ldaf,*)
complex(sp), intent(out) :: work(*)
real(sp), intent(in) :: c(*)
real(sp), intent(out) :: rwork(*)
! =====================================================================
! Local Scalars
integer(${ik}$) :: kase
real(sp) :: ainvnm, anorm, tmp
integer(${ik}$) :: i, j
logical(lk) :: up, upper
complex(sp) :: zdum
! Local Arrays
integer(${ik}$) :: isave(3_${ik}$)
! Intrinsic Functions
! Statement Functions
real(sp) :: cabs1
! Statement Function Definitions
cabs1( zdum ) = abs( real( zdum,KIND=sp) ) + abs( aimag( zdum ) )
! Executable Statements
stdlib${ii}$_cla_syrcond_c = zero
info = 0_${ik}$
upper = stdlib_lsame( uplo, 'U' )
if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then
info = -1_${ik}$
else if( n<0_${ik}$ ) then
info = -2_${ik}$
else if( lda<max( 1_${ik}$, n ) ) then
info = -4_${ik}$
else if( ldaf<max( 1_${ik}$, n ) ) then
info = -6_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'CLA_SYRCOND_C', -info )
return
end if
up = .false.
if ( stdlib_lsame( uplo, 'U' ) ) up = .true.
! compute norm of op(a)*op2(c).
anorm = zero
if ( up ) then
do i = 1, n
tmp = zero
if ( capply ) then
do j = 1, i
tmp = tmp + cabs1( a( j, i ) ) / c( j )
end do
do j = i+1, n
tmp = tmp + cabs1( a( i, j ) ) / c( j )
end do
else
do j = 1, i
tmp = tmp + cabs1( a( j, i ) )
end do
do j = i+1, n
tmp = tmp + cabs1( a( i, j ) )
end do
end if
rwork( i ) = tmp
anorm = max( anorm, tmp )
end do
else
do i = 1, n
tmp = zero
if ( capply ) then
do j = 1, i
tmp = tmp + cabs1( a( i, j ) ) / c( j )
end do
do j = i+1, n
tmp = tmp + cabs1( a( j, i ) ) / c( j )
end do
else
do j = 1, i
tmp = tmp + cabs1( a( i, j ) )
end do
do j = i+1, n
tmp = tmp + cabs1( a( j, i ) )
end do
end if
rwork( i ) = tmp
anorm = max( anorm, tmp )
end do
end if
! quick return if possible.
if( n==0_${ik}$ ) then
stdlib${ii}$_cla_syrcond_c = one
return
else if( anorm == zero ) then
return
end if
! estimate the norm of inv(op(a)).
ainvnm = zero
kase = 0_${ik}$
10 continue
call stdlib${ii}$_clacn2( n, work( n+1 ), work, ainvnm, kase, isave )
if( kase/=0_${ik}$ ) then
if( kase==2_${ik}$ ) then
! multiply by r.
do i = 1, n
work( i ) = work( i ) * rwork( i )
end do
if ( up ) then
call stdlib${ii}$_csytrs( 'U', n, 1_${ik}$, af, ldaf, ipiv,work, n, info )
else
call stdlib${ii}$_csytrs( 'L', n, 1_${ik}$, af, ldaf, ipiv,work, n, info )
endif
! multiply by inv(c).
if ( capply ) then
do i = 1, n
work( i ) = work( i ) * c( i )
end do
end if
else
! multiply by inv(c**t).
if ( capply ) then
do i = 1, n
work( i ) = work( i ) * c( i )
end do
end if
if ( up ) then
call stdlib${ii}$_csytrs( 'U', n, 1_${ik}$, af, ldaf, ipiv,work, n, info )
else
call stdlib${ii}$_csytrs( 'L', n, 1_${ik}$, af, ldaf, ipiv,work, n, info )
end if
! multiply by r.
do i = 1, n
work( i ) = work( i ) * rwork( i )
end do
end if
go to 10
end if
! compute the estimate of the reciprocal condition number.
if( ainvnm /= zero )stdlib${ii}$_cla_syrcond_c = one / ainvnm
return
end function stdlib${ii}$_cla_syrcond_c
real(dp) module function stdlib${ii}$_zla_syrcond_c( uplo, n, a, lda, af,ldaf, ipiv, c, capply,info, work,&
!! ZLA_SYRCOND_C Computes the infinity norm condition number of
!! op(A) * inv(diag(C)) where C is a DOUBLE PRECISION vector.
rwork )
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: uplo
logical(lk), intent(in) :: capply
integer(${ik}$), intent(in) :: n, lda, ldaf
integer(${ik}$), intent(out) :: info
! Array Arguments
integer(${ik}$), intent(in) :: ipiv(*)
complex(dp), intent(in) :: a(lda,*), af(ldaf,*)
complex(dp), intent(out) :: work(*)
real(dp), intent(in) :: c(*)
real(dp), intent(out) :: rwork(*)
! =====================================================================
! Local Scalars
integer(${ik}$) :: kase
real(dp) :: ainvnm, anorm, tmp
integer(${ik}$) :: i, j
logical(lk) :: up, upper
complex(dp) :: zdum
! Local Arrays
integer(${ik}$) :: isave(3_${ik}$)
! Intrinsic Functions
! Statement Functions
real(dp) :: cabs1
! Statement Function Definitions
cabs1( zdum ) = abs( real( zdum,KIND=dp) ) + abs( aimag( zdum ) )
! Executable Statements
stdlib${ii}$_zla_syrcond_c = zero
info = 0_${ik}$
upper = stdlib_lsame( uplo, 'U' )
if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then
info = -1_${ik}$
else if( n<0_${ik}$ ) then
info = -2_${ik}$
else if( lda<max( 1_${ik}$, n ) ) then
info = -4_${ik}$
else if( ldaf<max( 1_${ik}$, n ) ) then
info = -6_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'ZLA_SYRCOND_C', -info )
return
end if
up = .false.
if ( stdlib_lsame( uplo, 'U' ) ) up = .true.
! compute norm of op(a)*op2(c).
anorm = zero
if ( up ) then
do i = 1, n
tmp = zero
if ( capply ) then
do j = 1, i
tmp = tmp + cabs1( a( j, i ) ) / c( j )
end do
do j = i+1, n
tmp = tmp + cabs1( a( i, j ) ) / c( j )
end do
else
do j = 1, i
tmp = tmp + cabs1( a( j, i ) )
end do
do j = i+1, n
tmp = tmp + cabs1( a( i, j ) )
end do
end if
rwork( i ) = tmp
anorm = max( anorm, tmp )
end do
else
do i = 1, n
tmp = zero
if ( capply ) then
do j = 1, i
tmp = tmp + cabs1( a( i, j ) ) / c( j )
end do
do j = i+1, n
tmp = tmp + cabs1( a( j, i ) ) / c( j )
end do
else
do j = 1, i
tmp = tmp + cabs1( a( i, j ) )
end do
do j = i+1, n
tmp = tmp + cabs1( a( j, i ) )
end do
end if
rwork( i ) = tmp
anorm = max( anorm, tmp )
end do
end if
! quick return if possible.
if( n==0_${ik}$ ) then
stdlib${ii}$_zla_syrcond_c = one
return
else if( anorm == zero ) then
return
end if
! estimate the norm of inv(op(a)).
ainvnm = zero
kase = 0_${ik}$
10 continue
call stdlib${ii}$_zlacn2( n, work( n+1 ), work, ainvnm, kase, isave )
if( kase/=0_${ik}$ ) then
if( kase==2_${ik}$ ) then
! multiply by r.
do i = 1, n
work( i ) = work( i ) * rwork( i )
end do
if ( up ) then
call stdlib${ii}$_zsytrs( 'U', n, 1_${ik}$, af, ldaf, ipiv,work, n, info )
else
call stdlib${ii}$_zsytrs( 'L', n, 1_${ik}$, af, ldaf, ipiv,work, n, info )
endif
! multiply by inv(c).
if ( capply ) then
do i = 1, n
work( i ) = work( i ) * c( i )
end do
end if
else
! multiply by inv(c**t).
if ( capply ) then
do i = 1, n
work( i ) = work( i ) * c( i )
end do
end if
if ( up ) then
call stdlib${ii}$_zsytrs( 'U', n, 1_${ik}$, af, ldaf, ipiv,work, n, info )
else
call stdlib${ii}$_zsytrs( 'L', n, 1_${ik}$, af, ldaf, ipiv,work, n, info )
end if
! multiply by r.
do i = 1, n
work( i ) = work( i ) * rwork( i )
end do
end if
go to 10
end if
! compute the estimate of the reciprocal condition number.
if( ainvnm /= zero )stdlib${ii}$_zla_syrcond_c = one / ainvnm
return
end function stdlib${ii}$_zla_syrcond_c
#:for ck,ct,ci in CMPLX_KINDS_TYPES
#:if not ck in ["sp","dp"]
real(${ck}$) module function stdlib${ii}$_${ci}$la_syrcond_c( uplo, n, a, lda, af,ldaf, ipiv, c, capply,info, work,&
!! ZLA_SYRCOND_C: Computes the infinity norm condition number of
!! op(A) * inv(diag(C)) where C is a DOUBLE PRECISION vector.
rwork )
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_${ck}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: uplo
logical(lk), intent(in) :: capply
integer(${ik}$), intent(in) :: n, lda, ldaf
integer(${ik}$), intent(out) :: info
! Array Arguments
integer(${ik}$), intent(in) :: ipiv(*)
complex(${ck}$), intent(in) :: a(lda,*), af(ldaf,*)
complex(${ck}$), intent(out) :: work(*)
real(${ck}$), intent(in) :: c(*)
real(${ck}$), intent(out) :: rwork(*)
! =====================================================================
! Local Scalars
integer(${ik}$) :: kase
real(${ck}$) :: ainvnm, anorm, tmp
integer(${ik}$) :: i, j
logical(lk) :: up, upper
complex(${ck}$) :: zdum
! Local Arrays
integer(${ik}$) :: isave(3_${ik}$)
! Intrinsic Functions
! Statement Functions
real(${ck}$) :: cabs1
! Statement Function Definitions
cabs1( zdum ) = abs( real( zdum,KIND=${ck}$) ) + abs( aimag( zdum ) )
! Executable Statements
stdlib${ii}$_${ci}$la_syrcond_c = zero
info = 0_${ik}$
upper = stdlib_lsame( uplo, 'U' )
if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then
info = -1_${ik}$
else if( n<0_${ik}$ ) then
info = -2_${ik}$
else if( lda<max( 1_${ik}$, n ) ) then
info = -4_${ik}$
else if( ldaf<max( 1_${ik}$, n ) ) then
info = -6_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'ZLA_SYRCOND_C', -info )
return
end if
up = .false.
if ( stdlib_lsame( uplo, 'U' ) ) up = .true.
! compute norm of op(a)*op2(c).
anorm = zero
if ( up ) then
do i = 1, n
tmp = zero
if ( capply ) then
do j = 1, i
tmp = tmp + cabs1( a( j, i ) ) / c( j )
end do
do j = i+1, n
tmp = tmp + cabs1( a( i, j ) ) / c( j )
end do
else
do j = 1, i
tmp = tmp + cabs1( a( j, i ) )
end do
do j = i+1, n
tmp = tmp + cabs1( a( i, j ) )
end do
end if
rwork( i ) = tmp
anorm = max( anorm, tmp )
end do
else
do i = 1, n
tmp = zero
if ( capply ) then
do j = 1, i
tmp = tmp + cabs1( a( i, j ) ) / c( j )
end do
do j = i+1, n
tmp = tmp + cabs1( a( j, i ) ) / c( j )
end do
else
do j = 1, i
tmp = tmp + cabs1( a( i, j ) )
end do
do j = i+1, n
tmp = tmp + cabs1( a( j, i ) )
end do
end if
rwork( i ) = tmp
anorm = max( anorm, tmp )
end do
end if
! quick return if possible.
if( n==0_${ik}$ ) then
stdlib${ii}$_${ci}$la_syrcond_c = one
return
else if( anorm == zero ) then
return
end if
! estimate the norm of inv(op(a)).
ainvnm = zero
kase = 0_${ik}$
10 continue
call stdlib${ii}$_${ci}$lacn2( n, work( n+1 ), work, ainvnm, kase, isave )
if( kase/=0_${ik}$ ) then
if( kase==2_${ik}$ ) then
! multiply by r.
do i = 1, n
work( i ) = work( i ) * rwork( i )
end do
if ( up ) then
call stdlib${ii}$_${ci}$sytrs( 'U', n, 1_${ik}$, af, ldaf, ipiv,work, n, info )
else
call stdlib${ii}$_${ci}$sytrs( 'L', n, 1_${ik}$, af, ldaf, ipiv,work, n, info )
endif
! multiply by inv(c).
if ( capply ) then
do i = 1, n
work( i ) = work( i ) * c( i )
end do
end if
else
! multiply by inv(c**t).
if ( capply ) then
do i = 1, n
work( i ) = work( i ) * c( i )
end do
end if
if ( up ) then
call stdlib${ii}$_${ci}$sytrs( 'U', n, 1_${ik}$, af, ldaf, ipiv,work, n, info )
else
call stdlib${ii}$_${ci}$sytrs( 'L', n, 1_${ik}$, af, ldaf, ipiv,work, n, info )
end if
! multiply by r.
do i = 1, n
work( i ) = work( i ) * rwork( i )
end do
end if
go to 10
end if
! compute the estimate of the reciprocal condition number.
if( ainvnm /= zero )stdlib${ii}$_${ci}$la_syrcond_c = one / ainvnm
return
end function stdlib${ii}$_${ci}$la_syrcond_c
#:endif
#:endfor
real(sp) module function stdlib${ii}$_cla_porcond_c( uplo, n, a, lda, af, ldaf, c, capply,info, work, &
!! CLA_PORCOND_C Computes the infinity norm condition number of
!! op(A) * inv(diag(C)) where C is a REAL vector
rwork )
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_sp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: uplo
logical(lk), intent(in) :: capply
integer(${ik}$), intent(in) :: n, lda, ldaf
integer(${ik}$), intent(out) :: info
! Array Arguments
complex(sp), intent(in) :: a(lda,*), af(ldaf,*)
complex(sp), intent(out) :: work(*)
real(sp), intent(in) :: c(*)
real(sp), intent(out) :: rwork(*)
! =====================================================================
! Local Scalars
integer(${ik}$) :: kase
real(sp) :: ainvnm, anorm, tmp
integer(${ik}$) :: i, j
logical(lk) :: up, upper
complex(sp) :: zdum
! Local Arrays
integer(${ik}$) :: isave(3_${ik}$)
! Intrinsic Functions
! Statement Functions
real(sp) :: cabs1
! Statement Function Definitions
cabs1( zdum ) = abs( real( zdum,KIND=sp) ) + abs( aimag( zdum ) )
! Executable Statements
stdlib${ii}$_cla_porcond_c = zero
info = 0_${ik}$
upper = stdlib_lsame( uplo, 'U' )
if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then
info = -1_${ik}$
else if( n<0_${ik}$ ) then
info = -2_${ik}$
else if( lda<max( 1_${ik}$, n ) ) then
info = -4_${ik}$
else if( ldaf<max( 1_${ik}$, n ) ) then
info = -6_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'CLA_PORCOND_C', -info )
return
end if
up = .false.
if ( stdlib_lsame( uplo, 'U' ) ) up = .true.
! compute norm of op(a)*op2(c).
anorm = zero
if ( up ) then
do i = 1, n
tmp = zero
if ( capply ) then
do j = 1, i
tmp = tmp + cabs1( a( j, i ) ) / c( j )
end do
do j = i+1, n
tmp = tmp + cabs1( a( i, j ) ) / c( j )
end do
else
do j = 1, i
tmp = tmp + cabs1( a( j, i ) )
end do
do j = i+1, n
tmp = tmp + cabs1( a( i, j ) )
end do
end if
rwork( i ) = tmp
anorm = max( anorm, tmp )
end do
else
do i = 1, n
tmp = zero
if ( capply ) then
do j = 1, i
tmp = tmp + cabs1( a( i, j ) ) / c( j )
end do
do j = i+1, n
tmp = tmp + cabs1( a( j, i ) ) / c( j )
end do
else
do j = 1, i
tmp = tmp + cabs1( a( i, j ) )
end do
do j = i+1, n
tmp = tmp + cabs1( a( j, i ) )
end do
end if
rwork( i ) = tmp
anorm = max( anorm, tmp )
end do
end if
! quick return if possible.
if( n==0_${ik}$ ) then
stdlib${ii}$_cla_porcond_c = one
return
else if( anorm == zero ) then
return
end if
! estimate the norm of inv(op(a)).
ainvnm = zero
kase = 0_${ik}$
10 continue
call stdlib${ii}$_clacn2( n, work( n+1 ), work, ainvnm, kase, isave )
if( kase/=0_${ik}$ ) then
if( kase==2_${ik}$ ) then
! multiply by r.
do i = 1, n
work( i ) = work( i ) * rwork( i )
end do
if ( up ) then
call stdlib${ii}$_cpotrs( 'U', n, 1_${ik}$, af, ldaf,work, n, info )
else
call stdlib${ii}$_cpotrs( 'L', n, 1_${ik}$, af, ldaf,work, n, info )
endif
! multiply by inv(c).
if ( capply ) then
do i = 1, n
work( i ) = work( i ) * c( i )
end do
end if
else
! multiply by inv(c**h).
if ( capply ) then
do i = 1, n
work( i ) = work( i ) * c( i )
end do
end if
if ( up ) then
call stdlib${ii}$_cpotrs( 'U', n, 1_${ik}$, af, ldaf,work, n, info )
else
call stdlib${ii}$_cpotrs( 'L', n, 1_${ik}$, af, ldaf,work, n, info )
end if
! multiply by r.
do i = 1, n
work( i ) = work( i ) * rwork( i )
end do
end if
go to 10
end if
! compute the estimate of the reciprocal condition number.
if( ainvnm /= zero )stdlib${ii}$_cla_porcond_c = one / ainvnm
return
end function stdlib${ii}$_cla_porcond_c
real(dp) module function stdlib${ii}$_zla_porcond_c( uplo, n, a, lda, af,ldaf, c, capply, info,work, &
!! ZLA_PORCOND_C Computes the infinity norm condition number of
!! op(A) * inv(diag(C)) where C is a DOUBLE PRECISION vector
rwork )
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_dp, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: uplo
logical(lk), intent(in) :: capply
integer(${ik}$), intent(in) :: n, lda, ldaf
integer(${ik}$), intent(out) :: info
! Array Arguments
complex(dp), intent(in) :: a(lda,*), af(ldaf,*)
complex(dp), intent(out) :: work(*)
real(dp), intent(in) :: c(*)
real(dp), intent(out) :: rwork(*)
! =====================================================================
! Local Scalars
integer(${ik}$) :: kase
real(dp) :: ainvnm, anorm, tmp
integer(${ik}$) :: i, j
logical(lk) :: up, upper
complex(dp) :: zdum
! Local Arrays
integer(${ik}$) :: isave(3_${ik}$)
! Intrinsic Functions
! Statement Functions
real(dp) :: cabs1
! Statement Function Definitions
cabs1( zdum ) = abs( real( zdum,KIND=dp) ) + abs( aimag( zdum ) )
! Executable Statements
stdlib${ii}$_zla_porcond_c = zero
info = 0_${ik}$
upper = stdlib_lsame( uplo, 'U' )
if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then
info = -1_${ik}$
else if( n<0_${ik}$ ) then
info = -2_${ik}$
else if( lda<max( 1_${ik}$, n ) ) then
info = -4_${ik}$
else if( ldaf<max( 1_${ik}$, n ) ) then
info = -6_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'ZLA_PORCOND_C', -info )
return
end if
up = .false.
if ( stdlib_lsame( uplo, 'U' ) ) up = .true.
! compute norm of op(a)*op2(c).
anorm = zero
if ( up ) then
do i = 1, n
tmp = zero
if ( capply ) then
do j = 1, i
tmp = tmp + cabs1( a( j, i ) ) / c( j )
end do
do j = i+1, n
tmp = tmp + cabs1( a( i, j ) ) / c( j )
end do
else
do j = 1, i
tmp = tmp + cabs1( a( j, i ) )
end do
do j = i+1, n
tmp = tmp + cabs1( a( i, j ) )
end do
end if
rwork( i ) = tmp
anorm = max( anorm, tmp )
end do
else
do i = 1, n
tmp = zero
if ( capply ) then
do j = 1, i
tmp = tmp + cabs1( a( i, j ) ) / c( j )
end do
do j = i+1, n
tmp = tmp + cabs1( a( j, i ) ) / c( j )
end do
else
do j = 1, i
tmp = tmp + cabs1( a( i, j ) )
end do
do j = i+1, n
tmp = tmp + cabs1( a( j, i ) )
end do
end if
rwork( i ) = tmp
anorm = max( anorm, tmp )
end do
end if
! quick return if possible.
if( n==0_${ik}$ ) then
stdlib${ii}$_zla_porcond_c = one
return
else if( anorm == zero ) then
return
end if
! estimate the norm of inv(op(a)).
ainvnm = zero
kase = 0_${ik}$
10 continue
call stdlib${ii}$_zlacn2( n, work( n+1 ), work, ainvnm, kase, isave )
if( kase/=0_${ik}$ ) then
if( kase==2_${ik}$ ) then
! multiply by r.
do i = 1, n
work( i ) = work( i ) * rwork( i )
end do
if ( up ) then
call stdlib${ii}$_zpotrs( 'U', n, 1_${ik}$, af, ldaf,work, n, info )
else
call stdlib${ii}$_zpotrs( 'L', n, 1_${ik}$, af, ldaf,work, n, info )
endif
! multiply by inv(c).
if ( capply ) then
do i = 1, n
work( i ) = work( i ) * c( i )
end do
end if
else
! multiply by inv(c**h).
if ( capply ) then
do i = 1, n
work( i ) = work( i ) * c( i )
end do
end if
if ( up ) then
call stdlib${ii}$_zpotrs( 'U', n, 1_${ik}$, af, ldaf,work, n, info )
else
call stdlib${ii}$_zpotrs( 'L', n, 1_${ik}$, af, ldaf,work, n, info )
end if
! multiply by r.
do i = 1, n
work( i ) = work( i ) * rwork( i )
end do
end if
go to 10
end if
! compute the estimate of the reciprocal condition number.
if( ainvnm /= zero )stdlib${ii}$_zla_porcond_c = one / ainvnm
return
end function stdlib${ii}$_zla_porcond_c
#:for ck,ct,ci in CMPLX_KINDS_TYPES
#:if not ck in ["sp","dp"]
real(${ck}$) module function stdlib${ii}$_${ci}$la_porcond_c( uplo, n, a, lda, af,ldaf, c, capply, info,work, &
!! ZLA_PORCOND_C: Computes the infinity norm condition number of
!! op(A) * inv(diag(C)) where C is a DOUBLE PRECISION vector
rwork )
! -- lapack computational routine --
! -- lapack is a software package provided by univ. of tennessee, --
! -- univ. of california berkeley, univ. of colorado denver and nag ltd..--
use stdlib_blas_constants_${ck}$, only: negone, zero, half, one, two, three, four, eight, ten, czero, chalf, cone, cnegone
! Scalar Arguments
character, intent(in) :: uplo
logical(lk), intent(in) :: capply
integer(${ik}$), intent(in) :: n, lda, ldaf
integer(${ik}$), intent(out) :: info
! Array Arguments
complex(${ck}$), intent(in) :: a(lda,*), af(ldaf,*)
complex(${ck}$), intent(out) :: work(*)
real(${ck}$), intent(in) :: c(*)
real(${ck}$), intent(out) :: rwork(*)
! =====================================================================
! Local Scalars
integer(${ik}$) :: kase
real(${ck}$) :: ainvnm, anorm, tmp
integer(${ik}$) :: i, j
logical(lk) :: up, upper
complex(${ck}$) :: zdum
! Local Arrays
integer(${ik}$) :: isave(3_${ik}$)
! Intrinsic Functions
! Statement Functions
real(${ck}$) :: cabs1
! Statement Function Definitions
cabs1( zdum ) = abs( real( zdum,KIND=${ck}$) ) + abs( aimag( zdum ) )
! Executable Statements
stdlib${ii}$_${ci}$la_porcond_c = zero
info = 0_${ik}$
upper = stdlib_lsame( uplo, 'U' )
if( .not.upper .and. .not.stdlib_lsame( uplo, 'L' ) ) then
info = -1_${ik}$
else if( n<0_${ik}$ ) then
info = -2_${ik}$
else if( lda<max( 1_${ik}$, n ) ) then
info = -4_${ik}$
else if( ldaf<max( 1_${ik}$, n ) ) then
info = -6_${ik}$
end if
if( info/=0_${ik}$ ) then
call stdlib${ii}$_xerbla( 'ZLA_PORCOND_C', -info )
return
end if
up = .false.
if ( stdlib_lsame( uplo, 'U' ) ) up = .true.
! compute norm of op(a)*op2(c).
anorm = zero
if ( up ) then
do i = 1, n
tmp = zero
if ( capply ) then
do j = 1, i
tmp = tmp + cabs1( a( j, i ) ) / c( j )
end do
do j = i+1, n
tmp = tmp + cabs1( a( i, j ) ) / c( j )
end do
else
do j = 1, i
tmp = tmp + cabs1( a( j, i ) )
end do
do j = i+1, n
tmp = tmp + cabs1( a( i, j ) )
end do
end if
rwork( i ) = tmp
anorm = max( anorm, tmp )
end do
else
do i = 1, n
tmp = zero
if ( capply ) then
do j = 1, i
tmp = tmp + cabs1( a( i, j ) ) / c( j )
end do
do j = i+1, n
tmp = tmp + cabs1( a( j, i ) ) / c( j )
end do
else
do j = 1, i
tmp = tmp + cabs1( a( i, j ) )
end do
do j = i+1, n
tmp = tmp + cabs1( a( j, i ) )
end do
end if
rwork( i ) = tmp
anorm = max( anorm, tmp )
end do
end if
! quick return if possible.
if( n==0_${ik}$ ) then
stdlib${ii}$_${ci}$la_porcond_c = one
return
else if( anorm == zero ) then
return
end if
! estimate the norm of inv(op(a)).
ainvnm = zero
kase = 0_${ik}$
10 continue
call stdlib${ii}$_${ci}$lacn2( n, work( n+1 ), work, ainvnm, kase, isave )
if( kase/=0_${ik}$ ) then
if( kase==2_${ik}$ ) then
! multiply by r.
do i = 1, n
work( i ) = work( i ) * rwork( i )
end do
if ( up ) then
call stdlib${ii}$_${ci}$potrs( 'U', n, 1_${ik}$, af, ldaf,work, n, info )
else
call stdlib${ii}$_${ci}$potrs( 'L', n, 1_${ik}$, af, ldaf,work, n, info )
endif
! multiply by inv(c).
if ( capply ) then
do i = 1, n
work( i ) = work( i ) * c( i )
end do
end if
else
! multiply by inv(c**h).
if ( capply ) then
do i = 1, n
work( i ) = work( i ) * c( i )
end do
end if
if ( up ) then
call stdlib${ii}$_${ci}$potrs( 'U', n, 1_${ik}$, af, ldaf,work, n, info )
else
call stdlib${ii}$_${ci}$potrs( 'L', n, 1_${ik}$, af, ldaf,work, n, info )
end if
! multiply by r.
do i = 1, n
work( i ) = work( i ) * rwork( i )
end do
end if
go to 10
end if
! compute the estimate of the reciprocal condition number.
if( ainvnm /= zero )stdlib${ii}$_${ci}$la_porcond_c = one / ainvnm
return
end function stdlib${ii}$_${ci}$la_porcond_c
#:endif
#:endfor
#:endfor
end submodule stdlib_lapack_others_sm
