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
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 dgqrts(N, M, P, A, AF, Q, R, LDA, TAUA, B, BF, Z, T, BWK, LDB, TAUB, WORK, LWORK, RWORK, RESULT)
DGQRTS
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 dggqrf(N, M, P, A, LDA, TAUA, B, LDB, TAUB, WORK, LWORK, INFO)
DGGQRF
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 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