88 REAL,
ALLOCATABLE :: AF(:,:), Q(:,:),
89 $ R(:,:), RWORK(:), WORK( : ), T(:,:),
90 $ CF(:,:), DF(:,:), A(:,:), C(:,:), D(:,:)
94 parameter( zero = 0.0, one = 1.0 )
97 INTEGER INFO, J, K, L, LWORK
98 REAL ANORM, EPS, RESID, CNORM, DNORM
110 EXTERNAL slamch, slange, slansy, lsame
116 DATA iseed / 1988, 1989, 1990, 1991 /
118 eps = slamch(
'Epsilon' )
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)