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
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 slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
REAL function slansy(NORM, UPLO, N, A, LDA, WORK)
SLANSY 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.
REAL function slamch(CMACH)
SLAMCH
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
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real b(3) integer i
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.
REAL function slange(NORM, M, N, A, LDA, WORK)
SLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
subroutine sgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
SGEMM
subroutine sggrqf(M, P, N, A, LDA, TAUA, B, LDB, TAUB, WORK, LWORK, INFO)
SGGRQF
subroutine sorgqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
SORGQR