219 SUBROUTINE dspt21( ITYPE, UPLO, N, KBAND, AP, D, E, U, LDU, VP,
220 $ tau, work, result )
229 INTEGER itype, kband, ldu, n
232 DOUBLE PRECISION ap( * ), d( * ), e( * ), result( 2 ), tau( * ),
233 $ u( ldu, * ), vp( * ), work( * )
239 DOUBLE PRECISION zero, one, ten
240 parameter( zero = 0.0d0, one = 1.0d0, ten = 10.0d0 )
241 DOUBLE PRECISION half
242 parameter( half = 1.0d+0 / 2.0d+0 )
247 INTEGER iinfo,
j, jp, jp1, jr, lap
248 DOUBLE PRECISION anorm, temp, ulp, unfl, vsave, wnorm
260 INTRINSIC dble, max, min
272 lap = ( n*( n+1 ) ) / 2
274 IF(
lsame( uplo,
'U' ) )
THEN
282 unfl =
dlamch(
'Safe minimum' )
287 IF( itype.LT.1 .OR. itype.GT.3 )
THEN
288 result( 1 ) = ten / ulp
296 IF( itype.EQ.3 )
THEN
299 anorm = max(
dlansp(
'1', cuplo, n, ap, work ), unfl )
304 IF( itype.EQ.1 )
THEN
308 CALL
dlaset(
'Full', n, n, zero, zero, work, n )
309 CALL
dcopy( lap, ap, 1, work, 1 )
312 CALL
dspr( cuplo, n, -d(
j ), u( 1,
j ), 1, work )
315 IF( n.GT.1 .AND. kband.EQ.1 )
THEN
317 CALL
dspr2( cuplo, n, -e(
j ), u( 1,
j ), 1, u( 1,
j+1 ),
321 wnorm =
dlansp(
'1', cuplo, n, work, work( n**2+1 ) )
323 ELSE IF( itype.EQ.2 )
THEN
327 CALL
dlaset(
'Full', n, n, zero, zero, work, n )
331 DO 40
j = n - 1, 1, -1
332 jp = ( ( 2*n-
j )*(
j-1 ) ) / 2
334 IF( kband.EQ.1 )
THEN
335 work( jp+
j+1 ) = ( one-tau(
j ) )*e(
j )
337 work( jp+jr ) = -tau(
j )*e(
j )*vp( jp+jr )
341 IF( tau(
j ).NE.zero )
THEN
344 CALL
dspmv(
'L', n-
j, one, work( jp1+
j+1 ),
345 $ vp( jp+
j+1 ), 1, zero, work( lap+1 ), 1 )
346 temp = -half*tau(
j )*
ddot( n-
j, work( lap+1 ), 1,
348 CALL
daxpy( n-
j, temp, vp( jp+
j+1 ), 1, work( lap+1 ),
350 CALL
dspr2(
'L', n-
j, -tau(
j ), vp( jp+
j+1 ), 1,
351 $ work( lap+1 ), 1, work( jp1+
j+1 ) )
354 work( jp+
j ) = d(
j )
359 jp = (
j*(
j-1 ) ) / 2
361 IF( kband.EQ.1 )
THEN
362 work( jp1+
j ) = ( one-tau(
j ) )*e(
j )
364 work( jp1+jr ) = -tau(
j )*e(
j )*vp( jp1+jr )
368 IF( tau(
j ).NE.zero )
THEN
371 CALL
dspmv(
'U',
j, one, work, vp( jp1+1 ), 1, zero,
373 temp = -half*tau(
j )*
ddot(
j, work( lap+1 ), 1,
375 CALL
daxpy(
j, temp, vp( jp1+1 ), 1, work( lap+1 ),
377 CALL
dspr2(
'U',
j, -tau(
j ), vp( jp1+1 ), 1,
378 $ work( lap+1 ), 1, work )
381 work( jp1+
j+1 ) = d(
j+1 )
386 work(
j ) = work(
j ) - ap(
j )
388 wnorm =
dlansp(
'1', cuplo, n, work, work( lap+1 ) )
390 ELSE IF( itype.EQ.3 )
THEN
396 CALL
dlacpy(
' ', n, n, u, ldu, work, n )
397 CALL
dopmtr(
'R', cuplo,
'T', n, n, vp, tau, work, n,
398 $ work( n**2+1 ), iinfo )
399 IF( iinfo.NE.0 )
THEN
400 result( 1 ) = ten / ulp
405 work( ( n+1 )*(
j-1 )+1 ) = work( ( n+1 )*(
j-1 )+1 ) - one
408 wnorm =
dlange(
'1', n, n, work, n, work( n**2+1 ) )
411 IF( anorm.GT.wnorm )
THEN
412 result( 1 ) = ( wnorm / anorm ) / ( n*ulp )
414 IF( anorm.LT.one )
THEN
415 result( 1 ) = ( min( wnorm, n*anorm ) / anorm ) / ( n*ulp )
417 result( 1 ) = min( wnorm / anorm, dble( n ) ) / ( n*ulp )
425 IF( itype.EQ.1 )
THEN
426 CALL
dgemm(
'N',
'C', n, n, n, one, u, ldu, u, ldu, zero, work,
430 work( ( n+1 )*(
j-1 )+1 ) = work( ( n+1 )*(
j-1 )+1 ) - one
433 result( 2 ) = min(
dlange(
'1', n, n, work, n,
434 $ work( n**2+1 ) ), dble( n ) ) / ( n*ulp )
LOGICAL function lsame(CA, CB)
LSAME
subroutine dspr(UPLO, N, ALPHA, X, INCX, AP)
DSPR
subroutine dspmv(UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY)
DSPMV
subroutine dopmtr(SIDE, UPLO, TRANS, M, N, AP, TAU, C, LDC, WORK, INFO)
DOPMTR
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
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 ...
subroutine dspt21(ITYPE, UPLO, N, KBAND, AP, D, E, U, LDU, VP, TAU, WORK, RESULT)
DSPT21
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
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.
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
DOUBLE PRECISION function dlamch(CMACH)
DLAMCH
DOUBLE PRECISION function ddot(N, DX, INCX, DY, INCY)
DDOT
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...