126 INTEGER M, N, MB1, NB1, NB2
128 DOUBLE PRECISION RESULT(6)
134 COMPLEX*16 ,
ALLOCATABLE :: A(:,:), AF(:,:), Q(:,:), R(:,:),
135 $ WORK( : ), T1(:,:), T2(:,:), DIAG(:),
136 $ C(:,:), CF(:,:), D(:,:), DF(:,:)
137 DOUBLE PRECISION,
ALLOCATABLE :: RWORK(:)
140 DOUBLE PRECISION ZERO
141 parameter( zero = 0.0d+0 )
142 COMPLEX*16 CONE, CZERO
143 parameter( cone = ( 1.0d+0, 0.0d+0 ),
144 $ czero = ( 0.0d+0, 0.0d+0 ) )
148 INTEGER INFO, I, J, K, L, LWORK, NB1_UB, NB2_UB, NRB
149 DOUBLE PRECISION ANORM, EPS, RESID, CNORM, DNORM
153 COMPLEX*16 WORKQUERY( 1 )
156 DOUBLE PRECISION DLAMCH, ZLANGE, ZLANSY
164 INTRINSIC ceiling, dble, max, min
167 CHARACTER(LEN=32) SRNAMT
170 COMMON / srmnamc / srnamt
173 DATA iseed / 1988, 1989, 1990, 1991 /
185 ALLOCATE ( a(m,n), af(m,n), q(l,l), r(m,l), rwork(l),
192 CALL zlarnv( 2, iseed, m, a( 1, j ) )
197 CALL zlarnv( 2, iseed, m/2, a( m/4, j ) )
201 CALL zlacpy(
'Full', m, n, a, m, af, m )
205 nrb = max( 1, ceiling( dble( m - n ) / dble( mb1 - n ) ) )
207 ALLOCATE ( t1( nb1, n * nrb ) )
208 ALLOCATE ( t2( nb2, n ) )
209 ALLOCATE ( diag( n ) )
215 nb1_ub = min( nb1, n)
219 nb2_ub = min( nb2, n)
221 CALL zlatsqr( m, n, mb1, nb1_ub, af, m, t1, nb1,
222 $ workquery, -1, info )
223 lwork = int( workquery( 1 ) )
224 CALL zungtsqr( m, n, mb1, nb1, af, m, t1, nb1, workquery, -1,
227 lwork = max( lwork, int( workquery( 1 ) ) )
232 lwork = max( lwork, nb2_ub * n, nb2_ub * m )
234 ALLOCATE ( work( lwork ) )
244 CALL zlatsqr( m, n, mb1, nb1_ub, af, m, t1, nb1, work, lwork,
250 CALL zlacpy(
'U', n, n, af, m, r, m )
255 CALL zungtsqr( m, n, mb1, nb1, af, m, t1, nb1, work, lwork,
262 CALL zunhr_col( m, n, nb2, af, m, t2, nb2, diag, info )
272 CALL zlacpy(
'U', n, n, r, m, af, m )
275 IF( diag( i ).EQ.-cone )
THEN
276 CALL zscal( n+1-i, -cone, af( i, i ), m )
285 CALL zlaset(
'Full', m, m, czero, cone, q, m )
288 CALL zgemqrt(
'L',
'N', m, m, k, nb2_ub, af, m, t2, nb2, q, m,
293 CALL zlaset(
'Full', m, n, czero, czero, r, m )
295 CALL zlacpy(
'Upper', m, n, af, m, r, m )
300 CALL zgemm(
'C',
'N', m, n, m, -cone, q, m, a, m, cone, r, m )
302 anorm =
zlange(
'1', m, n, a, m, rwork )
303 resid =
zlange(
'1', m, n, r, m, rwork )
304 IF( anorm.GT.zero )
THEN
305 result( 1 ) = resid / ( eps * max( 1, m ) * anorm )
313 CALL zlaset(
'Full', m, m, czero, cone, r, m )
314 CALL zherk(
'U',
'C', m, m, -cone, q, m, cone, r, m )
315 resid =
zlansy(
'1',
'Upper', m, r, m, rwork )
316 result( 2 ) = resid / ( eps * max( 1, m ) )
321 CALL zlarnv( 2, iseed, m, c( 1, j ) )
323 cnorm =
zlange(
'1', m, n, c, m, rwork )
324 CALL zlacpy(
'Full', m, n, c, m, cf, m )
329 CALL zgemqrt(
'L',
'N', m, n, k, nb2_ub, af, m, t2, nb2, cf, m,
335 CALL zgemm(
'N',
'N', m, n, m, -cone, q, m, c, m, cone, cf, m )
336 resid =
zlange(
'1', m, n, cf, m, rwork )
337 IF( cnorm.GT.zero )
THEN
338 result( 3 ) = resid / ( eps * max( 1, m ) * cnorm )
345 CALL zlacpy(
'Full', m, n, c, m, cf, m )
350 CALL zgemqrt(
'L',
'C', m, n, k, nb2_ub, af, m, t2, nb2, cf, m,
356 CALL zgemm(
'C',
'N', m, n, m, -cone, q, m, c, m, cone, cf, m )
357 resid =
zlange(
'1', m, n, cf, m, rwork )
358 IF( cnorm.GT.zero )
THEN
359 result( 4 ) = resid / ( eps * max( 1, m ) * cnorm )
367 CALL zlarnv( 2, iseed, n, d( 1, j ) )
369 dnorm =
zlange(
'1', n, m, d, n, rwork )
370 CALL zlacpy(
'Full', n, m, d, n, df, n )
375 CALL zgemqrt(
'R',
'N', n, m, k, nb2_ub, af, m, t2, nb2, df, n,
381 CALL zgemm(
'N',
'N', n, m, m, -cone, d, n, q, m, cone, df, n )
382 resid =
zlange(
'1', n, m, df, n, rwork )
383 IF( dnorm.GT.zero )
THEN
384 result( 5 ) = resid / ( eps * max( 1, m ) * dnorm )
391 CALL zlacpy(
'Full', n, m, d, n, df, n )
396 CALL zgemqrt(
'R',
'C', n, m, k, nb2_ub, af, m, t2, nb2, df, n,
402 CALL zgemm(
'N',
'C', n, m, m, -cone, d, n, q, m, cone, df, n )
403 resid =
zlange(
'1', n, m, df, n, rwork )
404 IF( dnorm.GT.zero )
THEN
405 result( 6 ) = resid / ( eps * max( 1, m ) * dnorm )
412 DEALLOCATE ( a, af, q, r, rwork, work, t1, t2, diag,
double precision function dlamch(CMACH)
DLAMCH
subroutine zscal(N, ZA, ZX, INCX)
ZSCAL
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
subroutine zherk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
ZHERK
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 ...
subroutine zgemqrt(SIDE, TRANS, M, N, K, NB, V, LDV, T, LDT, C, LDC, WORK, INFO)
ZGEMQRT
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
subroutine zlarnv(IDIST, ISEED, N, X)
ZLARNV returns a vector of random numbers from a uniform or normal distribution.
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 zungtsqr(M, N, MB, NB, A, LDA, T, LDT, WORK, LWORK, INFO)
ZUNGTSQR
subroutine zunhr_col(M, N, NB, A, LDA, T, LDT, D, INFO)
ZUNHR_COL
double precision function zlansy(NORM, UPLO, N, A, LDA, WORK)
ZLANSY returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
subroutine zlatsqr(M, N, MB, NB, A, LDA, T, LDT, WORK, LWORK, INFO)
ZLATSQR