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 ssyrk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
SSYRK
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
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 sggqrf(N, M, P, A, LDA, TAUA, B, LDB, TAUB, WORK, LWORK, INFO)
SGGQRF
subroutine sgqrts(N, M, P, A, AF, Q, R, LDA, TAUA, B, BF, Z, T, BWK, LDB, TAUB, WORK, LWORK, RWORK, RESULT)
SGQRTS
subroutine sorgqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
SORGQR
subroutine sorgrq(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
SORGRQ