176 SUBROUTINE cgqrts( 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 rwork( * ), result( 4 )
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', n, m, a, lda, af, lda )
229 CALL
clacpy(
'Full', n, p,
b, ldb, bf, ldb )
231 anorm = max(
clange(
'1', n, m, a, lda, rwork ), unfl )
232 bnorm = max(
clange(
'1', n, p,
b, ldb, rwork ), unfl )
236 CALL
cggqrf( n, m, p, af, lda, taua, bf, ldb, taub, work,
241 CALL
claset(
'Full', n, n, crogue, crogue, q, lda )
242 CALL
clacpy(
'Lower', n-1, m, af( 2,1 ), lda, q( 2,1 ), lda )
243 CALL
cungqr( n, n, min( n, m ), q, lda, taua, work, lwork, info )
247 CALL
claset(
'Full', p, p, crogue, crogue, z, ldb )
249 IF( n.GT.0 .AND. n.LT.p )
250 $ CALL
clacpy(
'Full', n, p-n, bf, ldb, z( p-n+1, 1 ), ldb )
252 $ CALL
clacpy(
'Lower', n-1, n-1, bf( 2, p-n+1 ), ldb,
253 $ z( p-n+2, p-n+1 ), ldb )
256 $ CALL
clacpy(
'Lower', p-1, p-1, bf( n-p+2, 1 ), ldb,
259 CALL
cungrq( p, p, min( n, p ), z, ldb, taub, work, lwork, info )
263 CALL
claset(
'Full', n, m, czero, czero, r, lda )
264 CALL
clacpy(
'Upper', n, m, af, lda, r, lda )
268 CALL
claset(
'Full', n, p, czero, czero, t, ldb )
270 CALL
clacpy(
'Upper', n, n, bf( 1, p-n+1 ), ldb, t( 1, p-n+1 ),
273 CALL
clacpy(
'Full', n-p, p, bf, ldb, t, ldb )
274 CALL
clacpy(
'Upper', p, p, bf( n-p+1, 1 ), ldb, t( n-p+1, 1 ),
280 CALL
cgemm(
'Conjugate transpose',
'No transpose', n, m, n, -cone,
281 $ q, lda, a, lda, cone, r, lda )
285 resid =
clange(
'1', n, m, r, lda, rwork )
286 IF( anorm.GT.zero )
THEN
287 result( 1 ) = ( ( resid /
REAL( MAX(1,M,N) ) ) / anorm ) / ulp
294 CALL
cgemm(
'No Transpose',
'No transpose', n, p, p, cone, t, ldb,
295 $ z, ldb, czero, bwk, ldb )
296 CALL
cgemm(
'Conjugate transpose',
'No transpose', n, p, n, -cone,
297 $ q, lda,
b, ldb, cone, bwk, ldb )
301 resid =
clange(
'1', n, p, bwk, ldb, rwork )
302 IF( bnorm.GT.zero )
THEN
303 result( 2 ) = ( ( resid /
REAL( MAX(1,P,N ) ) )/bnorm ) / ulp
310 CALL
claset(
'Full', n, n, czero, cone, r, lda )
311 CALL
cherk(
'Upper',
'Conjugate transpose', n, n, -one, q, lda,
316 resid =
clanhe(
'1',
'Upper', n, r, lda, rwork )
317 result( 3 ) = ( resid /
REAL( MAX( 1, N ) ) ) / ulp
321 CALL
claset(
'Full', p, p, czero, cone, t, ldb )
322 CALL
cherk(
'Upper',
'Conjugate transpose', p, p, -one, z, ldb,
327 resid =
clanhe(
'1',
'Upper', p, t, ldb, rwork )
328 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 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 cggqrf(N, M, P, A, LDA, TAUA, B, LDB, TAUB, WORK, LWORK, INFO)
CGGQRF
REAL function slamch(CMACH)
SLAMCH
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
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 ...
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
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
subroutine cgqrts(N, M, P, A, AF, Q, R, LDA, TAUA, B, BF, Z, T, BWK, LDB, TAUB, WORK, LWORK, RWORK, RESULT)
CGQRTS