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
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