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 zggrqf(M, P, N, A, LDA, TAUA, B, LDB, TAUB, WORK, LWORK, INFO)
ZGGRQF
subroutine zungrq(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGRQ
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
subroutine zungqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGQR
subroutine zgrqts(M, P, N, A, AF, Q, R, LDA, TAUA, B, BF, Z, T, BWK, LDB, TAUB, WORK, LWORK, RWORK, RESULT)
ZGRQTS
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 zherk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
ZHERK