176 SUBROUTINE dgrqts( M, P, N, A, AF, Q, R, LDA, TAUA, B, BF, Z, T,
177 $ bwk, ldb, taub, work, lwork, rwork, result )
185 INTEGER lda, ldb, lwork, m, n, p
188 DOUBLE PRECISION a( lda, * ), af( lda, * ), b( ldb, * ),
189 $ bf( ldb, * ), bwk( ldb, * ), q( lda, * ),
190 $ r( lda, * ), result( 4 ), rwork( * ),
191 $ t( ldb, * ), taua( * ), taub( * ),
192 $ work( lwork ), z( ldb, * )
198 DOUBLE PRECISION zero, one
199 parameter( zero = 0.0d+0, one = 1.0d+0 )
200 DOUBLE PRECISION rogue
201 parameter( rogue = -1.0d+10 )
205 DOUBLE PRECISION anorm, bnorm, resid, ulp, unfl
216 INTRINSIC dble, max, min
220 ulp =
dlamch(
'Precision' )
221 unfl =
dlamch(
'Safe minimum' )
225 CALL
dlacpy(
'Full', m, n, a, lda, af, lda )
226 CALL
dlacpy(
'Full', p, n, b, ldb, bf, ldb )
228 anorm = max(
dlange(
'1', m, n, a, lda, rwork ), unfl )
229 bnorm = max(
dlange(
'1', p, n, b, ldb, rwork ), unfl )
233 CALL
dggrqf( m, p, n, af, lda, taua, bf, ldb, taub, work, lwork,
238 CALL
dlaset(
'Full', n, n, rogue, rogue, q, lda )
240 IF( m.GT.0 .AND. m.LT.n )
241 $ CALL
dlacpy(
'Full', m, n-m, af, lda, q( n-m+1, 1 ), lda )
243 $ CALL
dlacpy(
'Lower', m-1, m-1, af( 2, n-m+1 ), lda,
244 $ q( n-m+2, n-m+1 ), lda )
247 $ CALL
dlacpy(
'Lower', n-1, n-1, af( m-n+2, 1 ), lda,
250 CALL
dorgrq( n, n, min( m, n ), q, lda, taua, work, lwork, info )
254 CALL
dlaset(
'Full', p, p, rogue, rogue, z, ldb )
256 $ CALL
dlacpy(
'Lower', p-1, n, bf( 2, 1 ), ldb, z( 2, 1 ), ldb )
257 CALL
dorgqr( p, p, min( p, n ), z, ldb, taub, work, lwork, info )
261 CALL
dlaset(
'Full', m, n, zero, zero, r, lda )
263 CALL
dlacpy(
'Upper', m, m, af( 1, n-m+1 ), lda, r( 1, n-m+1 ),
266 CALL
dlacpy(
'Full', m-n, n, af, lda, r, lda )
267 CALL
dlacpy(
'Upper', n, n, af( m-n+1, 1 ), lda, r( m-n+1, 1 ),
273 CALL
dlaset(
'Full', p, n, zero, zero, t, ldb )
274 CALL
dlacpy(
'Upper', p, n, bf, ldb, t, ldb )
278 CALL
dgemm(
'No transpose',
'Transpose', m, n, n, -one, a, lda, q,
283 resid =
dlange(
'1', m, n, r, lda, rwork )
284 IF( anorm.GT.zero )
THEN
285 result( 1 ) = ( ( resid / dble( max( 1, m, n ) ) ) / anorm ) /
293 CALL
dgemm(
'Transpose',
'No transpose', p, n, p, one, z, ldb, b,
294 $ ldb, zero, bwk, ldb )
295 CALL
dgemm(
'No transpose',
'No transpose', p, n, n, one, t, ldb,
296 $ q, lda, -one, bwk, ldb )
300 resid =
dlange(
'1', p, n, bwk, ldb, rwork )
301 IF( bnorm.GT.zero )
THEN
302 result( 2 ) = ( ( resid / dble( max( 1, p, m ) ) ) / bnorm ) /
310 CALL
dlaset(
'Full', n, n, zero, one, r, lda )
311 CALL
dsyrk(
'Upper',
'No Transpose', n, n, -one, q, lda, one, r,
316 resid =
dlansy(
'1',
'Upper', n, r, lda, rwork )
317 result( 3 ) = ( resid / dble( max( 1, n ) ) ) / ulp
321 CALL
dlaset(
'Full', p, p, zero, one, t, ldb )
322 CALL
dsyrk(
'Upper',
'Transpose', p, p, -one, z, ldb, one, t,
327 resid =
dlansy(
'1',
'Upper', p, t, ldb, rwork )
328 result( 4 ) = ( resid / dble( max( 1, p ) ) ) / ulp