176 SUBROUTINE zgqrts( 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, 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', n, m, a, lda, af, lda )
229 CALL zlacpy(
'Full', n, p, b, ldb, bf, ldb )
231 anorm = max( zlange(
'1', n, m, a, lda, rwork ), unfl )
232 bnorm = max( zlange(
'1', n, p, b, ldb, rwork ), unfl )
236 CALL zggqrf( n, m, p, af, lda, taua, bf, ldb, taub, work, lwork,
241 CALL zlaset(
'Full', n, n, crogue, crogue, q, lda )
242 CALL zlacpy(
'Lower', n-1, m, af( 2, 1 ), lda, q( 2, 1 ), lda )
243 CALL zungqr( n, n, min( n, m ), q, lda, taua, work, lwork, info )
247 CALL zlaset(
'Full', p, p, crogue, crogue, z, ldb )
249 IF( n.GT.0 .AND. n.LT.p )
250 $
CALL zlacpy(
'Full', n, p-n, bf, ldb, z( p-n+1, 1 ), ldb )
252 $
CALL zlacpy(
'Lower', n-1, n-1, bf( 2, p-n+1 ), ldb,
253 $ z( p-n+2, p-n+1 ), ldb )
256 $
CALL zlacpy(
'Lower', p-1, p-1, bf( n-p+2, 1 ), ldb,
259 CALL zungrq( p, p, min( n, p ), z, ldb, taub, work, lwork, info )
263 CALL zlaset(
'Full', n, m, czero, czero, r, lda )
264 CALL zlacpy(
'Upper', n, m, af, lda, r, lda )
268 CALL zlaset(
'Full', n, p, czero, czero, t, ldb )
270 CALL zlacpy(
'Upper', n, n, bf( 1, p-n+1 ), ldb, t( 1, p-n+1 ),
273 CALL zlacpy(
'Full', n-p, p, bf, ldb, t, ldb )
274 CALL zlacpy(
'Upper', p, p, bf( n-p+1, 1 ), ldb, t( n-p+1, 1 ),
280 CALL zgemm(
'Conjugate transpose',
'No transpose', n, m, n, -cone,
281 $ q, lda, a, lda, cone, r, lda )
285 resid = zlange(
'1', n, m, r, lda, rwork )
286 IF( anorm.GT.zero )
THEN
287 result( 1 ) = ( ( resid / dble( max( 1, m, n ) ) ) / anorm ) /
295 CALL zgemm(
'No Transpose',
'No transpose', n, p, p, cone, t, ldb,
296 $ z, ldb, czero, bwk, ldb )
297 CALL zgemm(
'Conjugate transpose',
'No transpose', n, p, n, -cone,
298 $ q, lda, b, ldb, cone, bwk, ldb )
302 resid = zlange(
'1', n, p, bwk, ldb, rwork )
303 IF( bnorm.GT.zero )
THEN
304 result( 2 ) = ( ( resid / dble( max( 1, p, n ) ) ) / bnorm ) /
312 CALL zlaset(
'Full', n, n, czero, cone, r, lda )
313 CALL zherk(
'Upper',
'Conjugate transpose', n, n, -one, q, lda,
318 resid = zlanhe(
'1',
'Upper', n, r, lda, rwork )
319 result( 3 ) = ( resid / dble( max( 1, n ) ) ) / ulp
323 CALL zlaset(
'Full', p, p, czero, cone, t, ldb )
324 CALL zherk(
'Upper',
'Conjugate transpose', p, p, -one, z, ldb,
329 resid = zlanhe(
'1',
'Upper', p, t, ldb, rwork )
330 result( 4 ) = ( resid / dble( max( 1, p ) ) ) / ulp
subroutine zggqrf(N, M, P, A, LDA, TAUA, B, LDB, TAUB, WORK, LWORK, INFO)
ZGGQRF
subroutine zungrq(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGRQ
subroutine zgqrts(N, M, P, A, AF, Q, R, LDA, TAUA, B, BF, Z, T, BWK, LDB, TAUB, WORK, LWORK, RWORK, RESULT)
ZGQRTS
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 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