176 SUBROUTINE sgqrts( 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, P, N
188 REAL A( lda, * ), AF( lda, * ), R( lda, * ),
189 $ q( lda, * ), b( ldb, * ), bf( ldb, * ),
190 $ t( ldb, * ), z( ldb, * ), bwk( ldb, * ),
191 $ taua( * ), taub( * ), result( 4 ),
192 $ rwork( * ), work( lwork )
199 parameter( zero = 0.0e+0, one = 1.0e+0 )
201 parameter( rogue = -1.0e+10 )
205 REAL ANORM, BNORM, ULP, UNFL, RESID
208 REAL SLAMCH, SLANGE, SLANSY
209 EXTERNAL slamch, slange, slansy
216 INTRINSIC max, min, real
220 ulp = slamch(
'Precision' )
221 unfl = slamch(
'Safe minimum' )
225 CALL
slacpy(
'Full', n, m, a, lda, af, lda )
226 CALL
slacpy(
'Full', n, p, b, ldb, bf, ldb )
228 anorm = max( slange(
'1', n, m, a, lda, rwork ), unfl )
229 bnorm = max( slange(
'1', n, p, b, ldb, rwork ), unfl )
233 CALL
sggqrf( n, m, p, af, lda, taua, bf, ldb, taub, work,
238 CALL
slaset(
'Full', n, n, rogue, rogue, q, lda )
239 CALL
slacpy(
'Lower', n-1, m, af( 2,1 ), lda, q( 2,1 ), lda )
240 CALL
sorgqr( n, n, min( n, m ), q, lda, taua, work, lwork, info )
244 CALL
slaset(
'Full', p, p, rogue, rogue, z, ldb )
246 IF( n.GT.0 .AND. n.LT.p )
247 $ CALL
slacpy(
'Full', n, p-n, bf, ldb, z( p-n+1, 1 ), ldb )
249 $ CALL
slacpy(
'Lower', n-1, n-1, bf( 2, p-n+1 ), ldb,
250 $ z( p-n+2, p-n+1 ), ldb )
253 $ CALL
slacpy(
'Lower', p-1, p-1, bf( n-p+2, 1 ), ldb,
256 CALL
sorgrq( p, p, min( n, p ), z, ldb, taub, work, lwork, info )
260 CALL
slaset(
'Full', n, m, zero, zero, r, lda )
261 CALL
slacpy(
'Upper', n, m, af, lda, r, lda )
265 CALL
slaset(
'Full', n, p, zero, zero, t, ldb )
267 CALL
slacpy(
'Upper', n, n, bf( 1, p-n+1 ), ldb, t( 1, p-n+1 ),
270 CALL
slacpy(
'Full', n-p, p, bf, ldb, t, ldb )
271 CALL
slacpy(
'Upper', p, p, bf( n-p+1, 1 ), ldb, t( n-p+1, 1 ),
277 CALL
sgemm(
'Transpose',
'No transpose', n, m, n, -one, q, lda, a,
282 resid = slange(
'1', n, m, r, lda, rwork )
283 IF( anorm.GT.zero )
THEN
284 result( 1 ) = ( ( resid /
REAL( MAX(1,M,N) ) ) / anorm ) / ulp
291 CALL
sgemm(
'No Transpose',
'No transpose', n, p, p, one, t, ldb,
292 $ z, ldb, zero, bwk, ldb )
293 CALL
sgemm(
'Transpose',
'No transpose', n, p, n, -one, q, lda,
294 $ b, ldb, one, bwk, ldb )
298 resid = slange(
'1', n, p, bwk, ldb, rwork )
299 IF( bnorm.GT.zero )
THEN
300 result( 2 ) = ( ( resid /
REAL( MAX(1,P,N ) ) )/bnorm ) / ulp
307 CALL
slaset(
'Full', n, n, zero, one, r, lda )
308 CALL
ssyrk(
'Upper',
'Transpose', n, n, -one, q, lda, one, r,
313 resid = slansy(
'1',
'Upper', n, r, lda, rwork )
314 result( 3 ) = ( resid /
REAL( MAX( 1, N ) ) ) / ulp
318 CALL
slaset(
'Full', p, p, zero, one, t, ldb )
319 CALL
ssyrk(
'Upper',
'Transpose', p, p, -one, z, ldb, one, t,
324 resid = slansy(
'1',
'Upper', p, t, ldb, rwork )
325 result( 4 ) = ( resid /
REAL( MAX( 1, P ) ) ) / ulp
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine ssyrk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
SSYRK
subroutine sorgrq(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
SORGRQ
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
subroutine sggqrf(N, M, P, A, LDA, TAUA, B, LDB, TAUB, WORK, LWORK, INFO)
SGGQRF
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
subroutine sorgqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
SORGQR
subroutine sgqrts(N, M, P, A, AF, Q, R, LDA, TAUA, B, BF, Z, T, BWK, LDB, TAUB, WORK, LWORK, RWORK, RESULT)
SGQRTS