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 sspr2(UPLO, N, ALPHA, X, INCX, Y, INCY, AP)
SSPR2
subroutine sspmv(UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY)
SSPMV
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
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 sspt21(ITYPE, UPLO, N, KBAND, AP, D, E, U, LDU, VP, TAU, WORK, RESULT)
SSPT21
subroutine saxpy(N, SA, SX, INCX, SY, INCY)
SAXPY
subroutine sspr(UPLO, N, ALPHA, X, INCX, AP)
SSPR
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY