85 DOUBLE PRECISION result(6)
91 COMPLEX*16,
ALLOCATABLE :: af(:,:), q(:,:),
92 $ l(:,:), rwork(:), work( : ), t(:,:),
93 $ cf(:,:), df(:,:), a(:,:), c(:,:), d(:,:)
98 parameter( zero = 0.0)
99 parameter( one = (1.0,0.0), czero=(0.0,0.0) )
102 INTEGER info, j, k, ll, lwork, ldt
103 DOUBLE PRECISION anorm, eps, resid, cnorm, dnorm
118 DATA iseed / 1988, 1989, 1990, 1991 /
123 lwork = max(2,ll)*max(2,ll)*nb
127 ALLOCATE ( a(m,n), af(m,n), q(n,n), l(ll,n), rwork(ll),
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 zgelqt( m, n, nb, af, m, t, ldt, work, info )
145 CALL zlaset(
'Full', n, n, czero, one, q, n )
146 CALL zgemlqt(
'R',
'N', n, n, k, nb, af, m, t, ldt, q, n,
151 CALL zlaset(
'Full', ll, n, czero, czero, l, ll )
152 CALL zlacpy(
'Lower', m, n, af, m, l, ll )
156 CALL zgemm(
'N',
'C', m, n, n, -one, a, m, q, n, one, l, ll )
157 anorm =
zlange(
'1', m, n, a, m, rwork )
158 resid =
zlange(
'1', m, n, l, ll, rwork )
159 IF( anorm.GT.zero )
THEN
160 result( 1 ) = resid / (eps*max(1,m)*anorm)
167 CALL zlaset(
'Full', n, n, czero, one, l, ll )
168 CALL zherk(
'U',
'C', n, n, dreal(-one), q, n, dreal(one), l, ll)
169 resid =
zlansy(
'1',
'Upper', n, l, ll, rwork )
170 result( 2 ) = resid / (eps*max(1,n))
175 CALL zlarnv( 2, iseed, n, d( 1, j ) )
177 dnorm =
zlange(
'1', n, m, d, n, rwork)
178 CALL zlacpy(
'Full', n, m, d, n, df, n )
182 CALL zgemlqt(
'L',
'N', n, m, k, nb, af, m, t, nb, df, n,
187 CALL zgemm(
'N',
'N', n, m, n, -one, q, n, d, n, one, df, n )
188 resid =
zlange(
'1', n, m, df, n, rwork )
189 IF( dnorm.GT.zero )
THEN
190 result( 3 ) = resid / (eps*max(1,m)*dnorm)
197 CALL zlacpy(
'Full', n, m, d, n, df, n )
201 CALL zgemlqt(
'L',
'C', n, m, k, nb, af, m, t, nb, df, n,
206 CALL zgemm(
'C',
'N', n, m, n, -one, q, n, d, n, one, df, n )
207 resid =
zlange(
'1', n, m, df, n, rwork )
208 IF( dnorm.GT.zero )
THEN
209 result( 4 ) = resid / (eps*max(1,m)*dnorm)
217 CALL zlarnv( 2, iseed, m, c( 1, j ) )
219 cnorm =
zlange(
'1', m, n, c, m, rwork)
220 CALL zlacpy(
'Full', m, n, c, m, cf, m )
224 CALL zgemlqt(
'R',
'N', m, n, k, nb, af, m, t, nb, cf, m,
229 CALL zgemm(
'N',
'N', m, n, n, -one, c, m, q, n, one, cf, m )
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', m, n, c, m, cf, m )
243 CALL zgemlqt(
'R',
'C', m, n, k, nb, af, m, t, nb, cf, m,
248 CALL zgemm(
'N',
'C', m, n, n, -one, c, m, q, n, one, cf, m )
249 resid =
zlange(
'1', m, n, cf, m, rwork )
250 IF( cnorm.GT.zero )
THEN
251 result( 6 ) = resid / (eps*max(1,m)*dnorm)
258 DEALLOCATE ( a, af, q, l, rwork, work, t, c, d, cf, df)
subroutine zgelqt(M, N, MB, A, LDA, T, LDT, WORK, INFO)
ZGELQT
double precision function dlamch(CMACH)
DLAMCH
subroutine zgemlqt(SIDE, TRANS, M, N, K, MB, V, LDV, T, LDT, C, LDC, WORK, INFO)
ZGEMLQT
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 zherk(UPLO, TRANS, N, K, ALPHA, A, LDA, BETA, C, LDC)
ZHERK
subroutine zlarnv(IDIST, ISEED, N, X)
ZLARNV returns a vector of random numbers from a uniform or normal distribution.