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
252 DOUBLE PRECISION DDOT, DLAMCH, DLANGE, DLANSP
253 EXTERNAL lsame, ddot, dlamch, dlange, dlansp
260 INTRINSIC dble, max, min
272 lap = ( n*( n+1 ) ) / 2
274 IF( lsame( uplo,
'U' ) )
THEN
282 unfl = dlamch(
'Safe minimum' )
283 ulp = dlamch(
'Epsilon' )*dlamch(
'Base' )
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 )
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
subroutine dspmv(UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY)
DSPMV
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
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...
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
subroutine dspt21(ITYPE, UPLO, N, KBAND, AP, D, E, U, LDU, VP, TAU, WORK, RESULT)
DSPT21
subroutine dspr(UPLO, N, ALPHA, X, INCX, AP)
DSPR
subroutine dspr2(UPLO, N, ALPHA, X, INCX, Y, INCY, AP)
DSPR2
subroutine dopmtr(SIDE, UPLO, TRANS, M, N, AP, TAU, C, LDC, WORK, INFO)
DOPMTR