219 SUBROUTINE sspt21( ITYPE, UPLO, N, KBAND, AP, D, E, U, LDU, VP,
220 $ tau, work, result )
229 INTEGER ITYPE, KBAND, LDU, N
232 REAL AP( * ), D( * ), E( * ), RESULT( 2 ), TAU( * ),
233 $ u( ldu, * ), vp( * ), work( * )
240 parameter( zero = 0.0e0, one = 1.0e0, ten = 10.0e0 )
242 parameter( half = 1.0e+0 / 2.0e+0 )
247 INTEGER IINFO, J, JP, JP1, JR, LAP
248 REAL ANORM, TEMP, ULP, UNFL, VSAVE, WNORM
252 REAL SDOT, SLAMCH, SLANGE, SLANSP
253 EXTERNAL lsame, sdot, slamch, slange, slansp
260 INTRINSIC max, min, real
272 lap = ( n*( n+1 ) ) / 2
274 IF( lsame( uplo,
'U' ) )
THEN
282 unfl = slamch(
'Safe minimum' )
283 ulp = slamch(
'Epsilon' )*slamch(
'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( slansp(
'1', cuplo, n, ap, work ), unfl )
304 IF( itype.EQ.1 )
THEN
308 CALL
slaset(
'Full', n, n, zero, zero, work, n )
309 CALL
scopy( lap, ap, 1, work, 1 )
312 CALL
sspr( cuplo, n, -d( j ), u( 1, j ), 1, work )
315 IF( n.GT.1 .AND. kband.EQ.1 )
THEN
317 CALL
sspr2( cuplo, n, -e( j ), u( 1, j ), 1, u( 1, j+1 ),
321 wnorm = slansp(
'1', cuplo, n, work, work( n**2+1 ) )
323 ELSE IF( itype.EQ.2 )
THEN
327 CALL
slaset(
'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
sspmv(
'L', n-j, one, work( jp1+j+1 ),
345 $ vp( jp+j+1 ), 1, zero, work( lap+1 ), 1 )
346 temp = -half*tau( j )*sdot( n-j, work( lap+1 ), 1,
348 CALL
saxpy( n-j, temp, vp( jp+j+1 ), 1, work( lap+1 ),
350 CALL
sspr2(
'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
sspmv(
'U', j, one, work, vp( jp1+1 ), 1, zero,
373 temp = -half*tau( j )*sdot( j, work( lap+1 ), 1,
375 CALL
saxpy( j, temp, vp( jp1+1 ), 1, work( lap+1 ),
377 CALL
sspr2(
'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 = slansp(
'1', cuplo, n, work, work( lap+1 ) )
390 ELSE IF( itype.EQ.3 )
THEN
396 CALL
slacpy(
' ', n, n, u, ldu, work, n )
397 CALL
sopmtr(
'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 = slange(
'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,
REAL( N ) ) / ( N*ULP )
425 IF( itype.EQ.1 )
THEN
426 CALL
sgemm(
'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( slange(
'1', n, n, work, n,
434 $ work( n**2+1 ) ),
REAL( N ) ) / ( N*ULP )
subroutine sopmtr(SIDE, UPLO, TRANS, M, N, AP, TAU, C, LDC, WORK, INFO)
SOPMTR
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine sspmv(UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY)
SSPMV
subroutine sspt21(ITYPE, UPLO, N, KBAND, AP, D, E, U, LDU, VP, TAU, WORK, RESULT)
SSPT21
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
subroutine saxpy(N, SA, SX, INCX, SY, INCY)
SAXPY
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
subroutine sspr2(UPLO, N, ALPHA, X, INCX, Y, INCY, AP)
SSPR2
subroutine sspr(UPLO, N, ALPHA, X, INCX, AP)
SSPR