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
DOUBLE PRECISION function dlansy(NORM, UPLO, N, A, LDA, WORK)
DLANSY returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a real symmetric matrix.
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
subroutine dggrqf(M, P, N, A, LDA, TAUA, B, LDB, TAUB, WORK, LWORK, INFO)
DGGRQF
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real b(3) integer i
DOUBLE PRECISION function dlange(NORM, M, N, A, LDA, WORK)
DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
subroutine dorgqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
DORGQR
subroutine dsyrk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
DSYRK
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
DOUBLE PRECISION function dlamch(CMACH)
DLAMCH
subroutine dgrqts(M, P, N, A, AF, Q, R, LDA, TAUA, B, BF, Z, T, BWK, LDB, TAUB, WORK, LWORK, RWORK, RESULT)
DGRQTS
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 dorgrq(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
DORGRQ