91 REAL,
ALLOCATABLE :: af(:,:), q(:,:),
92 $ r(:,:), rwork(:), work( : ), t(:,:),
93 $ cf(:,:), df(:,:), a(:,:), c(:,:), d(:,:)
97 parameter( zero = 0.0, one = 1.0 )
100 INTEGER info, j, k, l, lwork
101 REAL anorm, eps, resid, cnorm, dnorm
116 DATA iseed / 1988, 1989, 1990, 1991 /
121 lwork = max(2,l)*max(2,l)*nb
125 ALLOCATE ( a(m,n), af(m,n), q(m,m), r(m,l), rwork(l),
126 $ work(lwork), t(nb,n), c(m,n), cf(m,n),
133 CALL slarnv( 2, iseed, m, a( 1, j ) )
135 CALL slacpy(
'Full', m, n, a, m, af, m )
139 CALL sgeqrt( m, n, nb, af, m, t, ldt, work, info )
143 CALL slaset(
'Full', m, m, zero, one, q, m )
144 CALL sgemqrt(
'R',
'N', m, m, k, nb, af, m, t, ldt, q, m,
149 CALL slaset(
'Full', m, n, zero, zero, r, m )
150 CALL slacpy(
'Upper', m, n, af, m, r, m )
154 CALL sgemm(
'T',
'N', m, n, m, -one, q, m, a, m, one, r, m )
155 anorm =
slange(
'1', m, n, a, m, rwork )
156 resid =
slange(
'1', m, n, r, m, rwork )
157 IF( anorm.GT.zero )
THEN
158 result( 1 ) = resid / (eps*max(1,m)*anorm)
165 CALL slaset(
'Full', m, m, zero, one, r, m )
166 CALL ssyrk(
'U',
'C', m, m, -one, q, m, one, r, m )
167 resid =
slansy(
'1',
'Upper', m, r, m, rwork )
168 result( 2 ) = resid / (eps*max(1,m))
173 CALL slarnv( 2, iseed, m, c( 1, j ) )
175 cnorm =
slange(
'1', m, n, c, m, rwork)
176 CALL slacpy(
'Full', m, n, c, m, cf, m )
180 CALL sgemqrt(
'L',
'N', m, n, k, nb, af, m, t, nb, cf, m,
185 CALL sgemm(
'N',
'N', m, n, m, -one, q, m, c, m, one, cf, m )
186 resid =
slange(
'1', m, n, cf, m, rwork )
187 IF( cnorm.GT.zero )
THEN
188 result( 3 ) = resid / (eps*max(1,m)*cnorm)
195 CALL slacpy(
'Full', m, n, c, m, cf, m )
199 CALL sgemqrt(
'L',
'T', m, n, k, nb, af, m, t, nb, cf, m,
204 CALL sgemm(
'T',
'N', m, n, m, -one, q, m, c, m, one, cf, m )
205 resid =
slange(
'1', m, n, cf, m, rwork )
206 IF( cnorm.GT.zero )
THEN
207 result( 4 ) = resid / (eps*max(1,m)*cnorm)
215 CALL slarnv( 2, iseed, n, d( 1, j ) )
217 dnorm =
slange(
'1', n, m, d, n, rwork)
218 CALL slacpy(
'Full', n, m, d, n, df, n )
222 CALL sgemqrt(
'R',
'N', n, m, k, nb, af, m, t, nb, df, n,
227 CALL sgemm(
'N',
'N', n, m, m, -one, d, n, q, m, one, df, n )
228 resid =
slange(
'1', n, m, df, n, rwork )
229 IF( cnorm.GT.zero )
THEN
230 result( 5 ) = resid / (eps*max(1,m)*dnorm)
237 CALL slacpy(
'Full', n, m, d, n, df, n )
241 CALL sgemqrt(
'R',
'T', n, m, k, nb, af, m, t, nb, df, n,
246 CALL sgemm(
'N',
'T', n, m, m, -one, d, n, q, m, one, df, n )
247 resid =
slange(
'1', n, m, df, n, rwork )
248 IF( cnorm.GT.zero )
THEN
249 result( 6 ) = resid / (eps*max(1,m)*dnorm)
256 DEALLOCATE ( a, af, q, r, rwork, work, t, c, d, cf, df)
subroutine sgeqrt(M, N, NB, A, LDA, T, LDT, WORK, INFO)
SGEQRT
subroutine ssyrk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
SSYRK
subroutine slarnv(IDIST, ISEED, N, X)
SLARNV returns a vector of random numbers from a uniform or normal distribution.
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
subroutine sgemqrt(SIDE, TRANS, M, N, K, NB, V, LDV, T, LDT, C, LDC, WORK, INFO)
SGEMQRT
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 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 slamch(CMACH)
SLAMCH
logical function lsame(CA, CB)
LSAME
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.