176 SUBROUTINE dgqrts( N, M, P, 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', n, m, a, lda, af, lda )
226 CALL
dlacpy(
'Full', n, p, b, ldb, bf, ldb )
228 anorm = max(
dlange(
'1', n, m, a, lda, rwork ), unfl )
229 bnorm = max(
dlange(
'1', n, p, b, ldb, rwork ), unfl )
233 CALL
dggqrf( n, m, p, af, lda, taua, bf, ldb, taub, work, lwork,
238 CALL
dlaset(
'Full', n, n, rogue, rogue, q, lda )
239 CALL
dlacpy(
'Lower', n-1, m, af( 2, 1 ), lda, q( 2, 1 ), lda )
240 CALL
dorgqr( n, n, min( n, m ), q, lda, taua, work, lwork, info )
244 CALL
dlaset(
'Full', p, p, rogue, rogue, z, ldb )
246 IF( n.GT.0 .AND. n.LT.p )
247 $ CALL
dlacpy(
'Full', n, p-n, bf, ldb, z( p-n+1, 1 ), ldb )
249 $ CALL
dlacpy(
'Lower', n-1, n-1, bf( 2, p-n+1 ), ldb,
250 $ z( p-n+2, p-n+1 ), ldb )
253 $ CALL
dlacpy(
'Lower', p-1, p-1, bf( n-p+2, 1 ), ldb,
256 CALL
dorgrq( p, p, min( n, p ), z, ldb, taub, work, lwork, info )
260 CALL
dlaset(
'Full', n, m, zero, zero, r, lda )
261 CALL
dlacpy(
'Upper', n, m, af, lda, r, lda )
265 CALL
dlaset(
'Full', n, p, zero, zero, t, ldb )
267 CALL
dlacpy(
'Upper', n, n, bf( 1, p-n+1 ), ldb, t( 1, p-n+1 ),
270 CALL
dlacpy(
'Full', n-p, p, bf, ldb, t, ldb )
271 CALL
dlacpy(
'Upper', p, p, bf( n-p+1, 1 ), ldb, t( n-p+1, 1 ),
277 CALL
dgemm(
'Transpose',
'No transpose', n, m, n, -one, q, lda, a,
282 resid =
dlange(
'1', n, m, r, lda, rwork )
283 IF( anorm.GT.zero )
THEN
284 result( 1 ) = ( ( resid / dble( max( 1, m, n ) ) ) / anorm ) /
292 CALL
dgemm(
'No Transpose',
'No transpose', n, p, p, one, t, ldb,
293 $ z, ldb, zero, bwk, ldb )
294 CALL
dgemm(
'Transpose',
'No transpose', n, p, n, -one, q, lda, b,
295 $ ldb, one, bwk, ldb )
299 resid =
dlange(
'1', n, p, bwk, ldb, rwork )
300 IF( bnorm.GT.zero )
THEN
301 result( 2 ) = ( ( resid / dble( max( 1, p, n ) ) ) / bnorm ) /
309 CALL
dlaset(
'Full', n, n, zero, one, r, lda )
310 CALL
dsyrk(
'Upper',
'Transpose', n, n, -one, q, lda, one, r,
315 resid =
dlansy(
'1',
'Upper', n, r, lda, rwork )
316 result( 3 ) = ( resid / dble( max( 1, n ) ) ) / ulp
320 CALL
dlaset(
'Full', p, p, zero, one, t, ldb )
321 CALL
dsyrk(
'Upper',
'Transpose', p, p, -one, z, ldb, one, t,
326 resid =
dlansy(
'1',
'Upper', p, t, ldb, rwork )
327 result( 4 ) = ( resid / dble( max( 1, p ) ) ) / ulp