82 IMPLICIT NONE
83
84
85
86
87
88
89 CHARACTER TSSW
90 INTEGER M, N, MB, NB
91
92 REAL RESULT(6)
93
94
95
96
97
98 COMPLEX, ALLOCATABLE :: AF(:,:), Q(:,:),
99 $ R(:,:), WORK( : ), T(:),
100 $ CF(:,:), DF(:,:), A(:,:), C(:,:), D(:,:), LQ(:,:)
101 REAL, ALLOCATABLE :: RWORK(:)
102
103
104 REAL ZERO
105 COMPLEX 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 REAL ANORM, EPS, RESID, CNORM, DNORM
112
113
114 INTEGER ISEED( 4 )
115 COMPLEX TQUERY( 5 ), WORKQUERY( 1 )
116
117
118 REAL SLAMCH, CLANGE, CLANSY
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 clarnv( 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 clarnv( 2, iseed, m/2, a( m/4, j ) )
163 END DO
164 END IF
165 END IF
166 CALL clacpy(
'Full', m, n, a, m, af, m )
167
168 IF (ts) THEN
169
170
171
172 CALL cgeqr( m, n, af, m, tquery, -1, workquery, -1, info )
173 tsize = int( tquery( 1 ) )
174 lwork = int( workquery( 1 ) )
175 CALL cgemqr(
'L',
'N', m, m, k, af, m, tquery, tsize, cf, m,
176 $ workquery, -1, info)
177 lwork = max( lwork, int( workquery( 1 ) ) )
178 CALL cgemqr(
'L',
'N', m, n, k, af, m, tquery, tsize, cf, m,
179 $ workquery, -1, info)
180 lwork = max( lwork, int( workquery( 1 ) ) )
181 CALL cgemqr(
'L',
'C', m, n, k, af, m, tquery, tsize, cf, m,
182 $ workquery, -1, info)
183 lwork = max( lwork, int( workquery( 1 ) ) )
184 CALL cgemqr(
'R',
'N', n, m, k, af, m, tquery, tsize, df, n,
185 $ workquery, -1, info)
186 lwork = max( lwork, int( workquery( 1 ) ) )
187 CALL cgemqr(
'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 = 'CGEQR'
193 CALL cgeqr( m, n, af, m, t, tsize, work, lwork, info )
194
195
196
197 CALL claset(
'Full', m, m, czero, one, q, m )
198 srnamt = 'CGEMQR'
199 CALL cgemqr(
'L',
'N', m, m, k, af, m, t, tsize, q, m,
200 $ work, lwork, info )
201
202
203
204 CALL claset(
'Full', m, n, czero, czero, r, m )
205 CALL clacpy(
'Upper', m, n, af, m, r, m )
206
207
208
209 CALL cgemm(
'C',
'N', m, n, m, -one, q, m, a, m, one, r, m )
210 anorm =
clange(
'1', m, n, a, m, rwork )
211 resid =
clange(
'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 claset(
'Full', m, m, czero, one, r, m )
221 CALL cherk(
'U',
'C', m, m, real(-one), q, m, real(one), r, m )
222 resid =
clansy(
'1',
'Upper', m, r, m, rwork )
223 result( 2 ) = resid / (eps*max(1,m))
224
225
226
227 DO j=1,n
228 CALL clarnv( 2, iseed, m, c( 1, j ) )
229 END DO
230 cnorm =
clange(
'1', m, n, c, m, rwork)
231 CALL clacpy(
'Full', m, n, c, m, cf, m )
232
233
234
235 srnamt = 'CGEMQR'
236 CALL cgemqr(
'L',
'N', m, n, k, af, m, t, tsize, cf, m,
237 $ work, lwork, info)
238
239
240
241 CALL cgemm(
'N',
'N', m, n, m, -one, q, m, c, m, one, cf, m )
242 resid =
clange(
'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 clacpy(
'Full', m, n, c, m, cf, m )
252
253
254
255 srnamt = 'CGEMQR'
256 CALL cgemqr(
'L',
'C', m, n, k, af, m, t, tsize, cf, m,
257 $ work, lwork, info)
258
259
260
261 CALL cgemm(
'C',
'N', m, n, m, -one, q, m, c, m, one, cf, m )
262 resid =
clange(
'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 clarnv( 2, iseed, n, d( 1, j ) )
273 END DO
274 dnorm =
clange(
'1', n, m, d, n, rwork)
275 CALL clacpy(
'Full', n, m, d, n, df, n )
276
277
278
279 srnamt = 'CGEMQR'
280 CALL cgemqr(
'R',
'N', n, m, k, af, m, t, tsize, df, n,
281 $ work, lwork, info)
282
283
284
285 CALL cgemm(
'N',
'N', n, m, m, -one, d, n, q, m, one, df, n )
286 resid =
clange(
'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 clacpy(
'Full', n, m, d, n, df, n )
296
297
298
299 CALL cgemqr(
'R',
'C', n, m, k, af, m, t, tsize, df, n,
300 $ work, lwork, info)
301
302
303
304 CALL cgemm(
'N',
'C', n, m, m, -one, d, n, q, m, one, df, n )
305 resid =
clange(
'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 cgelq( m, n, af, m, tquery, -1, workquery, -1, info )
316 tsize = int( tquery( 1 ) )
317 lwork = int( workquery( 1 ) )
318 CALL cgemlq(
'R',
'N', n, n, k, af, m, tquery, tsize, q, n,
319 $ workquery, -1, info )
320 lwork = max( lwork, int( workquery( 1 ) ) )
321 CALL cgemlq(
'L',
'N', n, m, k, af, m, tquery, tsize, df, n,
322 $ workquery, -1, info)
323 lwork = max( lwork, int( workquery( 1 ) ) )
324 CALL cgemlq(
'L',
'C', n, m, k, af, m, tquery, tsize, df, n,
325 $ workquery, -1, info)
326 lwork = max( lwork, int( workquery( 1 ) ) )
327 CALL cgemlq(
'R',
'N', m, n, k, af, m, tquery, tsize, cf, m,
328 $ workquery, -1, info)
329 lwork = max( lwork, int( workquery( 1 ) ) )
330 CALL cgemlq(
'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 = 'CGELQ'
336 CALL cgelq( m, n, af, m, t, tsize, work, lwork, info )
337
338
339
340
341 CALL claset(
'Full', n, n, czero, one, q, n )
342 srnamt = 'CGEMLQ'
343 CALL cgemlq(
'R',
'N', n, n, k, af, m, t, tsize, q, n,
344 $ work, lwork, info )
345
346
347
348 CALL claset(
'Full', m, n, czero, czero, lq, l )
349 CALL clacpy(
'Lower', m, n, af, m, lq, l )
350
351
352
353 CALL cgemm(
'N',
'C', m, n, n, -one, a, m, q, n, one, lq, l )
354 anorm =
clange(
'1', m, n, a, m, rwork )
355 resid =
clange(
'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 claset(
'Full', n, n, czero, one, lq, l )
365 CALL cherk(
'U',
'C', n, n, real(-one), q, n, real(one), lq, l)
366 resid =
clansy(
'1',
'Upper', n, lq, l, rwork )
367 result( 2 ) = resid / (eps*max(1,n))
368
369
370
371 DO j=1,m
372 CALL clarnv( 2, iseed, n, d( 1, j ) )
373 END DO
374 dnorm =
clange(
'1', n, m, d, n, rwork)
375 CALL clacpy(
'Full', n, m, d, n, df, n )
376
377
378
379 CALL cgemlq(
'L',
'N', n, m, k, af, m, t, tsize, df, n,
380 $ work, lwork, info)
381
382
383
384 CALL cgemm(
'N',
'N', n, m, n, -one, q, n, d, n, one, df, n )
385 resid =
clange(
'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 clacpy(
'Full', n, m, d, n, df, n )
395
396
397
398 CALL cgemlq(
'L',
'C', n, m, k, af, m, t, tsize, df, n,
399 $ work, lwork, info)
400
401
402
403 CALL cgemm(
'C',
'N', n, m, n, -one, q, n, d, n, one, df, n )
404 resid =
clange(
'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 clarnv( 2, iseed, m, c( 1, j ) )
415 END DO
416 cnorm =
clange(
'1', m, n, c, m, rwork)
417 CALL clacpy(
'Full', m, n, c, m, cf, m )
418
419
420
421 CALL cgemlq(
'R',
'N', m, n, k, af, m, t, tsize, cf, m,
422 $ work, lwork, info)
423
424
425
426 CALL cgemm(
'N',
'N', m, n, n, -one, c, m, q, n, one, cf, m )
427 resid =
clange(
'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 clacpy(
'Full', m, n, c, m, cf, m )
437
438
439
440 CALL cgemlq(
'R',
'C', m, n, k, af, m, t, tsize, cf, m,
441 $ work, lwork, info)
442
443
444
445 CALL cgemm(
'N',
'C', m, n, n, -one, c, m, q, n, one, cf, m )
446 resid =
clange(
'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 cgelq(m, n, a, lda, t, tsize, work, lwork, info)
CGELQ
subroutine cgemlq(side, trans, m, n, k, a, lda, t, tsize, c, ldc, work, lwork, info)
CGEMLQ
subroutine cgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
CGEMM
subroutine cgemqr(side, trans, m, n, k, a, lda, t, tsize, c, ldc, work, lwork, info)
CGEMQR
subroutine cgeqr(m, n, a, lda, t, tsize, work, lwork, info)
CGEQR
subroutine cherk(uplo, trans, n, k, alpha, a, lda, beta, c, ldc)
CHERK
integer function ilaenv(ispec, name, opts, n1, n2, n3, n4)
ILAENV
subroutine clacpy(uplo, m, n, a, lda, b, ldb)
CLACPY copies all or part of one two-dimensional array to another.
real function slamch(cmach)
SLAMCH
real function clange(norm, m, n, a, lda, work)
CLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
real function clansy(norm, uplo, n, a, lda, work)
CLANSY returns the value of the 1-norm, or the Frobenius norm, or the infinity norm,...
subroutine clarnv(idist, iseed, n, x)
CLARNV returns a vector of random numbers from a uniform or normal distribution.
subroutine claset(uplo, m, n, alpha, beta, a, lda)
CLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
logical function lsame(ca, cb)
LSAME