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
211 REAL SLAMCH, CLANGE, CLANHE
212 EXTERNAL slamch, clange, clanhe
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 cggrqf(M, P, N, A, LDA, TAUA, B, LDB, TAUB, WORK, LWORK, INFO)
CGGRQF
subroutine cungrq(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
CUNGRQ
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...
subroutine cungqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
CUNGQR
subroutine cherk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
CHERK
subroutine clacpy(UPLO, M, N, A, LDA, B, LDB)
CLACPY copies all or part of one two-dimensional array to another.
subroutine cgrqts(M, P, N, A, AF, Q, R, LDA, TAUA, B, BF, Z, T, BWK, LDB, TAUB, WORK, LWORK, RWORK, RESULT)
CGRQTS
subroutine cgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
CGEMM