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
236 REAL SLAMCH, SLANGE, SLANSY
237 EXTERNAL lsame, slamch, slange, slansy
244 INTRINSIC max, min, real
254 IF( lsame( uplo,
'U' ) )
THEN
262 unfl = slamch(
'Safe minimum' )
263 ulp = slamch(
'Epsilon' )*slamch(
'Base' )
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 )
subroutine sorm2r(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, INFO)
SORM2R multiplies a general matrix by the orthogonal matrix from a QR factorization determined by sge...
subroutine ssyr2(UPLO, N, ALPHA, X, INCX, Y, INCY, A, LDA)
SSYR2
subroutine sorm2l(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, INFO)
SORM2L multiplies a general matrix by the orthogonal matrix from a QL factorization determined by sge...
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
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 ssyt21(ITYPE, UPLO, N, KBAND, A, LDA, D, E, U, LDU, V, LDV, TAU, WORK, RESULT)
SSYT21
subroutine slarfy(UPLO, N, V, INCV, TAU, C, LDC, WORK)
SLARFY
subroutine ssyr(UPLO, N, ALPHA, X, INCX, A, LDA)
SSYR
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.