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