82 IMPLICIT NONE
83
84
85
86
87
88
89 CHARACTER TSSW
90 INTEGER M, N, MB, NB
91
92 DOUBLE PRECISION RESULT(6)
93
94
95
96
97
98 COMPLEX*16, ALLOCATABLE :: AF(:,:), Q(:,:),
99 $ R(:,:), WORK( : ), T(:),
100 $ CF(:,:), DF(:,:), A(:,:), C(:,:), D(:,:), LQ(:,:)
101 DOUBLE PRECISION, ALLOCATABLE :: RWORK(:)
102
103
104 DOUBLE PRECISION ZERO
105 COMPLEX*16 ONE, CZERO
106 parameter( zero = 0.0, one = (1.0,0.0), czero=(0.0,0.0) )
107
108
109 LOGICAL TESTZEROS, TS
110 INTEGER INFO, J, K, L, LWORK, TSIZE, MNB
111 DOUBLE PRECISION ANORM, EPS, RESID, CNORM, DNORM
112
113
114 INTEGER ISEED( 4 )
115 COMPLEX*16 TQUERY( 5 ), WORKQUERY( 1 )
116
117
118 DOUBLE PRECISION DLAMCH, ZLANGE, ZLANSY
119 LOGICAL LSAME
120 INTEGER ILAENV
122
123
124 INTRINSIC max, min
125
126 CHARACTER*32 srnamt
127
128
129 COMMON / srnamc / srnamt
130
131
132 DATA iseed / 1988, 1989, 1990, 1991 /
133
134
135
136 ts =
lsame(tssw,
'TS')
137
138
139
140 testzeros = .false.
141
143 k = min(m,n)
144 l = max(m,n,1)
145 mnb = max( mb, nb)
146 lwork = max(3,l)*mnb
147
148
149
150 ALLOCATE ( a(m,n), af(m,n), q(l,l), r(m,l), rwork(l),
151 $ c(m,n), cf(m,n),
152 $ d(n,m), df(n,m), lq(l,n) )
153
154
155
156 DO j=1,n
157 CALL zlarnv( 2, iseed, m, a( 1, j ) )
158 END DO
159 IF (testzeros) THEN
160 IF (m.GE.4) THEN
161 DO j=1,n
162 CALL zlarnv( 2, iseed, m/2, a( m/4, j ) )
163 END DO
164 END IF
165 END IF
166 CALL zlacpy(
'Full', m, n, a, m, af, m )
167
168 IF (ts) THEN
169
170
171
172 CALL zgeqr( m, n, af, m, tquery, -1, workquery, -1, info )
173 tsize = int( tquery( 1 ) )
174 lwork = int( workquery( 1 ) )
175 CALL zgemqr(
'L',
'N', m, m, k, af, m, tquery, tsize, cf, m,
176 $ workquery, -1, info)
177 lwork = max( lwork, int( workquery( 1 ) ) )
178 CALL zgemqr(
'L',
'N', m, n, k, af, m, tquery, tsize, cf, m,
179 $ workquery, -1, info)
180 lwork = max( lwork, int( workquery( 1 ) ) )
181 CALL zgemqr(
'L',
'C', m, n, k, af, m, tquery, tsize, cf, m,
182 $ workquery, -1, info)
183 lwork = max( lwork, int( workquery( 1 ) ) )
184 CALL zgemqr(
'R',
'N', n, m, k, af, m, tquery, tsize, df, n,
185 $ workquery, -1, info)
186 lwork = max( lwork, int( workquery( 1 ) ) )
187 CALL zgemqr(
'R',
'C', n, m, k, af, m, tquery, tsize, df, n,
188 $ workquery, -1, info)
189 lwork = max( lwork, int( workquery( 1 ) ) )
190 ALLOCATE ( t( tsize ) )
191 ALLOCATE ( work( lwork ) )
192 srnamt = 'ZGEQR'
193 CALL zgeqr( m, n, af, m, t, tsize, work, lwork, info )
194
195
196
197 CALL zlaset(
'Full', m, m, czero, one, q, m )
198 srnamt = 'ZGEMQR'
199 CALL zgemqr(
'L',
'N', m, m, k, af, m, t, tsize, q, m,
200 $ work, lwork, info )
201
202
203
204 CALL zlaset(
'Full', m, n, czero, czero, r, m )
205 CALL zlacpy(
'Upper', m, n, af, m, r, m )
206
207
208
209 CALL zgemm(
'C',
'N', m, n, m, -one, q, m, a, m, one, r, m )
210 anorm =
zlange(
'1', m, n, a, m, rwork )
211 resid =
zlange(
'1', m, n, r, m, rwork )
212 IF( anorm.GT.zero ) THEN
213 result( 1 ) = resid / (eps*max(1,m)*anorm)
214 ELSE
215 result( 1 ) = zero
216 END IF
217
218
219
220 CALL zlaset(
'Full', m, m, czero, one, r, m )
221 CALL zherk(
'U',
'C', m, m, dreal(-one), q, m, dreal(one), r, m )
222 resid =
zlansy(
'1',
'Upper', m, r, m, rwork )
223 result( 2 ) = resid / (eps*max(1,m))
224
225
226
227 DO j=1,n
228 CALL zlarnv( 2, iseed, m, c( 1, j ) )
229 END DO
230 cnorm =
zlange(
'1', m, n, c, m, rwork)
231 CALL zlacpy(
'Full', m, n, c, m, cf, m )
232
233
234
235 srnamt = 'ZGEMQR'
236 CALL zgemqr(
'L',
'N', m, n, k, af, m, t, tsize, cf, m,
237 $ work, lwork, info)
238
239
240
241 CALL zgemm(
'N',
'N', m, n, m, -one, q, m, c, m, one, cf, m )
242 resid =
zlange(
'1', m, n, cf, m, rwork )
243 IF( cnorm.GT.zero ) THEN
244 result( 3 ) = resid / (eps*max(1,m)*cnorm)
245 ELSE
246 result( 3 ) = zero
247 END IF
248
249
250
251 CALL zlacpy(
'Full', m, n, c, m, cf, m )
252
253
254
255 srnamt = 'ZGEMQR'
256 CALL zgemqr(
'L',
'C', m, n, k, af, m, t, tsize, cf, m,
257 $ work, lwork, info)
258
259
260
261 CALL zgemm(
'C',
'N', m, n, m, -one, q, m, c, m, one, cf, m )
262 resid =
zlange(
'1', m, n, cf, m, rwork )
263 IF( cnorm.GT.zero ) THEN
264 result( 4 ) = resid / (eps*max(1,m)*cnorm)
265 ELSE
266 result( 4 ) = zero
267 END IF
268
269
270
271 DO j=1,m
272 CALL zlarnv( 2, iseed, n, d( 1, j ) )
273 END DO
274 dnorm =
zlange(
'1', n, m, d, n, rwork)
275 CALL zlacpy(
'Full', n, m, d, n, df, n )
276
277
278
279 srnamt = 'ZGEMQR'
280 CALL zgemqr(
'R',
'N', n, m, k, af, m, t, tsize, df, n,
281 $ work, lwork, info)
282
283
284
285 CALL zgemm(
'N',
'N', n, m, m, -one, d, n, q, m, one, df, n )
286 resid =
zlange(
'1', n, m, df, n, rwork )
287 IF( dnorm.GT.zero ) THEN
288 result( 5 ) = resid / (eps*max(1,m)*dnorm)
289 ELSE
290 result( 5 ) = zero
291 END IF
292
293
294
295 CALL zlacpy(
'Full', n, m, d, n, df, n )
296
297
298
299 CALL zgemqr(
'R',
'C', n, m, k, af, m, t, tsize, df, n,
300 $ work, lwork, info)
301
302
303
304 CALL zgemm(
'N',
'C', n, m, m, -one, d, n, q, m, one, df, n )
305 resid =
zlange(
'1', n, m, df, n, rwork )
306 IF( cnorm.GT.zero ) THEN
307 result( 6 ) = resid / (eps*max(1,m)*dnorm)
308 ELSE
309 result( 6 ) = zero
310 END IF
311
312
313
314 ELSE
315 CALL zgelq( m, n, af, m, tquery, -1, workquery, -1, info )
316 tsize = int( tquery( 1 ) )
317 lwork = int( workquery( 1 ) )
318 CALL zgemlq(
'R',
'N', n, n, k, af, m, tquery, tsize, q, n,
319 $ workquery, -1, info )
320 lwork = max( lwork, int( workquery( 1 ) ) )
321 CALL zgemlq(
'L',
'N', n, m, k, af, m, tquery, tsize, df, n,
322 $ workquery, -1, info)
323 lwork = max( lwork, int( workquery( 1 ) ) )
324 CALL zgemlq(
'L',
'C', n, m, k, af, m, tquery, tsize, df, n,
325 $ workquery, -1, info)
326 lwork = max( lwork, int( workquery( 1 ) ) )
327 CALL zgemlq(
'R',
'N', m, n, k, af, m, tquery, tsize, cf, m,
328 $ workquery, -1, info)
329 lwork = max( lwork, int( workquery( 1 ) ) )
330 CALL zgemlq(
'R',
'C', m, n, k, af, m, tquery, tsize, cf, m,
331 $ workquery, -1, info)
332 lwork = max( lwork, int( workquery( 1 ) ) )
333 ALLOCATE ( t( tsize ) )
334 ALLOCATE ( work( lwork ) )
335 srnamt = 'ZGELQ'
336 CALL zgelq( m, n, af, m, t, tsize, work, lwork, info )
337
338
339
340
341 CALL zlaset(
'Full', n, n, czero, one, q, n )
342 srnamt = 'ZGEMLQ'
343 CALL zgemlq(
'R',
'N', n, n, k, af, m, t, tsize, q, n,
344 $ work, lwork, info )
345
346
347
348 CALL zlaset(
'Full', m, n, czero, czero, lq, l )
349 CALL zlacpy(
'Lower', m, n, af, m, lq, l )
350
351
352
353 CALL zgemm(
'N',
'C', m, n, n, -one, a, m, q, n, one, lq, l )
354 anorm =
zlange(
'1', m, n, a, m, rwork )
355 resid =
zlange(
'1', m, n, lq, l, rwork )
356 IF( anorm.GT.zero ) THEN
357 result( 1 ) = resid / (eps*max(1,n)*anorm)
358 ELSE
359 result( 1 ) = zero
360 END IF
361
362
363
364 CALL zlaset(
'Full', n, n, czero, one, lq, l )
365 CALL zherk(
'U',
'C', n, n, dreal(-one), q, n, dreal(one), lq, l)
366 resid =
zlansy(
'1',
'Upper', n, lq, l, rwork )
367 result( 2 ) = resid / (eps*max(1,n))
368
369
370
371 DO j=1,m
372 CALL zlarnv( 2, iseed, n, d( 1, j ) )
373 END DO
374 dnorm =
zlange(
'1', n, m, d, n, rwork)
375 CALL zlacpy(
'Full', n, m, d, n, df, n )
376
377
378
379 CALL zgemlq(
'L',
'N', n, m, k, af, m, t, tsize, df, n,
380 $ work, lwork, info)
381
382
383
384 CALL zgemm(
'N',
'N', n, m, n, -one, q, n, d, n, one, df, n )
385 resid =
zlange(
'1', n, m, df, n, rwork )
386 IF( dnorm.GT.zero ) THEN
387 result( 3 ) = resid / (eps*max(1,n)*dnorm)
388 ELSE
389 result( 3 ) = zero
390 END IF
391
392
393
394 CALL zlacpy(
'Full', n, m, d, n, df, n )
395
396
397
398 CALL zgemlq(
'L',
'C', n, m, k, af, m, t, tsize, df, n,
399 $ work, lwork, info)
400
401
402
403 CALL zgemm(
'C',
'N', n, m, n, -one, q, n, d, n, one, df, n )
404 resid =
zlange(
'1', n, m, df, n, rwork )
405 IF( dnorm.GT.zero ) THEN
406 result( 4 ) = resid / (eps*max(1,n)*dnorm)
407 ELSE
408 result( 4 ) = zero
409 END IF
410
411
412
413 DO j=1,n
414 CALL zlarnv( 2, iseed, m, c( 1, j ) )
415 END DO
416 cnorm =
zlange(
'1', m, n, c, m, rwork)
417 CALL zlacpy(
'Full', m, n, c, m, cf, m )
418
419
420
421 CALL zgemlq(
'R',
'N', m, n, k, af, m, t, tsize, cf, m,
422 $ work, lwork, info)
423
424
425
426 CALL zgemm(
'N',
'N', m, n, n, -one, c, m, q, n, one, cf, m )
427 resid =
zlange(
'1', n, m, df, n, rwork )
428 IF( cnorm.GT.zero ) THEN
429 result( 5 ) = resid / (eps*max(1,n)*cnorm)
430 ELSE
431 result( 5 ) = zero
432 END IF
433
434
435
436 CALL zlacpy(
'Full', m, n, c, m, cf, m )
437
438
439
440 CALL zgemlq(
'R',
'C', m, n, k, af, m, t, tsize, cf, m,
441 $ work, lwork, info)
442
443
444
445 CALL zgemm(
'N',
'C', m, n, n, -one, c, m, q, n, one, cf, m )
446 resid =
zlange(
'1', m, n, cf, m, rwork )
447 IF( cnorm.GT.zero ) THEN
448 result( 6 ) = resid / (eps*max(1,n)*cnorm)
449 ELSE
450 result( 6 ) = zero
451 END IF
452
453 END IF
454
455
456
457 DEALLOCATE ( a, af, q, r, rwork, work, t, c, d, cf, df)
458
459 RETURN
subroutine zgelq(m, n, a, lda, t, tsize, work, lwork, info)
ZGELQ
subroutine zgemlq(side, trans, m, n, k, a, lda, t, tsize, c, ldc, work, lwork, info)
ZGEMLQ
subroutine zgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
ZGEMM
subroutine zgemqr(side, trans, m, n, k, a, lda, t, tsize, c, ldc, work, lwork, info)
ZGEMQR
subroutine zgeqr(m, n, a, lda, t, tsize, work, lwork, info)
ZGEQR
subroutine zherk(uplo, trans, n, k, alpha, a, lda, beta, c, ldc)
ZHERK
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
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