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
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 zggqrf(N, M, P, A, LDA, TAUA, B, LDB, TAUB, WORK, LWORK, INFO)
ZGGQRF
subroutine zgqrts(N, M, P, A, AF, Q, R, LDA, TAUA, B, BF, Z, T, BWK, LDB, TAUB, WORK, LWORK, RWORK, RESULT)
ZGQRTS
double precision function zlanhe(NORM, UPLO, N, A, LDA, WORK)
ZLANHE returns the value of the 1-norm, or the Frobenius norm, or the infinity norm, or the element of largest absolute value of a complex Hermitian matrix.
subroutine zungqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGQR
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
set ue cd $ADTTMP cat<< EOF > tmp f Program LinearEquations Implicit none Real b(3) integer i
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 zungrq(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGRQ
double precision function zlange(NORM, M, N, A, LDA, WORK)
ZLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
double precision function dlamch(CMACH)
DLAMCH