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
260 INTRINSIC max, min, real
272 lap = ( n*( n+1 ) ) / 2
274 IF(
lsame( uplo,
'U' ) )
THEN
282 unfl =
slamch(
'Safe minimum' )
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...
real function slansp(NORM, UPLO, N, AP, WORK)
SLANSP returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a symmetric matrix supplied in packed form.
subroutine sspmv(UPLO, N, ALPHA, AP, X, INCX, BETA, Y, INCY)
SSPMV
real function slange(NORM, M, N, A, LDA, WORK)
SLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
real function sdot(N, SX, INCX, SY, INCY)
SDOT
subroutine sspt21(ITYPE, UPLO, N, KBAND, AP, D, E, U, LDU, VP, TAU, WORK, RESULT)
SSPT21
subroutine scopy(N, SX, INCX, SY, INCY)
SCOPY
logical function lsame(CA, CB)
LSAME
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
real function slamch(CMACH)
SLAMCH
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
subroutine sspr2(UPLO, N, ALPHA, X, INCX, Y, INCY, AP)
SSPR2
subroutine sspr(UPLO, N, ALPHA, X, INCX, AP)
SSPR