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
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...
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 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 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
real function slamch(CMACH)
SLAMCH
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.
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