73 IMPLICIT NONE
74
75
76
77
78
79
80 INTEGER M, N, NB, LDT
81
82 DOUBLE PRECISION RESULT(6)
83
84
85
86
87
88 COMPLEX*16, ALLOCATABLE :: AF(:,:), Q(:,:),
89 $ R(:,:), WORK( : ), T(:,:),
90 $ CF(:,:), DF(:,:), A(:,:), C(:,:), D(:,:)
91 DOUBLE PRECISION, ALLOCATABLE :: RWORK(:)
92
93
94 DOUBLE PRECISION ZERO
95 COMPLEX*16 ONE, CZERO
96 parameter( zero = 0.0, one = (1.0,0.0), czero=(0.0,0.0) )
97
98
99 INTEGER INFO, J, K, L, LWORK
100 DOUBLE PRECISION ANORM, EPS, RESID, CNORM, DNORM
101
102
103 INTEGER ISEED( 4 )
104
105
106 DOUBLE PRECISION DLAMCH
107 DOUBLE PRECISION ZLANGE, ZLANSY
108 LOGICAL LSAME
110
111
112 INTRINSIC max, min
113
114
115 DATA iseed / 1988, 1989, 1990, 1991 /
116
118 k = min(m,n)
119 l = max(m,n)
120 lwork = max(2,l)*max(2,l)*nb
121
122
123
124 ALLOCATE ( a(m,n), af(m,n), q(m,m), r(m,l), rwork(l),
125 $ work(lwork), t(nb,n), c(m,n), cf(m,n),
126 $ d(n,m), df(n,m) )
127
128
129
130 ldt=nb
131 DO j=1,n
132 CALL zlarnv( 2, iseed, m, a( 1, j ) )
133 END DO
134 CALL zlacpy(
'Full', m, n, a, m, af, m )
135
136
137
138 CALL zgeqrt( m, n, nb, af, m, t, ldt, work, info )
139
140
141
142 CALL zlaset(
'Full', m, m, czero, one, q, m )
143 CALL zgemqrt(
'R',
'N', m, m, k, nb, af, m, t, ldt, q, m,
144 $ work, info )
145
146
147
148 CALL zlaset(
'Full', m, n, czero, czero, r, m )
149 CALL zlacpy(
'Upper', m, n, af, m, r, m )
150
151
152
153 CALL zgemm(
'C',
'N', m, n, m, -one, q, m, a, m, one, r, m )
154 anorm =
zlange(
'1', m, n, a, m, rwork )
155 resid =
zlange(
'1', m, n, r, m, rwork )
156 IF( anorm.GT.zero ) THEN
157 result( 1 ) = resid / (eps*max(1,m)*anorm)
158 ELSE
159 result( 1 ) = zero
160 END IF
161
162
163
164 CALL zlaset(
'Full', m, m, czero, one, r, m )
165 CALL zherk(
'U',
'C', m, m, dreal(-one), q, m, dreal(one), r, m )
166 resid =
zlansy(
'1',
'Upper', m, r, m, rwork )
167 result( 2 ) = resid / (eps*max(1,m))
168
169
170
171 DO j=1,n
172 CALL zlarnv( 2, iseed, m, c( 1, j ) )
173 END DO
174 cnorm =
zlange(
'1', m, n, c, m, rwork)
175 CALL zlacpy(
'Full', m, n, c, m, cf, m )
176
177
178
179 CALL zgemqrt(
'L',
'N', m, n, k, nb, af, m, t, nb, cf, m,
180 $ work, info)
181
182
183
184 CALL zgemm(
'N',
'N', m, n, m, -one, q, m, c, m, one, cf, m )
185 resid =
zlange(
'1', m, n, cf, m, rwork )
186 IF( cnorm.GT.zero ) THEN
187 result( 3 ) = resid / (eps*max(1,m)*cnorm)
188 ELSE
189 result( 3 ) = zero
190 END IF
191
192
193
194 CALL zlacpy(
'Full', m, n, c, m, cf, m )
195
196
197
198 CALL zgemqrt(
'L',
'C', m, n, k, nb, af, m, t, nb, cf, m,
199 $ work, info)
200
201
202
203 CALL zgemm(
'C',
'N', m, n, m, -one, q, m, c, m, one, cf, m )
204 resid =
zlange(
'1', m, n, cf, m, rwork )
205 IF( cnorm.GT.zero ) THEN
206 result( 4 ) = resid / (eps*max(1,m)*cnorm)
207 ELSE
208 result( 4 ) = zero
209 END IF
210
211
212
213 DO j=1,m
214 CALL zlarnv( 2, iseed, n, d( 1, j ) )
215 END DO
216 dnorm =
zlange(
'1', n, m, d, n, rwork)
217 CALL zlacpy(
'Full', n, m, d, n, df, n )
218
219
220
221 CALL zgemqrt(
'R',
'N', n, m, k, nb, af, m, t, nb, df, n,
222 $ work, info)
223
224
225
226 CALL zgemm(
'N',
'N', n, m, m, -one, d, n, q, m, one, df, n )
227 resid =
zlange(
'1', n, m, df, n, rwork )
228 IF( cnorm.GT.zero ) THEN
229 result( 5 ) = resid / (eps*max(1,m)*dnorm)
230 ELSE
231 result( 5 ) = zero
232 END IF
233
234
235
236 CALL zlacpy(
'Full', n, m, d, n, df, n )
237
238
239
240 CALL zgemqrt(
'R',
'C', n, m, k, nb, af, m, t, nb, df, n,
241 $ work, info)
242
243
244
245 CALL zgemm(
'N',
'C', n, m, m, -one, d, n, q, m, one, df, n )
246 resid =
zlange(
'1', n, m, df, n, rwork )
247 IF( cnorm.GT.zero ) THEN
248 result( 6 ) = resid / (eps*max(1,m)*dnorm)
249 ELSE
250 result( 6 ) = zero
251 END IF
252
253
254
255 DEALLOCATE ( a, af, q, r, rwork, work, t, c, d, cf, df)
256
257 RETURN
subroutine zgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
ZGEMM
subroutine zgemqrt(side, trans, m, n, k, nb, v, ldv, t, ldt, c, ldc, work, info)
ZGEMQRT
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 zlacpy(uplo, m, n, a, lda, b, ldb)
ZLACPY copies all or part of one two-dimensional array to another.
double precision function dlamch(cmach)
DLAMCH
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,...
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.
logical function lsame(ca, cb)
LSAME