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
208 DOUBLE PRECISION DLAMCH, DLANGE, DLANSY
209 EXTERNAL dlamch, dlange, dlansy
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
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
subroutine dggrqf(M, P, N, A, LDA, TAUA, B, LDB, TAUB, WORK, LWORK, INFO)
DGGRQF
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
subroutine dsyrk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
DSYRK
subroutine dorgqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
DORGQR
subroutine dorgrq(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
DORGRQ
subroutine dgrqts(M, P, N, A, AF, Q, R, LDA, TAUA, B, BF, Z, T, BWK, LDB, TAUB, WORK, LWORK, RWORK, RESULT)
DGRQTS