176 SUBROUTINE zgrqts( 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, N, P
188 DOUBLE PRECISION RESULT( 4 ), RWORK( * )
189 COMPLEX*16 A( lda, * ), AF( lda, * ), B( ldb, * ),
190 $ bf( ldb, * ), bwk( ldb, * ), q( lda, * ),
191 $ r( lda, * ), t( ldb, * ), taua( * ), taub( * ),
192 $ work( lwork ), z( ldb, * )
198 DOUBLE PRECISION ZERO, ONE
199 parameter( zero = 0.0d+0, one = 1.0d+0 )
200 COMPLEX*16 CZERO, CONE
201 parameter( czero = ( 0.0d+0, 0.0d+0 ),
202 $ cone = ( 1.0d+0, 0.0d+0 ) )
204 parameter( crogue = ( -1.0d+10, 0.0d+0 ) )
208 DOUBLE PRECISION ANORM, BNORM, RESID, ULP, UNFL
211 DOUBLE PRECISION DLAMCH, ZLANGE, ZLANHE
212 EXTERNAL dlamch, zlange, zlanhe
219 INTRINSIC dble, max, min
223 ulp = dlamch(
'Precision' )
224 unfl = dlamch(
'Safe minimum' )
228 CALL
zlacpy(
'Full', m, n, a, lda, af, lda )
229 CALL
zlacpy(
'Full', p, n, b, ldb, bf, ldb )
231 anorm = max( zlange(
'1', m, n, a, lda, rwork ), unfl )
232 bnorm = max( zlange(
'1', p, n, b, ldb, rwork ), unfl )
236 CALL
zggrqf( m, p, n, af, lda, taua, bf, ldb, taub, work, lwork,
241 CALL
zlaset(
'Full', n, n, crogue, crogue, q, lda )
243 IF( m.GT.0 .AND. m.LT.n )
244 $ CALL
zlacpy(
'Full', m, n-m, af, lda, q( n-m+1, 1 ), lda )
246 $ CALL
zlacpy(
'Lower', m-1, m-1, af( 2, n-m+1 ), lda,
247 $ q( n-m+2, n-m+1 ), lda )
250 $ CALL
zlacpy(
'Lower', n-1, n-1, af( m-n+2, 1 ), lda,
253 CALL
zungrq( n, n, min( m, n ), q, lda, taua, work, lwork, info )
257 CALL
zlaset(
'Full', p, p, crogue, crogue, z, ldb )
259 $ CALL
zlacpy(
'Lower', p-1, n, bf( 2, 1 ), ldb, z( 2, 1 ), ldb )
260 CALL
zungqr( p, p, min( p, n ), z, ldb, taub, work, lwork, info )
264 CALL
zlaset(
'Full', m, n, czero, czero, r, lda )
266 CALL
zlacpy(
'Upper', m, m, af( 1, n-m+1 ), lda, r( 1, n-m+1 ),
269 CALL
zlacpy(
'Full', m-n, n, af, lda, r, lda )
270 CALL
zlacpy(
'Upper', n, n, af( m-n+1, 1 ), lda, r( m-n+1, 1 ),
276 CALL
zlaset(
'Full', p, n, czero, czero, t, ldb )
277 CALL
zlacpy(
'Upper', p, n, bf, ldb, t, ldb )
281 CALL
zgemm(
'No transpose',
'Conjugate transpose', m, n, n, -cone,
282 $ a, lda, q, lda, cone, r, lda )
286 resid = zlange(
'1', m, n, r, lda, rwork )
287 IF( anorm.GT.zero )
THEN
288 result( 1 ) = ( ( resid / dble( max( 1, m, n ) ) ) / anorm ) /
296 CALL
zgemm(
'Conjugate transpose',
'No transpose', p, n, p, cone,
297 $ z, ldb, b, ldb, czero, bwk, ldb )
298 CALL
zgemm(
'No transpose',
'No transpose', p, n, n, cone, t, ldb,
299 $ q, lda, -cone, bwk, ldb )
303 resid = zlange(
'1', p, n, bwk, ldb, rwork )
304 IF( bnorm.GT.zero )
THEN
305 result( 2 ) = ( ( resid / dble( max( 1, p, m ) ) ) / bnorm ) /
313 CALL
zlaset(
'Full', n, n, czero, cone, r, lda )
314 CALL
zherk(
'Upper',
'No Transpose', n, n, -one, q, lda, one, r,
319 resid = zlanhe(
'1',
'Upper', n, r, lda, rwork )
320 result( 3 ) = ( resid / dble( max( 1, n ) ) ) / ulp
324 CALL
zlaset(
'Full', p, p, czero, cone, t, ldb )
325 CALL
zherk(
'Upper',
'Conjugate transpose', p, p, -one, z, ldb,
330 resid = zlanhe(
'1',
'Upper', p, t, ldb, rwork )
331 result( 4 ) = ( resid / dble( max( 1, p ) ) ) / ulp
subroutine zherk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
ZHERK
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
subroutine zungqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGQR
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
subroutine zungrq(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGRQ
subroutine zgrqts(M, P, N, A, AF, Q, R, LDA, TAUA, B, BF, Z, T, BWK, LDB, TAUB, WORK, LWORK, RWORK, RESULT)
ZGRQTS
subroutine zggrqf(M, P, N, A, LDA, TAUA, B, LDB, TAUB, WORK, LWORK, INFO)
ZGGRQF