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
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', 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
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
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 dggqrf(N, M, P, A, LDA, TAUA, B, LDB, TAUB, WORK, LWORK, INFO)
DGGQRF
subroutine dgqrts(N, M, P, A, AF, Q, R, LDA, TAUA, B, BF, Z, T, BWK, LDB, TAUB, WORK, LWORK, RWORK, RESULT)
DGQRTS