155 SUBROUTINE ssyt22( ITYPE, UPLO, N, M, KBAND, A, LDA, D, E, U, LDU,
156 $ v, ldv, tau, work, result )
165 INTEGER itype, kband, lda, ldu, ldv, m, n
168 REAL a( lda, * ), d( * ), e( * ), result( 2 ),
169 $ tau( * ), u( ldu, * ), v( ldv, * ), work( * )
176 parameter( zero = 0.0e0, one = 1.0e0 )
179 INTEGER j, jj, jj1, jj2, nn, nnp1
180 REAL anorm, ulp, unfl, wnorm
190 INTRINSIC max, min, real
196 IF( n.LE.0 .OR. m.LE.0 )
199 unfl =
slamch(
'Safe minimum' )
200 ulp =
slamch(
'Precision' )
206 anorm = max(
slansy(
'1', uplo, n, a, lda, work ), unfl )
212 CALL
ssymm(
'L', uplo, n, m, one, a, lda, u, ldu, zero, work, n )
215 CALL
sgemm(
'T',
'N', m, m, n, one, u, ldu, work, n, zero,
218 jj = nn + (
j-1 )*n +
j
219 work( jj ) = work( jj ) - d(
j )
221 IF( kband.EQ.1 .AND. n.GT.1 )
THEN
223 jj1 = nn + (
j-1 )*n +
j - 1
224 jj2 = nn + (
j-2 )*n +
j
225 work( jj1 ) = work( jj1 ) - e(
j-1 )
226 work( jj2 ) = work( jj2 ) - e(
j-1 )
229 wnorm =
slansy(
'1', uplo, m, work( nnp1 ), n, work( 1 ) )
231 IF( anorm.GT.wnorm )
THEN
232 result( 1 ) = ( wnorm / anorm ) / ( m*ulp )
234 IF( anorm.LT.one )
THEN
235 result( 1 ) = ( min( wnorm, m*anorm ) / anorm ) / ( m*ulp )
237 result( 1 ) = min( wnorm / anorm,
REAL( M ) ) / ( m*ulp )
246 $ CALL
sort01(
'Columns', n, m, u, ldu, work, 2*n*n,
REAL function slansy(NORM, UPLO, N, A, LDA, WORK)
SLANSY returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a real symmetric matrix.
subroutine ssymm(SIDE, UPLO, M, N, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SSYMM
REAL function slamch(CMACH)
SLAMCH
subroutine sort01(ROWCOL, M, N, U, LDU, WORK, LWORK, RESID)
SORT01
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real j
subroutine ssyt22(ITYPE, UPLO, N, M, KBAND, A, LDA, D, E, U, LDU, V, LDV, TAU, WORK, RESULT)
SSYT22