174 SUBROUTINE zgqrts( N, M, P, A, AF, Q, R, LDA, TAUA, B, BF, Z, T,
175 $ BWK, LDB, TAUB, WORK, LWORK, RWORK, RESULT )
182 INTEGER LDA, LDB, LWORK, M, N, P
185 DOUBLE PRECISION RESULT( 4 ), RWORK( * )
186 COMPLEX*16 A( LDA, * ), AF( LDA, * ), B( LDB, * ),
187 $ bf( ldb, * ), bwk( ldb, * ), q( lda, * ),
188 $ r( lda, * ), t( ldb, * ), taua( * ), taub( * ),
189 $ work( lwork ), z( ldb, * )
195 DOUBLE PRECISION ZERO, ONE
196 parameter( zero = 0.0d+0, one = 1.0d+0 )
197 COMPLEX*16 CZERO, CONE
198 parameter( czero = ( 0.0d+0, 0.0d+0 ),
199 $ cone = ( 1.0d+0, 0.0d+0 ) )
201 parameter( crogue = ( -1.0d+10, 0.0d+0 ) )
205 DOUBLE PRECISION ANORM, BNORM, RESID, ULP, UNFL
208 DOUBLE PRECISION DLAMCH, ZLANGE, ZLANHE
209 EXTERNAL dlamch, zlange, zlanhe
216 INTRINSIC dble, max, min
220 ulp = dlamch(
'Precision' )
221 unfl = dlamch(
'Safe minimum' )
225 CALL zlacpy(
'Full', n, m, a, lda, af, lda )
226 CALL zlacpy(
'Full', n, p, b, ldb, bf, ldb )
228 anorm = max( zlange(
'1', n, m, a, lda, rwork ), unfl )
229 bnorm = max( zlange(
'1', n, p, b, ldb, rwork ), unfl )
233 CALL zggqrf( n, m, p, af, lda, taua, bf, ldb, taub, work, lwork,
238 CALL zlaset(
'Full', n, n, crogue, crogue, q, lda )
239 CALL zlacpy(
'Lower', n-1, m, af( 2, 1 ), lda, q( 2, 1 ), lda )
240 CALL zungqr( n, n, min( n, m ), q, lda, taua, work, lwork, info )
244 CALL zlaset(
'Full', p, p, crogue, crogue, z, ldb )
246 IF( n.GT.0 .AND. n.LT.p )
247 $
CALL zlacpy(
'Full', n, p-n, bf, ldb, z( p-n+1, 1 ), ldb )
249 $
CALL zlacpy(
'Lower', n-1, n-1, bf( 2, p-n+1 ), ldb,
250 $ z( p-n+2, p-n+1 ), ldb )
253 $
CALL zlacpy(
'Lower', p-1, p-1, bf( n-p+2, 1 ), ldb,
256 CALL zungrq( p, p, min( n, p ), z, ldb, taub, work, lwork, info )
260 CALL zlaset(
'Full', n, m, czero, czero, r, lda )
261 CALL zlacpy(
'Upper', n, m, af, lda, r, lda )
265 CALL zlaset(
'Full', n, p, czero, czero, t, ldb )
267 CALL zlacpy(
'Upper', n, n, bf( 1, p-n+1 ), ldb, t( 1, p-n+1 ),
270 CALL zlacpy(
'Full', n-p, p, bf, ldb, t, ldb )
271 CALL zlacpy(
'Upper', p, p, bf( n-p+1, 1 ), ldb, t( n-p+1, 1 ),
277 CALL zgemm(
'Conjugate transpose',
'No transpose', n, m, n, -cone,
278 $ q, lda, a, lda, cone, r, lda )
282 resid = zlange(
'1', n, m, r, lda, rwork )
283 IF( anorm.GT.zero )
THEN
284 result( 1 ) = ( ( resid / dble( max( 1, m, n ) ) ) / anorm ) /
292 CALL zgemm(
'No Transpose',
'No transpose', n, p, p, cone, t, ldb,
293 $ z, ldb, czero, bwk, ldb )
294 CALL zgemm(
'Conjugate transpose',
'No transpose', n, p, n, -cone,
295 $ q, lda, b, ldb, cone, bwk, ldb )
299 resid = zlange(
'1', n, p, bwk, ldb, rwork )
300 IF( bnorm.GT.zero )
THEN
301 result( 2 ) = ( ( resid / dble( max( 1, p, n ) ) ) / bnorm ) /
309 CALL zlaset(
'Full', n, n, czero, cone, r, lda )
310 CALL zherk(
'Upper',
'Conjugate transpose', n, n, -one, q, lda,
315 resid = zlanhe(
'1',
'Upper', n, r, lda, rwork )
316 result( 3 ) = ( resid / dble( max( 1, n ) ) ) / ulp
320 CALL zlaset(
'Full', p, p, czero, cone, t, ldb )
321 CALL zherk(
'Upper',
'Conjugate transpose', p, p, -one, z, ldb,
326 resid = zlanhe(
'1',
'Upper', p, t, ldb, rwork )
327 result( 4 ) = ( resid / dble( max( 1, p ) ) ) / ulp
subroutine zgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
ZGEMM
subroutine zggqrf(n, m, p, a, lda, taua, b, ldb, taub, work, lwork, info)
ZGGQRF
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 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 zungqr(m, n, k, a, lda, tau, work, lwork, info)
ZUNGQR
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