176 SUBROUTINE cgrqts( M, P, N, 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 result( 4 ), rwork( * )
189 COMPLEX a( lda, * ), af( lda, * ), r( lda, * ),
190 $ q( lda, * ),
b( ldb, * ), bf( ldb, * ),
191 $ t( ldb, * ), z( ldb, * ), bwk( ldb, * ),
192 $ taua( * ), taub( * ), work( lwork )
199 parameter( zero = 0.0e+0, one = 1.0e+0 )
201 parameter( czero = ( 0.0e+0, 0.0e+0 ),
202 $ cone = ( 1.0e+0, 0.0e+0 ) )
204 parameter( crogue = ( -1.0e+10, 0.0e+0 ) )
208 REAL anorm, bnorm, ulp, unfl, resid
219 INTRINSIC max, min, real
223 ulp =
slamch(
'Precision' )
224 unfl =
slamch(
'Safe minimum' )
228 CALL
clacpy(
'Full', m, n, a, lda, af, lda )
229 CALL
clacpy(
'Full', p, n,
b, ldb, bf, ldb )
231 anorm = max(
clange(
'1', m, n, a, lda, rwork ), unfl )
232 bnorm = max(
clange(
'1', p, n,
b, ldb, rwork ), unfl )
236 CALL
cggrqf( m, p, n, af, lda, taua, bf, ldb, taub, work,
241 CALL
claset(
'Full', n, n, crogue, crogue, q, lda )
243 IF( m.GT.0 .AND. m.LT.n )
244 $ CALL
clacpy(
'Full', m, n-m, af, lda, q( n-m+1, 1 ), lda )
246 $ CALL
clacpy(
'Lower', m-1, m-1, af( 2, n-m+1 ), lda,
247 $ q( n-m+2, n-m+1 ), lda )
250 $ CALL
clacpy(
'Lower', n-1, n-1, af( m-n+2, 1 ), lda,
253 CALL
cungrq( n, n, min( m, n ), q, lda, taua, work, lwork, info )
257 CALL
claset(
'Full', p, p, crogue, crogue, z, ldb )
259 $ CALL
clacpy(
'Lower', p-1, n, bf( 2,1 ), ldb, z( 2,1 ), ldb )
260 CALL
cungqr( p, p, min( p,n ), z, ldb, taub, work, lwork, info )
264 CALL
claset(
'Full', m, n, czero, czero, r, lda )
266 CALL
clacpy(
'Upper', m, m, af( 1, n-m+1 ), lda, r( 1, n-m+1 ),
269 CALL
clacpy(
'Full', m-n, n, af, lda, r, lda )
270 CALL
clacpy(
'Upper', n, n, af( m-n+1, 1 ), lda, r( m-n+1, 1 ),
276 CALL
claset(
'Full', p, n, czero, czero, t, ldb )
277 CALL
clacpy(
'Upper', p, n, bf, ldb, t, ldb )
281 CALL
cgemm(
'No transpose',
'Conjugate transpose', m, n, n, -cone,
282 $ a, lda, q, lda, cone, r, lda )
286 resid =
clange(
'1', m, n, r, lda, rwork )
287 IF( anorm.GT.zero )
THEN
288 result( 1 ) = ( ( resid /
REAL(MAX(1,M,N) ) ) / anorm ) / ulp
295 CALL
cgemm(
'Conjugate transpose',
'No transpose', p, n, p, cone,
296 $ z, ldb,
b, ldb, czero, bwk, ldb )
297 CALL
cgemm(
'No transpose',
'No transpose', p, n, n, cone, t, ldb,
298 $ q, lda, -cone, bwk, ldb )
302 resid =
clange(
'1', p, n, bwk, ldb, rwork )
303 IF( bnorm.GT.zero )
THEN
304 result( 2 ) = ( ( resid /
REAL( MAX( 1,P,M ) ) )/bnorm ) / ulp
311 CALL
claset(
'Full', n, n, czero, cone, r, lda )
312 CALL
cherk(
'Upper',
'No Transpose', n, n, -one, q, lda, one, r,
317 resid =
clanhe(
'1',
'Upper', n, r, lda, rwork )
318 result( 3 ) = ( resid /
REAL( MAX( 1,N ) ) ) / ulp
322 CALL
claset(
'Full', p, p, czero, cone, t, ldb )
323 CALL
cherk(
'Upper',
'Conjugate transpose', p, p, -one, z, ldb,
328 resid =
clanhe(
'1',
'Upper', p, t, ldb, rwork )
329 result( 4 ) = ( resid /
REAL( MAX( 1,P ) ) ) / ulp
subroutine claset(UPLO, M, N, ALPHA, BETA, A, LDA)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
real function clange(NORM, M, N, A, LDA, WORK)
CLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
real function clanhe(NORM, UPLO, N, A, LDA, WORK)
CLANHE returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a complex Hermitian matrix.
subroutine cggrqf(M, P, N, A, LDA, TAUA, B, LDB, TAUB, WORK, LWORK, INFO)
CGGRQF
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real b(3) integer i
subroutine cungrq(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
CUNGRQ
subroutine cgrqts(M, P, N, A, AF, Q, R, LDA, TAUA, B, BF, Z, T, BWK, LDB, TAUB, WORK, LWORK, RWORK, RESULT)
CGRQTS
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
real function slamch(CMACH)
SLAMCH
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM
subroutine cherk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
CHERK
subroutine cungqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
CUNGQR