205 SUBROUTINE ssyt21( ITYPE, UPLO, N, KBAND, A, LDA, D, E, U, LDU, V,
206 $ ldv, tau, work, result )
215 INTEGER itype, kband, lda, ldu, ldv, n
218 REAL a( lda, * ), d( * ), e( * ), result( 2 ),
219 $ tau( * ), u( ldu, * ), v( ldv, * ), work( * )
226 parameter( zero = 0.0e0, one = 1.0e0, ten = 10.0e0 )
231 INTEGER iinfo, j, jcol, jr, jrow
232 REAL anorm, ulp, unfl, vsave, wnorm
244 INTRINSIC max, min, real
254 IF(
lsame( uplo,
'U' ) )
THEN
262 unfl =
slamch(
'Safe minimum' )
267 IF( itype.LT.1 .OR. itype.GT.3 )
THEN
268 result( 1 ) = ten / ulp
276 IF( itype.EQ.3 )
THEN
279 anorm = max(
slansy(
'1', cuplo, n, a, lda, work ), unfl )
284 IF( itype.EQ.1 )
THEN
288 CALL
slaset(
'Full', n, n, zero, zero, work, n )
289 CALL
slacpy( cuplo, n, n, a, lda, work, n )
292 CALL
ssyr( cuplo, n, -d( j ), u( 1, j ), 1, work, n )
295 IF( n.GT.1 .AND. kband.EQ.1 )
THEN
297 CALL
ssyr2( cuplo, n, -e( j ), u( 1, j ), 1, u( 1, j+1 ),
301 wnorm =
slansy(
'1', cuplo, n, work, n, work( n**2+1 ) )
303 ELSE IF( itype.EQ.2 )
THEN
307 CALL
slaset(
'Full', n, n, zero, zero, work, n )
310 work( n**2 ) = d( n )
311 DO 40 j = n - 1, 1, -1
312 IF( kband.EQ.1 )
THEN
313 work( ( n+1 )*( j-1 )+2 ) = ( one-tau( j ) )*e( j )
315 work( ( j-1 )*n+jr ) = -tau( j )*e( j )*v( jr, j )
321 CALL
slarfy(
'L', n-j, v( j+1, j ), 1, tau( j ),
322 $ work( ( n+1 )*j+1 ), n, work( n**2+1 ) )
324 work( ( n+1 )*( j-1 )+1 ) = d( j )
329 IF( kband.EQ.1 )
THEN
330 work( ( n+1 )*j ) = ( one-tau( j ) )*e( j )
332 work( j*n+jr ) = -tau( j )*e( j )*v( jr, j+1 )
338 CALL
slarfy(
'U', j, v( 1, j+1 ), 1, tau( j ), work, n,
341 work( ( n+1 )*j+1 ) = d( j+1 )
348 work( jrow+n*( jcol-1 ) ) = work( jrow+n*( jcol-1 ) )
353 work( jrow+n*( jcol-1 ) ) = work( jrow+n*( jcol-1 ) )
358 wnorm =
slansy(
'1', cuplo, n, work, n, work( n**2+1 ) )
360 ELSE IF( itype.EQ.3 )
THEN
366 CALL
slacpy(
' ', n, n, u, ldu, work, n )
368 CALL
sorm2r(
'R',
'T', n, n-1, n-1, v( 2, 1 ), ldv, tau,
369 $ work( n+1 ), n, work( n**2+1 ), iinfo )
371 CALL
sorm2l(
'R',
'T', n, n-1, n-1, v( 1, 2 ), ldv, tau,
372 $ work, n, work( n**2+1 ), iinfo )
374 IF( iinfo.NE.0 )
THEN
375 result( 1 ) = ten / ulp
380 work( ( n+1 )*( j-1 )+1 ) = work( ( n+1 )*( j-1 )+1 ) - one
383 wnorm =
slange(
'1', n, n, work, n, work( n**2+1 ) )
386 IF( anorm.GT.wnorm )
THEN
387 result( 1 ) = ( wnorm / anorm ) / ( n*ulp )
389 IF( anorm.LT.one )
THEN
390 result( 1 ) = ( min( wnorm, n*anorm ) / anorm ) / ( n*ulp )
392 result( 1 ) = min( wnorm / anorm,
REAL( N ) ) / ( n*ulp )
400 IF( itype.EQ.1 )
THEN
401 CALL
sgemm(
'N',
'C', n, n, n, one, u, ldu, u, ldu, zero, work,
405 work( ( n+1 )*( j-1 )+1 ) = work( ( n+1 )*( j-1 )+1 ) - one
408 result( 2 ) = min(
slange(
'1', n, n, work, n,
409 $ work( n**2+1 ) ),
REAL( N ) ) / ( n*ulp )