146 SUBROUTINE dsbt21( UPLO, N, KA, KS, A, LDA, D, E, U, LDU, WORK,
156 INTEGER ka, ks, lda, ldu, n
159 DOUBLE PRECISION a( lda, * ), d( * ), e( * ), result( 2 ),
160 $ u( ldu, * ), work( * )
166 DOUBLE PRECISION zero, one
167 parameter( zero = 0.0d0, one = 1.0d0 )
172 INTEGER ika,
j, jc, jr, lw
173 DOUBLE PRECISION anorm, ulp, unfl, wnorm
184 INTRINSIC dble, max, min
195 ika = max( 0, min( n-1, ka ) )
196 lw = ( n*( n+1 ) ) / 2
198 IF(
lsame( uplo,
'U' ) )
THEN
206 unfl =
dlamch(
'Safe minimum' )
215 anorm = max(
dlansb(
'1', cuplo, n, ika, a, lda, work ), unfl )
224 DO 10 jr = 1, min( ika+1, n+1-jc )
226 work(
j ) = a( jr, jc )
228 DO 20 jr = ika + 2, n + 1 - jc
233 DO 30 jr = ika + 2, jc
237 DO 40 jr = min( ika, jc-1 ), 0, -1
239 work(
j ) = a( ika+1-jr, jc )
245 CALL
dspr( cuplo, n, -d(
j ), u( 1,
j ), 1, work )
248 IF( n.GT.1 .AND. ks.EQ.1 )
THEN
250 CALL
dspr2( cuplo, n, -e(
j ), u( 1,
j ), 1, u( 1,
j+1 ), 1,
254 wnorm =
dlansp(
'1', cuplo, n, work, work( lw+1 ) )
256 IF( anorm.GT.wnorm )
THEN
257 result( 1 ) = ( wnorm / anorm ) / ( n*ulp )
259 IF( anorm.LT.one )
THEN
260 result( 1 ) = ( min( wnorm, n*anorm ) / anorm ) / ( n*ulp )
262 result( 1 ) = min( wnorm / anorm, dble( n ) ) / ( n*ulp )
270 CALL
dgemm(
'N',
'C', n, n, n, one, u, ldu, u, ldu, zero, work,
274 work( ( n+1 )*(
j-1 )+1 ) = work( ( n+1 )*(
j-1 )+1 ) - one
277 result( 2 ) = min(
dlange(
'1', n, n, work, n, work( n**2+1 ) ),
278 $ dble( n ) ) / ( n*ulp )
LOGICAL function lsame(CA, CB)
LSAME
subroutine dspr(UPLO, N, ALPHA, X, INCX, AP)
DSPR
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
subroutine dspr2(UPLO, N, ALPHA, X, INCX, Y, INCY, AP)
DSPR2
DOUBLE PRECISION function dlange(NORM, M, N, A, LDA, WORK)
DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
DOUBLE PRECISION function dlansb(NORM, UPLO, N, K, AB, LDAB, WORK)
DLANSB returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a symmetric band matrix.
subroutine dsbt21(UPLO, N, KA, KS, A, LDA, D, E, U, LDU, WORK, RESULT)
DSBT21
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
DOUBLE PRECISION function dlansp(NORM, UPLO, N, AP, WORK)
DLANSP returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a symmetric matrix supplied in packed form.
DOUBLE PRECISION function dlamch(CMACH)
DLAMCH