174 SUBROUTINE zgrqts( M, P, N, 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', m, n, a, lda, af, lda )
226 CALL zlacpy(
'Full', p, n, b, ldb, bf, ldb )
228 anorm = max( zlange(
'1', m, n, a, lda, rwork ), unfl )
229 bnorm = max( zlange(
'1', p, n, b, ldb, rwork ), unfl )
233 CALL zggrqf( m, p, n, af, lda, taua, bf, ldb, taub, work, lwork,
238 CALL zlaset(
'Full', n, n, crogue, crogue, q, lda )
240 IF( m.GT.0 .AND. m.LT.n )
241 $
CALL zlacpy(
'Full', m, n-m, af, lda, q( n-m+1, 1 ), lda )
243 $
CALL zlacpy(
'Lower', m-1, m-1, af( 2, n-m+1 ), lda,
244 $ q( n-m+2, n-m+1 ), lda )
247 $
CALL zlacpy(
'Lower', n-1, n-1, af( m-n+2, 1 ), lda,
250 CALL zungrq( n, n, min( m, n ), q, lda, taua, work, lwork, info )
254 CALL zlaset(
'Full', p, p, crogue, crogue, z, ldb )
256 $
CALL zlacpy(
'Lower', p-1, n, bf( 2, 1 ), ldb, z( 2, 1 ), ldb )
257 CALL zungqr( p, p, min( p, n ), z, ldb, taub, work, lwork, info )
261 CALL zlaset(
'Full', m, n, czero, czero, r, lda )
263 CALL zlacpy(
'Upper', m, m, af( 1, n-m+1 ), lda, r( 1, n-m+1 ),
266 CALL zlacpy(
'Full', m-n, n, af, lda, r, lda )
267 CALL zlacpy(
'Upper', n, n, af( m-n+1, 1 ), lda, r( m-n+1, 1 ),
273 CALL zlaset(
'Full', p, n, czero, czero, t, ldb )
274 CALL zlacpy(
'Upper', p, n, bf, ldb, t, ldb )
278 CALL zgemm(
'No transpose',
'Conjugate transpose', m, n, n, -cone,
279 $ a, lda, q, lda, cone, r, lda )
283 resid = zlange(
'1', m, n, r, lda, rwork )
284 IF( anorm.GT.zero )
THEN
285 result( 1 ) = ( ( resid / dble( max( 1, m, n ) ) ) / anorm ) /
293 CALL zgemm(
'Conjugate transpose',
'No transpose', p, n, p, cone,
294 $ z, ldb, b, ldb, czero, bwk, ldb )
295 CALL zgemm(
'No transpose',
'No transpose', p, n, n, cone, t, ldb,
296 $ q, lda, -cone, bwk, ldb )
300 resid = zlange(
'1', p, n, bwk, ldb, rwork )
301 IF( bnorm.GT.zero )
THEN
302 result( 2 ) = ( ( resid / dble( max( 1, p, m ) ) ) / bnorm ) /
310 CALL zlaset(
'Full', n, n, czero, cone, r, lda )
311 CALL zherk(
'Upper',
'No Transpose', n, n, -one, q, lda, one, r,
316 resid = zlanhe(
'1',
'Upper', n, r, lda, rwork )
317 result( 3 ) = ( resid / dble( max( 1, n ) ) ) / ulp
321 CALL zlaset(
'Full', p, p, czero, cone, t, ldb )
322 CALL zherk(
'Upper',
'Conjugate transpose', p, p, -one, z, ldb,
327 resid = zlanhe(
'1',
'Upper', p, t, ldb, rwork )
328 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 zggrqf(m, p, n, a, lda, taua, b, ldb, taub, work, lwork, info)
ZGGRQF
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 zgrqts(m, p, n, a, af, q, r, lda, taua, b, bf, z, t, bwk, ldb, taub, work, lwork, rwork, result)
ZGRQTS