177 SUBROUTINE sgrqts( M, P, N, A, AF, Q, R, LDA, TAUA, B, BF, Z, T,
178 $ bwk, ldb, taub, work, lwork, rwork, result )
186 INTEGER LDA, LDB, LWORK, M, P, N
189 REAL A( lda, * ), AF( lda, * ), R( lda, * ),
191 $ b( ldb, * ), bf( ldb, * ), t( ldb, * ),
192 $ z( ldb, * ), bwk( ldb, * ),
193 $ taua( * ), taub( * ),
194 $ result( 4 ), rwork( * ), work( lwork )
201 parameter ( zero = 0.0e+0, one = 1.0e+0 )
203 parameter ( rogue = -1.0e+10 )
207 REAL ANORM, BNORM, ULP, UNFL, RESID
210 REAL SLAMCH, SLANGE, SLANSY
211 EXTERNAL slamch, slange, slansy
218 INTRINSIC max, min, real
222 ulp = slamch(
'Precision' )
223 unfl = slamch(
'Safe minimum' )
227 CALL slacpy(
'Full', m, n, a, lda, af, lda )
228 CALL slacpy(
'Full', p, n, b, ldb, bf, ldb )
230 anorm = max( slange(
'1', m, n, a, lda, rwork ), unfl )
231 bnorm = max( slange(
'1', p, n, b, ldb, rwork ), unfl )
235 CALL sggrqf( m, p, n, af, lda, taua, bf, ldb, taub, work,
240 CALL slaset(
'Full', n, n, rogue, rogue, q, lda )
242 IF( m.GT.0 .AND. m.LT.n )
243 $
CALL slacpy(
'Full', m, n-m, af, lda, q( n-m+1, 1 ), lda )
245 $
CALL slacpy(
'Lower', m-1, m-1, af( 2, n-m+1 ), lda,
246 $ q( n-m+2, n-m+1 ), lda )
249 $
CALL slacpy(
'Lower', n-1, n-1, af( m-n+2, 1 ), lda,
252 CALL sorgrq( n, n, min( m, n ), q, lda, taua, work, lwork, info )
256 CALL slaset(
'Full', p, p, rogue, rogue, z, ldb )
258 $
CALL slacpy(
'Lower', p-1, n, bf( 2,1 ), ldb, z( 2,1 ), ldb )
259 CALL sorgqr( p, p, min( p,n ), z, ldb, taub, work, lwork, info )
263 CALL slaset(
'Full', m, n, zero, zero, r, lda )
265 CALL slacpy(
'Upper', m, m, af( 1, n-m+1 ), lda, r( 1, n-m+1 ),
268 CALL slacpy(
'Full', m-n, n, af, lda, r, lda )
269 CALL slacpy(
'Upper', n, n, af( m-n+1, 1 ), lda, r( m-n+1, 1 ),
275 CALL slaset(
'Full', p, n, zero, zero, t, ldb )
276 CALL slacpy(
'Upper', p, n, bf, ldb, t, ldb )
280 CALL sgemm(
'No transpose',
'Transpose', m, n, n, -one, a, lda, q,
285 resid = slange(
'1', m, n, r, lda, rwork )
286 IF( anorm.GT.zero )
THEN
287 result( 1 ) = ( ( resid /
REAL(MAX(1,M,N) ) ) / anorm ) / ulp
294 CALL sgemm(
'Transpose',
'No transpose', p, n, p, one, z, ldb, b,
295 $ ldb, zero, bwk, ldb )
296 CALL sgemm(
'No transpose',
'No transpose', p, n, n, one, t, ldb,
297 $ q, lda, -one, bwk, ldb )
301 resid = slange(
'1', p, n, bwk, ldb, rwork )
302 IF( bnorm.GT.zero )
THEN
303 result( 2 ) = ( ( resid /
REAL( MAX( 1,P,M ) ) )/bnorm ) / ulp
310 CALL slaset(
'Full', n, n, zero, one, r, lda )
311 CALL ssyrk(
'Upper',
'No Transpose', n, n, -one, q, lda, one, r,
316 resid = slansy(
'1',
'Upper', n, r, lda, rwork )
317 result( 3 ) = ( resid /
REAL( MAX( 1,N ) ) ) / ulp
321 CALL slaset(
'Full', p, p, zero, one, t, ldb )
322 CALL ssyrk(
'Upper',
'Transpose', p, p, -one, z, ldb, one, t,
327 resid = slansy(
'1',
'Upper', p, t, ldb, rwork )
328 result( 4 ) = ( resid /
REAL( MAX( 1,P ) ) ) / ulp
subroutine ssyrk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
SSYRK
subroutine sggrqf(M, P, N, A, LDA, TAUA, B, LDB, TAUB, WORK, LWORK, INFO)
SGGRQF
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
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 sorgqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
SORGQR
subroutine sgrqts(M, P, N, A, AF, Q, R, LDA, TAUA, B, BF, Z, T, BWK, LDB, TAUB, WORK, LWORK, RWORK, RESULT)
SGRQTS
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
subroutine sorgrq(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
SORGRQ