85 DOUBLE PRECISION result(6)
91 COMPLEX*16,
ALLOCATABLE :: af(:,:), q(:,:),
92 $ r(:,:), work( : ), t(:,:),
93 $ cf(:,:), df(:,:), a(:,:), c(:,:), d(:,:)
94 DOUBLE PRECISION,
ALLOCATABLE :: rwork(:)
99 parameter( zero = 0.0, one = (1.0,0.0), czero=(0.0,0.0) )
102 INTEGER info, j, k, l, lwork
103 DOUBLE PRECISION anorm, eps, resid, cnorm, dnorm
118 DATA iseed / 1988, 1989, 1990, 1991 /
123 lwork = max(2,l)*max(2,l)*nb
127 ALLOCATE ( a(m,n), af(m,n), q(m,m), r(m,l), rwork(l),
128 $ work(lwork), t(nb,n), c(m,n), cf(m,n),
135 CALL zlarnv( 2, iseed, m, a( 1, j ) )
137 CALL zlacpy(
'Full', m, n, a, m, af, m )
141 CALL zgeqrt( m, n, nb, af, m, t, ldt, work, info )
145 CALL zlaset(
'Full', m, m, czero, one, q, m )
146 CALL zgemqrt(
'R',
'N', m, m, k, nb, af, m, t, ldt, q, m,
151 CALL zlaset(
'Full', m, n, czero, czero, r, m )
152 CALL zlacpy(
'Upper', m, n, af, m, r, m )
156 CALL zgemm(
'C',
'N', m, n, m, -one, q, m, a, m, one, r, m )
157 anorm =
zlange(
'1', m, n, a, m, rwork )
158 resid =
zlange(
'1', m, n, r, m, rwork )
159 IF( anorm.GT.zero )
THEN
160 result( 1 ) = resid / (eps*max(1,m)*anorm)
167 CALL zlaset(
'Full', m, m, czero, one, r, m )
168 CALL zherk(
'U',
'C', m, m, dreal(-one), q, m, dreal(one), r, m )
169 resid =
zlansy(
'1',
'Upper', m, r, m, rwork )
170 result( 2 ) = resid / (eps*max(1,m))
175 CALL zlarnv( 2, iseed, m, c( 1, j ) )
177 cnorm =
zlange(
'1', m, n, c, m, rwork)
178 CALL zlacpy(
'Full', m, n, c, m, cf, m )
182 CALL zgemqrt(
'L',
'N', m, n, k, nb, af, m, t, nb, cf, m,
187 CALL zgemm(
'N',
'N', m, n, m, -one, q, m, c, m, one, cf, m )
188 resid =
zlange(
'1', m, n, cf, m, rwork )
189 IF( cnorm.GT.zero )
THEN
190 result( 3 ) = resid / (eps*max(1,m)*cnorm)
197 CALL zlacpy(
'Full', m, n, c, m, cf, m )
201 CALL zgemqrt(
'L',
'C', m, n, k, nb, af, m, t, nb, cf, m,
206 CALL zgemm(
'C',
'N', m, n, m, -one, q, m, c, m, one, cf, m )
207 resid =
zlange(
'1', m, n, cf, m, rwork )
208 IF( cnorm.GT.zero )
THEN
209 result( 4 ) = resid / (eps*max(1,m)*cnorm)
217 CALL zlarnv( 2, iseed, n, d( 1, j ) )
219 dnorm =
zlange(
'1', n, m, d, n, rwork)
220 CALL zlacpy(
'Full', n, m, d, n, df, n )
224 CALL zgemqrt(
'R',
'N', n, m, k, nb, af, m, t, nb, df, n,
229 CALL zgemm(
'N',
'N', n, m, m, -one, d, n, q, m, one, df, n )
230 resid =
zlange(
'1', n, m, df, n, rwork )
231 IF( cnorm.GT.zero )
THEN
232 result( 5 ) = resid / (eps*max(1,m)*dnorm)
239 CALL zlacpy(
'Full', n, m, d, n, df, n )
243 CALL zgemqrt(
'R',
'C', n, m, k, nb, af, m, t, nb, df, n,
248 CALL zgemm(
'N',
'C', n, m, m, -one, d, n, q, m, one, df, n )
249 resid =
zlange(
'1', n, m, df, n, rwork )
250 IF( cnorm.GT.zero )
THEN
251 result( 6 ) = resid / (eps*max(1,m)*dnorm)
258 DEALLOCATE ( a, af, q, r, rwork, work, t, c, d, cf, df)
double precision function dlamch(CMACH)
DLAMCH
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
subroutine zgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
ZGEMM
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 zlansy(NORM, UPLO, N, A, LDA, WORK)
ZLANSY 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 symmetric matrix.
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...
logical function lsame(CA, CB)
LSAME
subroutine zgeqrt(M, N, NB, A, LDA, T, LDT, WORK, INFO)
ZGEQRT
subroutine zherk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
ZHERK
subroutine zgemqrt(SIDE, TRANS, M, N, K, NB, V, LDV, T, LDT, C, LDC, WORK, INFO)
ZGEMQRT
subroutine zlarnv(IDIST, ISEED, N, X)
ZLARNV returns a vector of random numbers from a uniform or normal distribution.