1 SUBROUTINE ZCHKGG( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
2 $ TSTDIF, THRSHN, NOUNIT, A, LDA, B, H, T, S1,
3 $ S2, P1, P2, U, LDU, V, Q, Z, ALPHA1, BETA1,
4 $ ALPHA3, BETA3, EVECTL, EVECTR, WORK, LWORK,
5 $ RWORK, LLWORK, RESULT, INFO )
6 *
7 * -- LAPACK test routine (version 3.1) --
8 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
9 * November 2006
10 *
11 * .. Scalar Arguments ..
12 LOGICAL TSTDIF
13 INTEGER INFO, LDA, LDU, LWORK, NOUNIT, NSIZES, NTYPES
14 DOUBLE PRECISION THRESH, THRSHN
15 * ..
16 * .. Array Arguments ..
17 LOGICAL DOTYPE( * ), LLWORK( * )
18 INTEGER ISEED( 4 ), NN( * )
19 DOUBLE PRECISION RESULT( 15 ), RWORK( * )
20 COMPLEX*16 A( LDA, * ), ALPHA1( * ), ALPHA3( * ),
21 $ B( LDA, * ), BETA1( * ), BETA3( * ),
22 $ EVECTL( LDU, * ), EVECTR( LDU, * ),
23 $ H( LDA, * ), P1( LDA, * ), P2( LDA, * ),
24 $ Q( LDU, * ), S1( LDA, * ), S2( LDA, * ),
25 $ T( LDA, * ), U( LDU, * ), V( LDU, * ),
26 $ WORK( * ), Z( LDU, * )
27 * ..
28 *
29 * Purpose
30 * =======
31 *
32 * ZCHKGG checks the nonsymmetric generalized eigenvalue problem
33 * routines.
34 * H H H
35 * ZGGHRD factors A and B as U H V and U T V , where means conjugate
36 * transpose, H is hessenberg, T is triangular and U and V are unitary.
37 *
38 * H H
39 * ZHGEQZ factors H and T as Q S Z and Q P Z , where P and S are upper
40 * triangular and Q and Z are unitary. It also computes the generalized
41 * eigenvalues (alpha(1),beta(1)),...,(alpha(n),beta(n)), where
42 * alpha(j)=S(j,j) and beta(j)=P(j,j) -- thus, w(j) = alpha(j)/beta(j)
43 * is a root of the generalized eigenvalue problem
44 *
45 * det( A - w(j) B ) = 0
46 *
47 * and m(j) = beta(j)/alpha(j) is a root of the essentially equivalent
48 * problem
49 *
50 * det( m(j) A - B ) = 0
51 *
52 * ZTGEVC computes the matrix L of left eigenvectors and the matrix R
53 * of right eigenvectors for the matrix pair ( S, P ). In the
54 * description below, l and r are left and right eigenvectors
55 * corresponding to the generalized eigenvalues (alpha,beta).
56 *
57 * When ZCHKGG is called, a number of matrix "sizes" ("n's") and a
58 * number of matrix "types" are specified. For each size ("n")
59 * and each type of matrix, one matrix will be generated and used
60 * to test the nonsymmetric eigenroutines. For each matrix, 13
61 * tests will be performed. The first twelve "test ratios" should be
62 * small -- O(1). They will be compared with the threshhold THRESH:
63 *
64 * H
65 * (1) | A - U H V | / ( |A| n ulp )
66 *
67 * H
68 * (2) | B - U T V | / ( |B| n ulp )
69 *
70 * H
71 * (3) | I - UU | / ( n ulp )
72 *
73 * H
74 * (4) | I - VV | / ( n ulp )
75 *
76 * H
77 * (5) | H - Q S Z | / ( |H| n ulp )
78 *
79 * H
80 * (6) | T - Q P Z | / ( |T| n ulp )
81 *
82 * H
83 * (7) | I - QQ | / ( n ulp )
84 *
85 * H
86 * (8) | I - ZZ | / ( n ulp )
87 *
88 * (9) max over all left eigenvalue/-vector pairs (beta/alpha,l) of
89 * H
90 * | (beta A - alpha B) l | / ( ulp max( |beta A|, |alpha B| ) )
91 *
92 * (10) max over all left eigenvalue/-vector pairs (beta/alpha,l') of
93 * H
94 * | (beta H - alpha T) l' | / ( ulp max( |beta H|, |alpha T| ) )
95 *
96 * where the eigenvectors l' are the result of passing Q to
97 * DTGEVC and back transforming (JOB='B').
98 *
99 * (11) max over all right eigenvalue/-vector pairs (beta/alpha,r) of
100 *
101 * | (beta A - alpha B) r | / ( ulp max( |beta A|, |alpha B| ) )
102 *
103 * (12) max over all right eigenvalue/-vector pairs (beta/alpha,r') of
104 *
105 * | (beta H - alpha T) r' | / ( ulp max( |beta H|, |alpha T| ) )
106 *
107 * where the eigenvectors r' are the result of passing Z to
108 * DTGEVC and back transforming (JOB='B').
109 *
110 * The last three test ratios will usually be small, but there is no
111 * mathematical requirement that they be so. They are therefore
112 * compared with THRESH only if TSTDIF is .TRUE.
113 *
114 * (13) | S(Q,Z computed) - S(Q,Z not computed) | / ( |S| ulp )
115 *
116 * (14) | P(Q,Z computed) - P(Q,Z not computed) | / ( |P| ulp )
117 *
118 * (15) max( |alpha(Q,Z computed) - alpha(Q,Z not computed)|/|S| ,
119 * |beta(Q,Z computed) - beta(Q,Z not computed)|/|P| ) / ulp
120 *
121 * In addition, the normalization of L and R are checked, and compared
122 * with the threshhold THRSHN.
123 *
124 * Test Matrices
125 * ---- --------
126 *
127 * The sizes of the test matrices are specified by an array
128 * NN(1:NSIZES); the value of each element NN(j) specifies one size.
129 * The "types" are specified by a logical array DOTYPE( 1:NTYPES ); if
130 * DOTYPE(j) is .TRUE., then matrix type "j" will be generated.
131 * Currently, the list of possible types is:
132 *
133 * (1) ( 0, 0 ) (a pair of zero matrices)
134 *
135 * (2) ( I, 0 ) (an identity and a zero matrix)
136 *
137 * (3) ( 0, I ) (an identity and a zero matrix)
138 *
139 * (4) ( I, I ) (a pair of identity matrices)
140 *
141 * t t
142 * (5) ( J , J ) (a pair of transposed Jordan blocks)
143 *
144 * t ( I 0 )
145 * (6) ( X, Y ) where X = ( J 0 ) and Y = ( t )
146 * ( 0 I ) ( 0 J )
147 * and I is a k x k identity and J a (k+1)x(k+1)
148 * Jordan block; k=(N-1)/2
149 *
150 * (7) ( D, I ) where D is P*D1, P is a random unitary diagonal
151 * matrix (i.e., with random magnitude 1 entries
152 * on the diagonal), and D1=diag( 0, 1,..., N-1 )
153 * (i.e., a diagonal matrix with D1(1,1)=0,
154 * D1(2,2)=1, ..., D1(N,N)=N-1.)
155 * (8) ( I, D )
156 *
157 * (9) ( big*D, small*I ) where "big" is near overflow and small=1/big
158 *
159 * (10) ( small*D, big*I )
160 *
161 * (11) ( big*I, small*D )
162 *
163 * (12) ( small*I, big*D )
164 *
165 * (13) ( big*D, big*I )
166 *
167 * (14) ( small*D, small*I )
168 *
169 * (15) ( D1, D2 ) where D1=P*diag( 0, 0, 1, ..., N-3, 0 ) and
170 * D2=Q*diag( 0, N-3, N-4,..., 1, 0, 0 ), and
171 * P and Q are random unitary diagonal matrices.
172 * t t
173 * (16) U ( J , J ) V where U and V are random unitary matrices.
174 *
175 * (17) U ( T1, T2 ) V where T1 and T2 are upper triangular matrices
176 * with random O(1) entries above the diagonal
177 * and diagonal entries diag(T1) =
178 * P*( 0, 0, 1, ..., N-3, 0 ) and diag(T2) =
179 * Q*( 0, N-3, N-4,..., 1, 0, 0 )
180 *
181 * (18) U ( T1, T2 ) V diag(T1) = ( 0, 0, 1, 1, s, ..., s, 0 )
182 * diag(T2) = ( 0, 1, 0, 1,..., 1, 0 )
183 * s = machine precision.
184 *
185 * (19) U ( T1, T2 ) V diag(T1)=( 0,0,1,1, 1-d, ..., 1-(N-5)*d=s, 0 )
186 * diag(T2) = ( 0, 1, 0, 1, ..., 1, 0 )
187 *
188 * N-5
189 * (20) U ( T1, T2 ) V diag(T1)=( 0, 0, 1, 1, a, ..., a =s, 0 )
190 * diag(T2) = ( 0, 1, 0, 1, ..., 1, 0, 0 )
191 *
192 * (21) U ( T1, T2 ) V diag(T1)=( 0, 0, 1, r1, r2, ..., r(N-4), 0 )
193 * diag(T2) = ( 0, 1, 0, 1, ..., 1, 0, 0 )
194 * where r1,..., r(N-4) are random.
195 *
196 * (22) U ( big*T1, small*T2 ) V diag(T1) = P*( 0, 0, 1, ..., N-3, 0 )
197 * diag(T2) = ( 0, 1, ..., 1, 0, 0 )
198 *
199 * (23) U ( small*T1, big*T2 ) V diag(T1) = P*( 0, 0, 1, ..., N-3, 0 )
200 * diag(T2) = ( 0, 1, ..., 1, 0, 0 )
201 *
202 * (24) U ( small*T1, small*T2 ) V diag(T1) = P*( 0, 0, 1, ..., N-3, 0 )
203 * diag(T2) = ( 0, 1, ..., 1, 0, 0 )
204 *
205 * (25) U ( big*T1, big*T2 ) V diag(T1) = P*( 0, 0, 1, ..., N-3, 0 )
206 * diag(T2) = ( 0, 1, ..., 1, 0, 0 )
207 *
208 * (26) U ( T1, T2 ) V where T1 and T2 are random upper-triangular
209 * matrices.
210 *
211 * Arguments
212 * =========
213 *
214 * NSIZES (input) INTEGER
215 * The number of sizes of matrices to use. If it is zero,
216 * ZCHKGG does nothing. It must be at least zero.
217 *
218 * NN (input) INTEGER array, dimension (NSIZES)
219 * An array containing the sizes to be used for the matrices.
220 * Zero values will be skipped. The values must be at least
221 * zero.
222 *
223 * NTYPES (input) INTEGER
224 * The number of elements in DOTYPE. If it is zero, ZCHKGG
225 * does nothing. It must be at least zero. If it is MAXTYP+1
226 * and NSIZES is 1, then an additional type, MAXTYP+1 is
227 * defined, which is to use whatever matrix is in A. This
228 * is only useful if DOTYPE(1:MAXTYP) is .FALSE. and
229 * DOTYPE(MAXTYP+1) is .TRUE. .
230 *
231 * DOTYPE (input) LOGICAL array, dimension (NTYPES)
232 * If DOTYPE(j) is .TRUE., then for each size in NN a
233 * matrix of that size and of type j will be generated.
234 * If NTYPES is smaller than the maximum number of types
235 * defined (PARAMETER MAXTYP), then types NTYPES+1 through
236 * MAXTYP will not be generated. If NTYPES is larger
237 * than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES)
238 * will be ignored.
239 *
240 * ISEED (input/output) INTEGER array, dimension (4)
241 * On entry ISEED specifies the seed of the random number
242 * generator. The array elements should be between 0 and 4095;
243 * if not they will be reduced mod 4096. Also, ISEED(4) must
244 * be odd. The random number generator uses a linear
245 * congruential sequence limited to small integers, and so
246 * should produce machine independent random numbers. The
247 * values of ISEED are changed on exit, and can be used in the
248 * next call to ZCHKGG to continue the same random number
249 * sequence.
250 *
251 * THRESH (input) DOUBLE PRECISION
252 * A test will count as "failed" if the "error", computed as
253 * described above, exceeds THRESH. Note that the error
254 * is scaled to be O(1), so THRESH should be a reasonably
255 * small multiple of 1, e.g., 10 or 100. In particular,
256 * it should not depend on the precision (single vs. double)
257 * or the size of the matrix. It must be at least zero.
258 *
259 * TSTDIF (input) LOGICAL
260 * Specifies whether test ratios 13-15 will be computed and
261 * compared with THRESH.
262 * = .FALSE.: Only test ratios 1-12 will be computed and tested.
263 * Ratios 13-15 will be set to zero.
264 * = .TRUE.: All the test ratios 1-15 will be computed and
265 * tested.
266 *
267 * THRSHN (input) DOUBLE PRECISION
268 * Threshhold for reporting eigenvector normalization error.
269 * If the normalization of any eigenvector differs from 1 by
270 * more than THRSHN*ulp, then a special error message will be
271 * printed. (This is handled separately from the other tests,
272 * since only a compiler or programming error should cause an
273 * error message, at least if THRSHN is at least 5--10.)
274 *
275 * NOUNIT (input) INTEGER
276 * The FORTRAN unit number for printing out error messages
277 * (e.g., if a routine returns IINFO not equal to 0.)
278 *
279 * A (input/workspace) COMPLEX*16 array, dimension (LDA, max(NN))
280 * Used to hold the original A matrix. Used as input only
281 * if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and
282 * DOTYPE(MAXTYP+1)=.TRUE.
283 *
284 * LDA (input) INTEGER
285 * The leading dimension of A, B, H, T, S1, P1, S2, and P2.
286 * It must be at least 1 and at least max( NN ).
287 *
288 * B (input/workspace) COMPLEX*16 array, dimension (LDA, max(NN))
289 * Used to hold the original B matrix. Used as input only
290 * if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and
291 * DOTYPE(MAXTYP+1)=.TRUE.
292 *
293 * H (workspace) COMPLEX*16 array, dimension (LDA, max(NN))
294 * The upper Hessenberg matrix computed from A by ZGGHRD.
295 *
296 * T (workspace) COMPLEX*16 array, dimension (LDA, max(NN))
297 * The upper triangular matrix computed from B by ZGGHRD.
298 *
299 * S1 (workspace) COMPLEX*16 array, dimension (LDA, max(NN))
300 * The Schur (upper triangular) matrix computed from H by ZHGEQZ
301 * when Q and Z are also computed.
302 *
303 * S2 (workspace) COMPLEX*16 array, dimension (LDA, max(NN))
304 * The Schur (upper triangular) matrix computed from H by ZHGEQZ
305 * when Q and Z are not computed.
306 *
307 * P1 (workspace) COMPLEX*16 array, dimension (LDA, max(NN))
308 * The upper triangular matrix computed from T by ZHGEQZ
309 * when Q and Z are also computed.
310 *
311 * P2 (workspace) COMPLEX*16 array, dimension (LDA, max(NN))
312 * The upper triangular matrix computed from T by ZHGEQZ
313 * when Q and Z are not computed.
314 *
315 * U (workspace) COMPLEX*16 array, dimension (LDU, max(NN))
316 * The (left) unitary matrix computed by ZGGHRD.
317 *
318 * LDU (input) INTEGER
319 * The leading dimension of U, V, Q, Z, EVECTL, and EVEZTR. It
320 * must be at least 1 and at least max( NN ).
321 *
322 * V (workspace) COMPLEX*16 array, dimension (LDU, max(NN))
323 * The (right) unitary matrix computed by ZGGHRD.
324 *
325 * Q (workspace) COMPLEX*16 array, dimension (LDU, max(NN))
326 * The (left) unitary matrix computed by ZHGEQZ.
327 *
328 * Z (workspace) COMPLEX*16 array, dimension (LDU, max(NN))
329 * The (left) unitary matrix computed by ZHGEQZ.
330 *
331 * ALPHA1 (workspace) COMPLEX*16 array, dimension (max(NN))
332 * BETA1 (workspace) COMPLEX*16 array, dimension (max(NN))
333 * The generalized eigenvalues of (A,B) computed by ZHGEQZ
334 * when Q, Z, and the full Schur matrices are computed.
335 *
336 * ALPHA3 (workspace) COMPLEX*16 array, dimension (max(NN))
337 * BETA3 (workspace) COMPLEX*16 array, dimension (max(NN))
338 * The generalized eigenvalues of (A,B) computed by ZHGEQZ
339 * when neither Q, Z, nor the Schur matrices are computed.
340 *
341 * EVECTL (workspace) COMPLEX*16 array, dimension (LDU, max(NN))
342 * The (lower triangular) left eigenvector matrix for the
343 * matrices in S1 and P1.
344 *
345 * EVEZTR (workspace) COMPLEX*16 array, dimension (LDU, max(NN))
346 * The (upper triangular) right eigenvector matrix for the
347 * matrices in S1 and P1.
348 *
349 * WORK (workspace) COMPLEX*16 array, dimension (LWORK)
350 *
351 * LWORK (input) INTEGER
352 * The number of entries in WORK. This must be at least
353 * max( 4*N, 2 * N**2, 1 ), for all N=NN(j).
354 *
355 * RWORK (workspace) DOUBLE PRECISION array, dimension (2*max(NN))
356 *
357 * LLWORK (workspace) LOGICAL array, dimension (max(NN))
358 *
359 * RESULT (output) DOUBLE PRECISION array, dimension (15)
360 * The values computed by the tests described above.
361 * The values are currently limited to 1/ulp, to avoid
362 * overflow.
363 *
364 * INFO (output) INTEGER
365 * = 0: successful exit.
366 * < 0: if INFO = -i, the i-th argument had an illegal value.
367 * > 0: A routine returned an error code. INFO is the
368 * absolute value of the INFO value returned.
369 *
370 * =====================================================================
371 *
372 * .. Parameters ..
373 DOUBLE PRECISION ZERO, ONE
374 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
375 COMPLEX*16 CZERO, CONE
376 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
377 $ CONE = ( 1.0D+0, 0.0D+0 ) )
378 INTEGER MAXTYP
379 PARAMETER ( MAXTYP = 26 )
380 * ..
381 * .. Local Scalars ..
382 LOGICAL BADNN
383 INTEGER I1, IADD, IINFO, IN, J, JC, JR, JSIZE, JTYPE,
384 $ LWKOPT, MTYPES, N, N1, NERRS, NMATS, NMAX,
385 $ NTEST, NTESTT
386 DOUBLE PRECISION ANORM, BNORM, SAFMAX, SAFMIN, TEMP1, TEMP2,
387 $ ULP, ULPINV
388 COMPLEX*16 CTEMP
389 * ..
390 * .. Local Arrays ..
391 LOGICAL LASIGN( MAXTYP ), LBSIGN( MAXTYP )
392 INTEGER IOLDSD( 4 ), KADD( 6 ), KAMAGN( MAXTYP ),
393 $ KATYPE( MAXTYP ), KAZERO( MAXTYP ),
394 $ KBMAGN( MAXTYP ), KBTYPE( MAXTYP ),
395 $ KBZERO( MAXTYP ), KCLASS( MAXTYP ),
396 $ KTRIAN( MAXTYP ), KZ1( 6 ), KZ2( 6 )
397 DOUBLE PRECISION DUMMA( 4 ), RMAGN( 0: 3 )
398 COMPLEX*16 CDUMMA( 4 )
399 * ..
400 * .. External Functions ..
401 DOUBLE PRECISION DLAMCH, ZLANGE
402 COMPLEX*16 ZLARND
403 EXTERNAL DLAMCH, ZLANGE, ZLARND
404 * ..
405 * .. External Subroutines ..
406 EXTERNAL DLABAD, DLASUM, XERBLA, ZGEQR2, ZGET51, ZGET52,
407 $ ZGGHRD, ZHGEQZ, ZLACPY, ZLARFG, ZLASET, ZLATM4,
408 $ ZTGEVC, ZUNM2R
409 * ..
410 * .. Intrinsic Functions ..
411 INTRINSIC ABS, DBLE, DCONJG, MAX, MIN, SIGN
412 * ..
413 * .. Data statements ..
414 DATA KCLASS / 15*1, 10*2, 1*3 /
415 DATA KZ1 / 0, 1, 2, 1, 3, 3 /
416 DATA KZ2 / 0, 0, 1, 2, 1, 1 /
417 DATA KADD / 0, 0, 0, 0, 3, 2 /
418 DATA KATYPE / 0, 1, 0, 1, 2, 3, 4, 1, 4, 4, 1, 1, 4,
419 $ 4, 4, 2, 4, 5, 8, 7, 9, 4*4, 0 /
420 DATA KBTYPE / 0, 0, 1, 1, 2, -3, 1, 4, 1, 1, 4, 4,
421 $ 1, 1, -4, 2, -4, 8*8, 0 /
422 DATA KAZERO / 6*1, 2, 1, 2*2, 2*1, 2*2, 3, 1, 3,
423 $ 4*5, 4*3, 1 /
424 DATA KBZERO / 6*1, 1, 2, 2*1, 2*2, 2*1, 4, 1, 4,
425 $ 4*6, 4*4, 1 /
426 DATA KAMAGN / 8*1, 2, 3, 2, 3, 2, 3, 7*1, 2, 3, 3,
427 $ 2, 1 /
428 DATA KBMAGN / 8*1, 3, 2, 3, 2, 2, 3, 7*1, 3, 2, 3,
429 $ 2, 1 /
430 DATA KTRIAN / 16*0, 10*1 /
431 DATA LASIGN / 6*.FALSE., .TRUE., .FALSE., 2*.TRUE.,
432 $ 2*.FALSE., 3*.TRUE., .FALSE., .TRUE.,
433 $ 3*.FALSE., 5*.TRUE., .FALSE. /
434 DATA LBSIGN / 7*.FALSE., .TRUE., 2*.FALSE.,
435 $ 2*.TRUE., 2*.FALSE., .TRUE., .FALSE., .TRUE.,
436 $ 9*.FALSE. /
437 * ..
438 * .. Executable Statements ..
439 *
440 * Check for errors
441 *
442 INFO = 0
443 *
444 BADNN = .FALSE.
445 NMAX = 1
446 DO 10 J = 1, NSIZES
447 NMAX = MAX( NMAX, NN( J ) )
448 IF( NN( J ).LT.0 )
449 $ BADNN = .TRUE.
450 10 CONTINUE
451 *
452 LWKOPT = MAX( 2*NMAX*NMAX, 4*NMAX, 1 )
453 *
454 * Check for errors
455 *
456 IF( NSIZES.LT.0 ) THEN
457 INFO = -1
458 ELSE IF( BADNN ) THEN
459 INFO = -2
460 ELSE IF( NTYPES.LT.0 ) THEN
461 INFO = -3
462 ELSE IF( THRESH.LT.ZERO ) THEN
463 INFO = -6
464 ELSE IF( LDA.LE.1 .OR. LDA.LT.NMAX ) THEN
465 INFO = -10
466 ELSE IF( LDU.LE.1 .OR. LDU.LT.NMAX ) THEN
467 INFO = -19
468 ELSE IF( LWKOPT.GT.LWORK ) THEN
469 INFO = -30
470 END IF
471 *
472 IF( INFO.NE.0 ) THEN
473 CALL XERBLA( 'ZCHKGG', -INFO )
474 RETURN
475 END IF
476 *
477 * Quick return if possible
478 *
479 IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 )
480 $ RETURN
481 *
482 SAFMIN = DLAMCH( 'Safe minimum' )
483 ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' )
484 SAFMIN = SAFMIN / ULP
485 SAFMAX = ONE / SAFMIN
486 CALL DLABAD( SAFMIN, SAFMAX )
487 ULPINV = ONE / ULP
488 *
489 * The values RMAGN(2:3) depend on N, see below.
490 *
491 RMAGN( 0 ) = ZERO
492 RMAGN( 1 ) = ONE
493 *
494 * Loop over sizes, types
495 *
496 NTESTT = 0
497 NERRS = 0
498 NMATS = 0
499 *
500 DO 240 JSIZE = 1, NSIZES
501 N = NN( JSIZE )
502 N1 = MAX( 1, N )
503 RMAGN( 2 ) = SAFMAX*ULP / DBLE( N1 )
504 RMAGN( 3 ) = SAFMIN*ULPINV*N1
505 *
506 IF( NSIZES.NE.1 ) THEN
507 MTYPES = MIN( MAXTYP, NTYPES )
508 ELSE
509 MTYPES = MIN( MAXTYP+1, NTYPES )
510 END IF
511 *
512 DO 230 JTYPE = 1, MTYPES
513 IF( .NOT.DOTYPE( JTYPE ) )
514 $ GO TO 230
515 NMATS = NMATS + 1
516 NTEST = 0
517 *
518 * Save ISEED in case of an error.
519 *
520 DO 20 J = 1, 4
521 IOLDSD( J ) = ISEED( J )
522 20 CONTINUE
523 *
524 * Initialize RESULT
525 *
526 DO 30 J = 1, 15
527 RESULT( J ) = ZERO
528 30 CONTINUE
529 *
530 * Compute A and B
531 *
532 * Description of control parameters:
533 *
534 * KZLASS: =1 means w/o rotation, =2 means w/ rotation,
535 * =3 means random.
536 * KATYPE: the "type" to be passed to ZLATM4 for computing A.
537 * KAZERO: the pattern of zeros on the diagonal for A:
538 * =1: ( xxx ), =2: (0, xxx ) =3: ( 0, 0, xxx, 0 ),
539 * =4: ( 0, xxx, 0, 0 ), =5: ( 0, 0, 1, xxx, 0 ),
540 * =6: ( 0, 1, 0, xxx, 0 ). (xxx means a string of
541 * non-zero entries.)
542 * KAMAGN: the magnitude of the matrix: =0: zero, =1: O(1),
543 * =2: large, =3: small.
544 * LASIGN: .TRUE. if the diagonal elements of A are to be
545 * multiplied by a random magnitude 1 number.
546 * KBTYPE, KBZERO, KBMAGN, LBSIGN: the same, but for B.
547 * KTRIAN: =0: don't fill in the upper triangle, =1: do.
548 * KZ1, KZ2, KADD: used to implement KAZERO and KBZERO.
549 * RMAGN: used to implement KAMAGN and KBMAGN.
550 *
551 IF( MTYPES.GT.MAXTYP )
552 $ GO TO 110
553 IINFO = 0
554 IF( KCLASS( JTYPE ).LT.3 ) THEN
555 *
556 * Generate A (w/o rotation)
557 *
558 IF( ABS( KATYPE( JTYPE ) ).EQ.3 ) THEN
559 IN = 2*( ( N-1 ) / 2 ) + 1
560 IF( IN.NE.N )
561 $ CALL ZLASET( 'Full', N, N, CZERO, CZERO, A, LDA )
562 ELSE
563 IN = N
564 END IF
565 CALL ZLATM4( KATYPE( JTYPE ), IN, KZ1( KAZERO( JTYPE ) ),
566 $ KZ2( KAZERO( JTYPE ) ), LASIGN( JTYPE ),
567 $ RMAGN( KAMAGN( JTYPE ) ), ULP,
568 $ RMAGN( KTRIAN( JTYPE )*KAMAGN( JTYPE ) ), 4,
569 $ ISEED, A, LDA )
570 IADD = KADD( KAZERO( JTYPE ) )
571 IF( IADD.GT.0 .AND. IADD.LE.N )
572 $ A( IADD, IADD ) = RMAGN( KAMAGN( JTYPE ) )
573 *
574 * Generate B (w/o rotation)
575 *
576 IF( ABS( KBTYPE( JTYPE ) ).EQ.3 ) THEN
577 IN = 2*( ( N-1 ) / 2 ) + 1
578 IF( IN.NE.N )
579 $ CALL ZLASET( 'Full', N, N, CZERO, CZERO, B, LDA )
580 ELSE
581 IN = N
582 END IF
583 CALL ZLATM4( KBTYPE( JTYPE ), IN, KZ1( KBZERO( JTYPE ) ),
584 $ KZ2( KBZERO( JTYPE ) ), LBSIGN( JTYPE ),
585 $ RMAGN( KBMAGN( JTYPE ) ), ONE,
586 $ RMAGN( KTRIAN( JTYPE )*KBMAGN( JTYPE ) ), 4,
587 $ ISEED, B, LDA )
588 IADD = KADD( KBZERO( JTYPE ) )
589 IF( IADD.NE.0 )
590 $ B( IADD, IADD ) = RMAGN( KBMAGN( JTYPE ) )
591 *
592 IF( KCLASS( JTYPE ).EQ.2 .AND. N.GT.0 ) THEN
593 *
594 * Include rotations
595 *
596 * Generate U, V as Householder transformations times a
597 * diagonal matrix. (Note that ZLARFG makes U(j,j) and
598 * V(j,j) real.)
599 *
600 DO 50 JC = 1, N - 1
601 DO 40 JR = JC, N
602 U( JR, JC ) = ZLARND( 3, ISEED )
603 V( JR, JC ) = ZLARND( 3, ISEED )
604 40 CONTINUE
605 CALL ZLARFG( N+1-JC, U( JC, JC ), U( JC+1, JC ), 1,
606 $ WORK( JC ) )
607 WORK( 2*N+JC ) = SIGN( ONE, DBLE( U( JC, JC ) ) )
608 U( JC, JC ) = CONE
609 CALL ZLARFG( N+1-JC, V( JC, JC ), V( JC+1, JC ), 1,
610 $ WORK( N+JC ) )
611 WORK( 3*N+JC ) = SIGN( ONE, DBLE( V( JC, JC ) ) )
612 V( JC, JC ) = CONE
613 50 CONTINUE
614 CTEMP = ZLARND( 3, ISEED )
615 U( N, N ) = CONE
616 WORK( N ) = CZERO
617 WORK( 3*N ) = CTEMP / ABS( CTEMP )
618 CTEMP = ZLARND( 3, ISEED )
619 V( N, N ) = CONE
620 WORK( 2*N ) = CZERO
621 WORK( 4*N ) = CTEMP / ABS( CTEMP )
622 *
623 * Apply the diagonal matrices
624 *
625 DO 70 JC = 1, N
626 DO 60 JR = 1, N
627 A( JR, JC ) = WORK( 2*N+JR )*
628 $ DCONJG( WORK( 3*N+JC ) )*
629 $ A( JR, JC )
630 B( JR, JC ) = WORK( 2*N+JR )*
631 $ DCONJG( WORK( 3*N+JC ) )*
632 $ B( JR, JC )
633 60 CONTINUE
634 70 CONTINUE
635 CALL ZUNM2R( 'L', 'N', N, N, N-1, U, LDU, WORK, A,
636 $ LDA, WORK( 2*N+1 ), IINFO )
637 IF( IINFO.NE.0 )
638 $ GO TO 100
639 CALL ZUNM2R( 'R', 'C', N, N, N-1, V, LDU, WORK( N+1 ),
640 $ A, LDA, WORK( 2*N+1 ), IINFO )
641 IF( IINFO.NE.0 )
642 $ GO TO 100
643 CALL ZUNM2R( 'L', 'N', N, N, N-1, U, LDU, WORK, B,
644 $ LDA, WORK( 2*N+1 ), IINFO )
645 IF( IINFO.NE.0 )
646 $ GO TO 100
647 CALL ZUNM2R( 'R', 'C', N, N, N-1, V, LDU, WORK( N+1 ),
648 $ B, LDA, WORK( 2*N+1 ), IINFO )
649 IF( IINFO.NE.0 )
650 $ GO TO 100
651 END IF
652 ELSE
653 *
654 * Random matrices
655 *
656 DO 90 JC = 1, N
657 DO 80 JR = 1, N
658 A( JR, JC ) = RMAGN( KAMAGN( JTYPE ) )*
659 $ ZLARND( 4, ISEED )
660 B( JR, JC ) = RMAGN( KBMAGN( JTYPE ) )*
661 $ ZLARND( 4, ISEED )
662 80 CONTINUE
663 90 CONTINUE
664 END IF
665 *
666 ANORM = ZLANGE( '1', N, N, A, LDA, RWORK )
667 BNORM = ZLANGE( '1', N, N, B, LDA, RWORK )
668 *
669 100 CONTINUE
670 *
671 IF( IINFO.NE.0 ) THEN
672 WRITE( NOUNIT, FMT = 9999 )'Generator', IINFO, N, JTYPE,
673 $ IOLDSD
674 INFO = ABS( IINFO )
675 RETURN
676 END IF
677 *
678 110 CONTINUE
679 *
680 * Call ZGEQR2, ZUNM2R, and ZGGHRD to compute H, T, U, and V
681 *
682 CALL ZLACPY( ' ', N, N, A, LDA, H, LDA )
683 CALL ZLACPY( ' ', N, N, B, LDA, T, LDA )
684 NTEST = 1
685 RESULT( 1 ) = ULPINV
686 *
687 CALL ZGEQR2( N, N, T, LDA, WORK, WORK( N+1 ), IINFO )
688 IF( IINFO.NE.0 ) THEN
689 WRITE( NOUNIT, FMT = 9999 )'ZGEQR2', IINFO, N, JTYPE,
690 $ IOLDSD
691 INFO = ABS( IINFO )
692 GO TO 210
693 END IF
694 *
695 CALL ZUNM2R( 'L', 'C', N, N, N, T, LDA, WORK, H, LDA,
696 $ WORK( N+1 ), IINFO )
697 IF( IINFO.NE.0 ) THEN
698 WRITE( NOUNIT, FMT = 9999 )'ZUNM2R', IINFO, N, JTYPE,
699 $ IOLDSD
700 INFO = ABS( IINFO )
701 GO TO 210
702 END IF
703 *
704 CALL ZLASET( 'Full', N, N, CZERO, CONE, U, LDU )
705 CALL ZUNM2R( 'R', 'N', N, N, N, T, LDA, WORK, U, LDU,
706 $ WORK( N+1 ), IINFO )
707 IF( IINFO.NE.0 ) THEN
708 WRITE( NOUNIT, FMT = 9999 )'ZUNM2R', IINFO, N, JTYPE,
709 $ IOLDSD
710 INFO = ABS( IINFO )
711 GO TO 210
712 END IF
713 *
714 CALL ZGGHRD( 'V', 'I', N, 1, N, H, LDA, T, LDA, U, LDU, V,
715 $ LDU, IINFO )
716 IF( IINFO.NE.0 ) THEN
717 WRITE( NOUNIT, FMT = 9999 )'ZGGHRD', IINFO, N, JTYPE,
718 $ IOLDSD
719 INFO = ABS( IINFO )
720 GO TO 210
721 END IF
722 NTEST = 4
723 *
724 * Do tests 1--4
725 *
726 CALL ZGET51( 1, N, A, LDA, H, LDA, U, LDU, V, LDU, WORK,
727 $ RWORK, RESULT( 1 ) )
728 CALL ZGET51( 1, N, B, LDA, T, LDA, U, LDU, V, LDU, WORK,
729 $ RWORK, RESULT( 2 ) )
730 CALL ZGET51( 3, N, B, LDA, T, LDA, U, LDU, U, LDU, WORK,
731 $ RWORK, RESULT( 3 ) )
732 CALL ZGET51( 3, N, B, LDA, T, LDA, V, LDU, V, LDU, WORK,
733 $ RWORK, RESULT( 4 ) )
734 *
735 * Call ZHGEQZ to compute S1, P1, S2, P2, Q, and Z, do tests.
736 *
737 * Compute T1 and UZ
738 *
739 * Eigenvalues only
740 *
741 CALL ZLACPY( ' ', N, N, H, LDA, S2, LDA )
742 CALL ZLACPY( ' ', N, N, T, LDA, P2, LDA )
743 NTEST = 5
744 RESULT( 5 ) = ULPINV
745 *
746 CALL ZHGEQZ( 'E', 'N', 'N', N, 1, N, S2, LDA, P2, LDA,
747 $ ALPHA3, BETA3, Q, LDU, Z, LDU, WORK, LWORK,
748 $ RWORK, IINFO )
749 IF( IINFO.NE.0 ) THEN
750 WRITE( NOUNIT, FMT = 9999 )'ZHGEQZ(E)', IINFO, N, JTYPE,
751 $ IOLDSD
752 INFO = ABS( IINFO )
753 GO TO 210
754 END IF
755 *
756 * Eigenvalues and Full Schur Form
757 *
758 CALL ZLACPY( ' ', N, N, H, LDA, S2, LDA )
759 CALL ZLACPY( ' ', N, N, T, LDA, P2, LDA )
760 *
761 CALL ZHGEQZ( 'S', 'N', 'N', N, 1, N, S2, LDA, P2, LDA,
762 $ ALPHA1, BETA1, Q, LDU, Z, LDU, WORK, LWORK,
763 $ RWORK, IINFO )
764 IF( IINFO.NE.0 ) THEN
765 WRITE( NOUNIT, FMT = 9999 )'ZHGEQZ(S)', IINFO, N, JTYPE,
766 $ IOLDSD
767 INFO = ABS( IINFO )
768 GO TO 210
769 END IF
770 *
771 * Eigenvalues, Schur Form, and Schur Vectors
772 *
773 CALL ZLACPY( ' ', N, N, H, LDA, S1, LDA )
774 CALL ZLACPY( ' ', N, N, T, LDA, P1, LDA )
775 *
776 CALL ZHGEQZ( 'S', 'I', 'I', N, 1, N, S1, LDA, P1, LDA,
777 $ ALPHA1, BETA1, Q, LDU, Z, LDU, WORK, LWORK,
778 $ RWORK, IINFO )
779 IF( IINFO.NE.0 ) THEN
780 WRITE( NOUNIT, FMT = 9999 )'ZHGEQZ(V)', IINFO, N, JTYPE,
781 $ IOLDSD
782 INFO = ABS( IINFO )
783 GO TO 210
784 END IF
785 *
786 NTEST = 8
787 *
788 * Do Tests 5--8
789 *
790 CALL ZGET51( 1, N, H, LDA, S1, LDA, Q, LDU, Z, LDU, WORK,
791 $ RWORK, RESULT( 5 ) )
792 CALL ZGET51( 1, N, T, LDA, P1, LDA, Q, LDU, Z, LDU, WORK,
793 $ RWORK, RESULT( 6 ) )
794 CALL ZGET51( 3, N, T, LDA, P1, LDA, Q, LDU, Q, LDU, WORK,
795 $ RWORK, RESULT( 7 ) )
796 CALL ZGET51( 3, N, T, LDA, P1, LDA, Z, LDU, Z, LDU, WORK,
797 $ RWORK, RESULT( 8 ) )
798 *
799 * Compute the Left and Right Eigenvectors of (S1,P1)
800 *
801 * 9: Compute the left eigenvector Matrix without
802 * back transforming:
803 *
804 NTEST = 9
805 RESULT( 9 ) = ULPINV
806 *
807 * To test "SELECT" option, compute half of the eigenvectors
808 * in one call, and half in another
809 *
810 I1 = N / 2
811 DO 120 J = 1, I1
812 LLWORK( J ) = .TRUE.
813 120 CONTINUE
814 DO 130 J = I1 + 1, N
815 LLWORK( J ) = .FALSE.
816 130 CONTINUE
817 *
818 CALL ZTGEVC( 'L', 'S', LLWORK, N, S1, LDA, P1, LDA, EVECTL,
819 $ LDU, CDUMMA, LDU, N, IN, WORK, RWORK, IINFO )
820 IF( IINFO.NE.0 ) THEN
821 WRITE( NOUNIT, FMT = 9999 )'ZTGEVC(L,S1)', IINFO, N,
822 $ JTYPE, IOLDSD
823 INFO = ABS( IINFO )
824 GO TO 210
825 END IF
826 *
827 I1 = IN
828 DO 140 J = 1, I1
829 LLWORK( J ) = .FALSE.
830 140 CONTINUE
831 DO 150 J = I1 + 1, N
832 LLWORK( J ) = .TRUE.
833 150 CONTINUE
834 *
835 CALL ZTGEVC( 'L', 'S', LLWORK, N, S1, LDA, P1, LDA,
836 $ EVECTL( 1, I1+1 ), LDU, CDUMMA, LDU, N, IN,
837 $ WORK, RWORK, IINFO )
838 IF( IINFO.NE.0 ) THEN
839 WRITE( NOUNIT, FMT = 9999 )'ZTGEVC(L,S2)', IINFO, N,
840 $ JTYPE, IOLDSD
841 INFO = ABS( IINFO )
842 GO TO 210
843 END IF
844 *
845 CALL ZGET52( .TRUE., N, S1, LDA, P1, LDA, EVECTL, LDU,
846 $ ALPHA1, BETA1, WORK, RWORK, DUMMA( 1 ) )
847 RESULT( 9 ) = DUMMA( 1 )
848 IF( DUMMA( 2 ).GT.THRSHN ) THEN
849 WRITE( NOUNIT, FMT = 9998 )'Left', 'ZTGEVC(HOWMNY=S)',
850 $ DUMMA( 2 ), N, JTYPE, IOLDSD
851 END IF
852 *
853 * 10: Compute the left eigenvector Matrix with
854 * back transforming:
855 *
856 NTEST = 10
857 RESULT( 10 ) = ULPINV
858 CALL ZLACPY( 'F', N, N, Q, LDU, EVECTL, LDU )
859 CALL ZTGEVC( 'L', 'B', LLWORK, N, S1, LDA, P1, LDA, EVECTL,
860 $ LDU, CDUMMA, LDU, N, IN, WORK, RWORK, IINFO )
861 IF( IINFO.NE.0 ) THEN
862 WRITE( NOUNIT, FMT = 9999 )'ZTGEVC(L,B)', IINFO, N,
863 $ JTYPE, IOLDSD
864 INFO = ABS( IINFO )
865 GO TO 210
866 END IF
867 *
868 CALL ZGET52( .TRUE., N, H, LDA, T, LDA, EVECTL, LDU, ALPHA1,
869 $ BETA1, WORK, RWORK, DUMMA( 1 ) )
870 RESULT( 10 ) = DUMMA( 1 )
871 IF( DUMMA( 2 ).GT.THRSHN ) THEN
872 WRITE( NOUNIT, FMT = 9998 )'Left', 'ZTGEVC(HOWMNY=B)',
873 $ DUMMA( 2 ), N, JTYPE, IOLDSD
874 END IF
875 *
876 * 11: Compute the right eigenvector Matrix without
877 * back transforming:
878 *
879 NTEST = 11
880 RESULT( 11 ) = ULPINV
881 *
882 * To test "SELECT" option, compute half of the eigenvectors
883 * in one call, and half in another
884 *
885 I1 = N / 2
886 DO 160 J = 1, I1
887 LLWORK( J ) = .TRUE.
888 160 CONTINUE
889 DO 170 J = I1 + 1, N
890 LLWORK( J ) = .FALSE.
891 170 CONTINUE
892 *
893 CALL ZTGEVC( 'R', 'S', LLWORK, N, S1, LDA, P1, LDA, CDUMMA,
894 $ LDU, EVECTR, LDU, N, IN, WORK, RWORK, IINFO )
895 IF( IINFO.NE.0 ) THEN
896 WRITE( NOUNIT, FMT = 9999 )'ZTGEVC(R,S1)', IINFO, N,
897 $ JTYPE, IOLDSD
898 INFO = ABS( IINFO )
899 GO TO 210
900 END IF
901 *
902 I1 = IN
903 DO 180 J = 1, I1
904 LLWORK( J ) = .FALSE.
905 180 CONTINUE
906 DO 190 J = I1 + 1, N
907 LLWORK( J ) = .TRUE.
908 190 CONTINUE
909 *
910 CALL ZTGEVC( 'R', 'S', LLWORK, N, S1, LDA, P1, LDA, CDUMMA,
911 $ LDU, EVECTR( 1, I1+1 ), LDU, N, IN, WORK,
912 $ RWORK, IINFO )
913 IF( IINFO.NE.0 ) THEN
914 WRITE( NOUNIT, FMT = 9999 )'ZTGEVC(R,S2)', IINFO, N,
915 $ JTYPE, IOLDSD
916 INFO = ABS( IINFO )
917 GO TO 210
918 END IF
919 *
920 CALL ZGET52( .FALSE., N, S1, LDA, P1, LDA, EVECTR, LDU,
921 $ ALPHA1, BETA1, WORK, RWORK, DUMMA( 1 ) )
922 RESULT( 11 ) = DUMMA( 1 )
923 IF( DUMMA( 2 ).GT.THRESH ) THEN
924 WRITE( NOUNIT, FMT = 9998 )'Right', 'ZTGEVC(HOWMNY=S)',
925 $ DUMMA( 2 ), N, JTYPE, IOLDSD
926 END IF
927 *
928 * 12: Compute the right eigenvector Matrix with
929 * back transforming:
930 *
931 NTEST = 12
932 RESULT( 12 ) = ULPINV
933 CALL ZLACPY( 'F', N, N, Z, LDU, EVECTR, LDU )
934 CALL ZTGEVC( 'R', 'B', LLWORK, N, S1, LDA, P1, LDA, CDUMMA,
935 $ LDU, EVECTR, LDU, N, IN, WORK, RWORK, IINFO )
936 IF( IINFO.NE.0 ) THEN
937 WRITE( NOUNIT, FMT = 9999 )'ZTGEVC(R,B)', IINFO, N,
938 $ JTYPE, IOLDSD
939 INFO = ABS( IINFO )
940 GO TO 210
941 END IF
942 *
943 CALL ZGET52( .FALSE., N, H, LDA, T, LDA, EVECTR, LDU,
944 $ ALPHA1, BETA1, WORK, RWORK, DUMMA( 1 ) )
945 RESULT( 12 ) = DUMMA( 1 )
946 IF( DUMMA( 2 ).GT.THRESH ) THEN
947 WRITE( NOUNIT, FMT = 9998 )'Right', 'ZTGEVC(HOWMNY=B)',
948 $ DUMMA( 2 ), N, JTYPE, IOLDSD
949 END IF
950 *
951 * Tests 13--15 are done only on request
952 *
953 IF( TSTDIF ) THEN
954 *
955 * Do Tests 13--14
956 *
957 CALL ZGET51( 2, N, S1, LDA, S2, LDA, Q, LDU, Z, LDU,
958 $ WORK, RWORK, RESULT( 13 ) )
959 CALL ZGET51( 2, N, P1, LDA, P2, LDA, Q, LDU, Z, LDU,
960 $ WORK, RWORK, RESULT( 14 ) )
961 *
962 * Do Test 15
963 *
964 TEMP1 = ZERO
965 TEMP2 = ZERO
966 DO 200 J = 1, N
967 TEMP1 = MAX( TEMP1, ABS( ALPHA1( J )-ALPHA3( J ) ) )
968 TEMP2 = MAX( TEMP2, ABS( BETA1( J )-BETA3( J ) ) )
969 200 CONTINUE
970 *
971 TEMP1 = TEMP1 / MAX( SAFMIN, ULP*MAX( TEMP1, ANORM ) )
972 TEMP2 = TEMP2 / MAX( SAFMIN, ULP*MAX( TEMP2, BNORM ) )
973 RESULT( 15 ) = MAX( TEMP1, TEMP2 )
974 NTEST = 15
975 ELSE
976 RESULT( 13 ) = ZERO
977 RESULT( 14 ) = ZERO
978 RESULT( 15 ) = ZERO
979 NTEST = 12
980 END IF
981 *
982 * End of Loop -- Check for RESULT(j) > THRESH
983 *
984 210 CONTINUE
985 *
986 NTESTT = NTESTT + NTEST
987 *
988 * Print out tests which fail.
989 *
990 DO 220 JR = 1, NTEST
991 IF( RESULT( JR ).GE.THRESH ) THEN
992 *
993 * If this is the first test to fail,
994 * print a header to the data file.
995 *
996 IF( NERRS.EQ.0 ) THEN
997 WRITE( NOUNIT, FMT = 9997 )'ZGG'
998 *
999 * Matrix types
1000 *
1001 WRITE( NOUNIT, FMT = 9996 )
1002 WRITE( NOUNIT, FMT = 9995 )
1003 WRITE( NOUNIT, FMT = 9994 )'Unitary'
1004 *
1005 * Tests performed
1006 *
1007 WRITE( NOUNIT, FMT = 9993 )'unitary', '*',
1008 $ 'conjugate transpose', ( '*', J = 1, 10 )
1009 *
1010 END IF
1011 NERRS = NERRS + 1
1012 IF( RESULT( JR ).LT.10000.0D0 ) THEN
1013 WRITE( NOUNIT, FMT = 9992 )N, JTYPE, IOLDSD, JR,
1014 $ RESULT( JR )
1015 ELSE
1016 WRITE( NOUNIT, FMT = 9991 )N, JTYPE, IOLDSD, JR,
1017 $ RESULT( JR )
1018 END IF
1019 END IF
1020 220 CONTINUE
1021 *
1022 230 CONTINUE
1023 240 CONTINUE
1024 *
1025 * Summary
1026 *
1027 CALL DLASUM( 'ZGG', NOUNIT, NERRS, NTESTT )
1028 RETURN
1029 *
1030 9999 FORMAT( ' ZCHKGG: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
1031 $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' )
1032 *
1033 9998 FORMAT( ' ZCHKGG: ', A, ' Eigenvectors from ', A, ' incorrectly ',
1034 $ 'normalized.', / ' Bits of error=', 0P, G10.3, ',', 9X,
1035 $ 'N=', I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5,
1036 $ ')' )
1037 *
1038 9997 FORMAT( 1X, A3, ' -- Complex Generalized eigenvalue problem' )
1039 *
1040 9996 FORMAT( ' Matrix types (see ZCHKGG for details): ' )
1041 *
1042 9995 FORMAT( ' Special Matrices:', 23X,
1043 $ '(J''=transposed Jordan block)',
1044 $ / ' 1=(0,0) 2=(I,0) 3=(0,I) 4=(I,I) 5=(J'',J'') ',
1045 $ '6=(diag(J'',I), diag(I,J''))', / ' Diagonal Matrices: ( ',
1046 $ 'D=diag(0,1,2,...) )', / ' 7=(D,I) 9=(large*D, small*I',
1047 $ ') 11=(large*I, small*D) 13=(large*D, large*I)', /
1048 $ ' 8=(I,D) 10=(small*D, large*I) 12=(small*I, large*D) ',
1049 $ ' 14=(small*D, small*I)', / ' 15=(D, reversed D)' )
1050 9994 FORMAT( ' Matrices Rotated by Random ', A, ' Matrices U, V:',
1051 $ / ' 16=Transposed Jordan Blocks 19=geometric ',
1052 $ 'alpha, beta=0,1', / ' 17=arithm. alpha&beta ',
1053 $ ' 20=arithmetic alpha, beta=0,1', / ' 18=clustered ',
1054 $ 'alpha, beta=0,1 21=random alpha, beta=0,1',
1055 $ / ' Large & Small Matrices:', / ' 22=(large, small) ',
1056 $ '23=(small,large) 24=(small,small) 25=(large,large)',
1057 $ / ' 26=random O(1) matrices.' )
1058 *
1059 9993 FORMAT( / ' Tests performed: (H is Hessenberg, S is Schur, B, ',
1060 $ 'T, P are triangular,', / 20X, 'U, V, Q, and Z are ', A,
1061 $ ', l and r are the', / 20X,
1062 $ 'appropriate left and right eigenvectors, resp., a is',
1063 $ / 20X, 'alpha, b is beta, and ', A, ' means ', A, '.)',
1064 $ / ' 1 = | A - U H V', A,
1065 $ ' | / ( |A| n ulp ) 2 = | B - U T V', A,
1066 $ ' | / ( |B| n ulp )', / ' 3 = | I - UU', A,
1067 $ ' | / ( n ulp ) 4 = | I - VV', A,
1068 $ ' | / ( n ulp )', / ' 5 = | H - Q S Z', A,
1069 $ ' | / ( |H| n ulp )', 6X, '6 = | T - Q P Z', A,
1070 $ ' | / ( |T| n ulp )', / ' 7 = | I - QQ', A,
1071 $ ' | / ( n ulp ) 8 = | I - ZZ', A,
1072 $ ' | / ( n ulp )', / ' 9 = max | ( b S - a P )', A,
1073 $ ' l | / const. 10 = max | ( b H - a T )', A,
1074 $ ' l | / const.', /
1075 $ ' 11= max | ( b S - a P ) r | / const. 12 = max | ( b H',
1076 $ ' - a T ) r | / const.', / 1X )
1077 *
1078 9992 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=',
1079 $ 4( I4, ',' ), ' result ', I2, ' is', 0P, F8.2 )
1080 9991 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=',
1081 $ 4( I4, ',' ), ' result ', I2, ' is', 1P, D10.3 )
1082 *
1083 * End of ZCHKGG
1084 *
1085 END
2 $ TSTDIF, THRSHN, NOUNIT, A, LDA, B, H, T, S1,
3 $ S2, P1, P2, U, LDU, V, Q, Z, ALPHA1, BETA1,
4 $ ALPHA3, BETA3, EVECTL, EVECTR, WORK, LWORK,
5 $ RWORK, LLWORK, RESULT, INFO )
6 *
7 * -- LAPACK test routine (version 3.1) --
8 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
9 * November 2006
10 *
11 * .. Scalar Arguments ..
12 LOGICAL TSTDIF
13 INTEGER INFO, LDA, LDU, LWORK, NOUNIT, NSIZES, NTYPES
14 DOUBLE PRECISION THRESH, THRSHN
15 * ..
16 * .. Array Arguments ..
17 LOGICAL DOTYPE( * ), LLWORK( * )
18 INTEGER ISEED( 4 ), NN( * )
19 DOUBLE PRECISION RESULT( 15 ), RWORK( * )
20 COMPLEX*16 A( LDA, * ), ALPHA1( * ), ALPHA3( * ),
21 $ B( LDA, * ), BETA1( * ), BETA3( * ),
22 $ EVECTL( LDU, * ), EVECTR( LDU, * ),
23 $ H( LDA, * ), P1( LDA, * ), P2( LDA, * ),
24 $ Q( LDU, * ), S1( LDA, * ), S2( LDA, * ),
25 $ T( LDA, * ), U( LDU, * ), V( LDU, * ),
26 $ WORK( * ), Z( LDU, * )
27 * ..
28 *
29 * Purpose
30 * =======
31 *
32 * ZCHKGG checks the nonsymmetric generalized eigenvalue problem
33 * routines.
34 * H H H
35 * ZGGHRD factors A and B as U H V and U T V , where means conjugate
36 * transpose, H is hessenberg, T is triangular and U and V are unitary.
37 *
38 * H H
39 * ZHGEQZ factors H and T as Q S Z and Q P Z , where P and S are upper
40 * triangular and Q and Z are unitary. It also computes the generalized
41 * eigenvalues (alpha(1),beta(1)),...,(alpha(n),beta(n)), where
42 * alpha(j)=S(j,j) and beta(j)=P(j,j) -- thus, w(j) = alpha(j)/beta(j)
43 * is a root of the generalized eigenvalue problem
44 *
45 * det( A - w(j) B ) = 0
46 *
47 * and m(j) = beta(j)/alpha(j) is a root of the essentially equivalent
48 * problem
49 *
50 * det( m(j) A - B ) = 0
51 *
52 * ZTGEVC computes the matrix L of left eigenvectors and the matrix R
53 * of right eigenvectors for the matrix pair ( S, P ). In the
54 * description below, l and r are left and right eigenvectors
55 * corresponding to the generalized eigenvalues (alpha,beta).
56 *
57 * When ZCHKGG is called, a number of matrix "sizes" ("n's") and a
58 * number of matrix "types" are specified. For each size ("n")
59 * and each type of matrix, one matrix will be generated and used
60 * to test the nonsymmetric eigenroutines. For each matrix, 13
61 * tests will be performed. The first twelve "test ratios" should be
62 * small -- O(1). They will be compared with the threshhold THRESH:
63 *
64 * H
65 * (1) | A - U H V | / ( |A| n ulp )
66 *
67 * H
68 * (2) | B - U T V | / ( |B| n ulp )
69 *
70 * H
71 * (3) | I - UU | / ( n ulp )
72 *
73 * H
74 * (4) | I - VV | / ( n ulp )
75 *
76 * H
77 * (5) | H - Q S Z | / ( |H| n ulp )
78 *
79 * H
80 * (6) | T - Q P Z | / ( |T| n ulp )
81 *
82 * H
83 * (7) | I - QQ | / ( n ulp )
84 *
85 * H
86 * (8) | I - ZZ | / ( n ulp )
87 *
88 * (9) max over all left eigenvalue/-vector pairs (beta/alpha,l) of
89 * H
90 * | (beta A - alpha B) l | / ( ulp max( |beta A|, |alpha B| ) )
91 *
92 * (10) max over all left eigenvalue/-vector pairs (beta/alpha,l') of
93 * H
94 * | (beta H - alpha T) l' | / ( ulp max( |beta H|, |alpha T| ) )
95 *
96 * where the eigenvectors l' are the result of passing Q to
97 * DTGEVC and back transforming (JOB='B').
98 *
99 * (11) max over all right eigenvalue/-vector pairs (beta/alpha,r) of
100 *
101 * | (beta A - alpha B) r | / ( ulp max( |beta A|, |alpha B| ) )
102 *
103 * (12) max over all right eigenvalue/-vector pairs (beta/alpha,r') of
104 *
105 * | (beta H - alpha T) r' | / ( ulp max( |beta H|, |alpha T| ) )
106 *
107 * where the eigenvectors r' are the result of passing Z to
108 * DTGEVC and back transforming (JOB='B').
109 *
110 * The last three test ratios will usually be small, but there is no
111 * mathematical requirement that they be so. They are therefore
112 * compared with THRESH only if TSTDIF is .TRUE.
113 *
114 * (13) | S(Q,Z computed) - S(Q,Z not computed) | / ( |S| ulp )
115 *
116 * (14) | P(Q,Z computed) - P(Q,Z not computed) | / ( |P| ulp )
117 *
118 * (15) max( |alpha(Q,Z computed) - alpha(Q,Z not computed)|/|S| ,
119 * |beta(Q,Z computed) - beta(Q,Z not computed)|/|P| ) / ulp
120 *
121 * In addition, the normalization of L and R are checked, and compared
122 * with the threshhold THRSHN.
123 *
124 * Test Matrices
125 * ---- --------
126 *
127 * The sizes of the test matrices are specified by an array
128 * NN(1:NSIZES); the value of each element NN(j) specifies one size.
129 * The "types" are specified by a logical array DOTYPE( 1:NTYPES ); if
130 * DOTYPE(j) is .TRUE., then matrix type "j" will be generated.
131 * Currently, the list of possible types is:
132 *
133 * (1) ( 0, 0 ) (a pair of zero matrices)
134 *
135 * (2) ( I, 0 ) (an identity and a zero matrix)
136 *
137 * (3) ( 0, I ) (an identity and a zero matrix)
138 *
139 * (4) ( I, I ) (a pair of identity matrices)
140 *
141 * t t
142 * (5) ( J , J ) (a pair of transposed Jordan blocks)
143 *
144 * t ( I 0 )
145 * (6) ( X, Y ) where X = ( J 0 ) and Y = ( t )
146 * ( 0 I ) ( 0 J )
147 * and I is a k x k identity and J a (k+1)x(k+1)
148 * Jordan block; k=(N-1)/2
149 *
150 * (7) ( D, I ) where D is P*D1, P is a random unitary diagonal
151 * matrix (i.e., with random magnitude 1 entries
152 * on the diagonal), and D1=diag( 0, 1,..., N-1 )
153 * (i.e., a diagonal matrix with D1(1,1)=0,
154 * D1(2,2)=1, ..., D1(N,N)=N-1.)
155 * (8) ( I, D )
156 *
157 * (9) ( big*D, small*I ) where "big" is near overflow and small=1/big
158 *
159 * (10) ( small*D, big*I )
160 *
161 * (11) ( big*I, small*D )
162 *
163 * (12) ( small*I, big*D )
164 *
165 * (13) ( big*D, big*I )
166 *
167 * (14) ( small*D, small*I )
168 *
169 * (15) ( D1, D2 ) where D1=P*diag( 0, 0, 1, ..., N-3, 0 ) and
170 * D2=Q*diag( 0, N-3, N-4,..., 1, 0, 0 ), and
171 * P and Q are random unitary diagonal matrices.
172 * t t
173 * (16) U ( J , J ) V where U and V are random unitary matrices.
174 *
175 * (17) U ( T1, T2 ) V where T1 and T2 are upper triangular matrices
176 * with random O(1) entries above the diagonal
177 * and diagonal entries diag(T1) =
178 * P*( 0, 0, 1, ..., N-3, 0 ) and diag(T2) =
179 * Q*( 0, N-3, N-4,..., 1, 0, 0 )
180 *
181 * (18) U ( T1, T2 ) V diag(T1) = ( 0, 0, 1, 1, s, ..., s, 0 )
182 * diag(T2) = ( 0, 1, 0, 1,..., 1, 0 )
183 * s = machine precision.
184 *
185 * (19) U ( T1, T2 ) V diag(T1)=( 0,0,1,1, 1-d, ..., 1-(N-5)*d=s, 0 )
186 * diag(T2) = ( 0, 1, 0, 1, ..., 1, 0 )
187 *
188 * N-5
189 * (20) U ( T1, T2 ) V diag(T1)=( 0, 0, 1, 1, a, ..., a =s, 0 )
190 * diag(T2) = ( 0, 1, 0, 1, ..., 1, 0, 0 )
191 *
192 * (21) U ( T1, T2 ) V diag(T1)=( 0, 0, 1, r1, r2, ..., r(N-4), 0 )
193 * diag(T2) = ( 0, 1, 0, 1, ..., 1, 0, 0 )
194 * where r1,..., r(N-4) are random.
195 *
196 * (22) U ( big*T1, small*T2 ) V diag(T1) = P*( 0, 0, 1, ..., N-3, 0 )
197 * diag(T2) = ( 0, 1, ..., 1, 0, 0 )
198 *
199 * (23) U ( small*T1, big*T2 ) V diag(T1) = P*( 0, 0, 1, ..., N-3, 0 )
200 * diag(T2) = ( 0, 1, ..., 1, 0, 0 )
201 *
202 * (24) U ( small*T1, small*T2 ) V diag(T1) = P*( 0, 0, 1, ..., N-3, 0 )
203 * diag(T2) = ( 0, 1, ..., 1, 0, 0 )
204 *
205 * (25) U ( big*T1, big*T2 ) V diag(T1) = P*( 0, 0, 1, ..., N-3, 0 )
206 * diag(T2) = ( 0, 1, ..., 1, 0, 0 )
207 *
208 * (26) U ( T1, T2 ) V where T1 and T2 are random upper-triangular
209 * matrices.
210 *
211 * Arguments
212 * =========
213 *
214 * NSIZES (input) INTEGER
215 * The number of sizes of matrices to use. If it is zero,
216 * ZCHKGG does nothing. It must be at least zero.
217 *
218 * NN (input) INTEGER array, dimension (NSIZES)
219 * An array containing the sizes to be used for the matrices.
220 * Zero values will be skipped. The values must be at least
221 * zero.
222 *
223 * NTYPES (input) INTEGER
224 * The number of elements in DOTYPE. If it is zero, ZCHKGG
225 * does nothing. It must be at least zero. If it is MAXTYP+1
226 * and NSIZES is 1, then an additional type, MAXTYP+1 is
227 * defined, which is to use whatever matrix is in A. This
228 * is only useful if DOTYPE(1:MAXTYP) is .FALSE. and
229 * DOTYPE(MAXTYP+1) is .TRUE. .
230 *
231 * DOTYPE (input) LOGICAL array, dimension (NTYPES)
232 * If DOTYPE(j) is .TRUE., then for each size in NN a
233 * matrix of that size and of type j will be generated.
234 * If NTYPES is smaller than the maximum number of types
235 * defined (PARAMETER MAXTYP), then types NTYPES+1 through
236 * MAXTYP will not be generated. If NTYPES is larger
237 * than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES)
238 * will be ignored.
239 *
240 * ISEED (input/output) INTEGER array, dimension (4)
241 * On entry ISEED specifies the seed of the random number
242 * generator. The array elements should be between 0 and 4095;
243 * if not they will be reduced mod 4096. Also, ISEED(4) must
244 * be odd. The random number generator uses a linear
245 * congruential sequence limited to small integers, and so
246 * should produce machine independent random numbers. The
247 * values of ISEED are changed on exit, and can be used in the
248 * next call to ZCHKGG to continue the same random number
249 * sequence.
250 *
251 * THRESH (input) DOUBLE PRECISION
252 * A test will count as "failed" if the "error", computed as
253 * described above, exceeds THRESH. Note that the error
254 * is scaled to be O(1), so THRESH should be a reasonably
255 * small multiple of 1, e.g., 10 or 100. In particular,
256 * it should not depend on the precision (single vs. double)
257 * or the size of the matrix. It must be at least zero.
258 *
259 * TSTDIF (input) LOGICAL
260 * Specifies whether test ratios 13-15 will be computed and
261 * compared with THRESH.
262 * = .FALSE.: Only test ratios 1-12 will be computed and tested.
263 * Ratios 13-15 will be set to zero.
264 * = .TRUE.: All the test ratios 1-15 will be computed and
265 * tested.
266 *
267 * THRSHN (input) DOUBLE PRECISION
268 * Threshhold for reporting eigenvector normalization error.
269 * If the normalization of any eigenvector differs from 1 by
270 * more than THRSHN*ulp, then a special error message will be
271 * printed. (This is handled separately from the other tests,
272 * since only a compiler or programming error should cause an
273 * error message, at least if THRSHN is at least 5--10.)
274 *
275 * NOUNIT (input) INTEGER
276 * The FORTRAN unit number for printing out error messages
277 * (e.g., if a routine returns IINFO not equal to 0.)
278 *
279 * A (input/workspace) COMPLEX*16 array, dimension (LDA, max(NN))
280 * Used to hold the original A matrix. Used as input only
281 * if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and
282 * DOTYPE(MAXTYP+1)=.TRUE.
283 *
284 * LDA (input) INTEGER
285 * The leading dimension of A, B, H, T, S1, P1, S2, and P2.
286 * It must be at least 1 and at least max( NN ).
287 *
288 * B (input/workspace) COMPLEX*16 array, dimension (LDA, max(NN))
289 * Used to hold the original B matrix. Used as input only
290 * if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and
291 * DOTYPE(MAXTYP+1)=.TRUE.
292 *
293 * H (workspace) COMPLEX*16 array, dimension (LDA, max(NN))
294 * The upper Hessenberg matrix computed from A by ZGGHRD.
295 *
296 * T (workspace) COMPLEX*16 array, dimension (LDA, max(NN))
297 * The upper triangular matrix computed from B by ZGGHRD.
298 *
299 * S1 (workspace) COMPLEX*16 array, dimension (LDA, max(NN))
300 * The Schur (upper triangular) matrix computed from H by ZHGEQZ
301 * when Q and Z are also computed.
302 *
303 * S2 (workspace) COMPLEX*16 array, dimension (LDA, max(NN))
304 * The Schur (upper triangular) matrix computed from H by ZHGEQZ
305 * when Q and Z are not computed.
306 *
307 * P1 (workspace) COMPLEX*16 array, dimension (LDA, max(NN))
308 * The upper triangular matrix computed from T by ZHGEQZ
309 * when Q and Z are also computed.
310 *
311 * P2 (workspace) COMPLEX*16 array, dimension (LDA, max(NN))
312 * The upper triangular matrix computed from T by ZHGEQZ
313 * when Q and Z are not computed.
314 *
315 * U (workspace) COMPLEX*16 array, dimension (LDU, max(NN))
316 * The (left) unitary matrix computed by ZGGHRD.
317 *
318 * LDU (input) INTEGER
319 * The leading dimension of U, V, Q, Z, EVECTL, and EVEZTR. It
320 * must be at least 1 and at least max( NN ).
321 *
322 * V (workspace) COMPLEX*16 array, dimension (LDU, max(NN))
323 * The (right) unitary matrix computed by ZGGHRD.
324 *
325 * Q (workspace) COMPLEX*16 array, dimension (LDU, max(NN))
326 * The (left) unitary matrix computed by ZHGEQZ.
327 *
328 * Z (workspace) COMPLEX*16 array, dimension (LDU, max(NN))
329 * The (left) unitary matrix computed by ZHGEQZ.
330 *
331 * ALPHA1 (workspace) COMPLEX*16 array, dimension (max(NN))
332 * BETA1 (workspace) COMPLEX*16 array, dimension (max(NN))
333 * The generalized eigenvalues of (A,B) computed by ZHGEQZ
334 * when Q, Z, and the full Schur matrices are computed.
335 *
336 * ALPHA3 (workspace) COMPLEX*16 array, dimension (max(NN))
337 * BETA3 (workspace) COMPLEX*16 array, dimension (max(NN))
338 * The generalized eigenvalues of (A,B) computed by ZHGEQZ
339 * when neither Q, Z, nor the Schur matrices are computed.
340 *
341 * EVECTL (workspace) COMPLEX*16 array, dimension (LDU, max(NN))
342 * The (lower triangular) left eigenvector matrix for the
343 * matrices in S1 and P1.
344 *
345 * EVEZTR (workspace) COMPLEX*16 array, dimension (LDU, max(NN))
346 * The (upper triangular) right eigenvector matrix for the
347 * matrices in S1 and P1.
348 *
349 * WORK (workspace) COMPLEX*16 array, dimension (LWORK)
350 *
351 * LWORK (input) INTEGER
352 * The number of entries in WORK. This must be at least
353 * max( 4*N, 2 * N**2, 1 ), for all N=NN(j).
354 *
355 * RWORK (workspace) DOUBLE PRECISION array, dimension (2*max(NN))
356 *
357 * LLWORK (workspace) LOGICAL array, dimension (max(NN))
358 *
359 * RESULT (output) DOUBLE PRECISION array, dimension (15)
360 * The values computed by the tests described above.
361 * The values are currently limited to 1/ulp, to avoid
362 * overflow.
363 *
364 * INFO (output) INTEGER
365 * = 0: successful exit.
366 * < 0: if INFO = -i, the i-th argument had an illegal value.
367 * > 0: A routine returned an error code. INFO is the
368 * absolute value of the INFO value returned.
369 *
370 * =====================================================================
371 *
372 * .. Parameters ..
373 DOUBLE PRECISION ZERO, ONE
374 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
375 COMPLEX*16 CZERO, CONE
376 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
377 $ CONE = ( 1.0D+0, 0.0D+0 ) )
378 INTEGER MAXTYP
379 PARAMETER ( MAXTYP = 26 )
380 * ..
381 * .. Local Scalars ..
382 LOGICAL BADNN
383 INTEGER I1, IADD, IINFO, IN, J, JC, JR, JSIZE, JTYPE,
384 $ LWKOPT, MTYPES, N, N1, NERRS, NMATS, NMAX,
385 $ NTEST, NTESTT
386 DOUBLE PRECISION ANORM, BNORM, SAFMAX, SAFMIN, TEMP1, TEMP2,
387 $ ULP, ULPINV
388 COMPLEX*16 CTEMP
389 * ..
390 * .. Local Arrays ..
391 LOGICAL LASIGN( MAXTYP ), LBSIGN( MAXTYP )
392 INTEGER IOLDSD( 4 ), KADD( 6 ), KAMAGN( MAXTYP ),
393 $ KATYPE( MAXTYP ), KAZERO( MAXTYP ),
394 $ KBMAGN( MAXTYP ), KBTYPE( MAXTYP ),
395 $ KBZERO( MAXTYP ), KCLASS( MAXTYP ),
396 $ KTRIAN( MAXTYP ), KZ1( 6 ), KZ2( 6 )
397 DOUBLE PRECISION DUMMA( 4 ), RMAGN( 0: 3 )
398 COMPLEX*16 CDUMMA( 4 )
399 * ..
400 * .. External Functions ..
401 DOUBLE PRECISION DLAMCH, ZLANGE
402 COMPLEX*16 ZLARND
403 EXTERNAL DLAMCH, ZLANGE, ZLARND
404 * ..
405 * .. External Subroutines ..
406 EXTERNAL DLABAD, DLASUM, XERBLA, ZGEQR2, ZGET51, ZGET52,
407 $ ZGGHRD, ZHGEQZ, ZLACPY, ZLARFG, ZLASET, ZLATM4,
408 $ ZTGEVC, ZUNM2R
409 * ..
410 * .. Intrinsic Functions ..
411 INTRINSIC ABS, DBLE, DCONJG, MAX, MIN, SIGN
412 * ..
413 * .. Data statements ..
414 DATA KCLASS / 15*1, 10*2, 1*3 /
415 DATA KZ1 / 0, 1, 2, 1, 3, 3 /
416 DATA KZ2 / 0, 0, 1, 2, 1, 1 /
417 DATA KADD / 0, 0, 0, 0, 3, 2 /
418 DATA KATYPE / 0, 1, 0, 1, 2, 3, 4, 1, 4, 4, 1, 1, 4,
419 $ 4, 4, 2, 4, 5, 8, 7, 9, 4*4, 0 /
420 DATA KBTYPE / 0, 0, 1, 1, 2, -3, 1, 4, 1, 1, 4, 4,
421 $ 1, 1, -4, 2, -4, 8*8, 0 /
422 DATA KAZERO / 6*1, 2, 1, 2*2, 2*1, 2*2, 3, 1, 3,
423 $ 4*5, 4*3, 1 /
424 DATA KBZERO / 6*1, 1, 2, 2*1, 2*2, 2*1, 4, 1, 4,
425 $ 4*6, 4*4, 1 /
426 DATA KAMAGN / 8*1, 2, 3, 2, 3, 2, 3, 7*1, 2, 3, 3,
427 $ 2, 1 /
428 DATA KBMAGN / 8*1, 3, 2, 3, 2, 2, 3, 7*1, 3, 2, 3,
429 $ 2, 1 /
430 DATA KTRIAN / 16*0, 10*1 /
431 DATA LASIGN / 6*.FALSE., .TRUE., .FALSE., 2*.TRUE.,
432 $ 2*.FALSE., 3*.TRUE., .FALSE., .TRUE.,
433 $ 3*.FALSE., 5*.TRUE., .FALSE. /
434 DATA LBSIGN / 7*.FALSE., .TRUE., 2*.FALSE.,
435 $ 2*.TRUE., 2*.FALSE., .TRUE., .FALSE., .TRUE.,
436 $ 9*.FALSE. /
437 * ..
438 * .. Executable Statements ..
439 *
440 * Check for errors
441 *
442 INFO = 0
443 *
444 BADNN = .FALSE.
445 NMAX = 1
446 DO 10 J = 1, NSIZES
447 NMAX = MAX( NMAX, NN( J ) )
448 IF( NN( J ).LT.0 )
449 $ BADNN = .TRUE.
450 10 CONTINUE
451 *
452 LWKOPT = MAX( 2*NMAX*NMAX, 4*NMAX, 1 )
453 *
454 * Check for errors
455 *
456 IF( NSIZES.LT.0 ) THEN
457 INFO = -1
458 ELSE IF( BADNN ) THEN
459 INFO = -2
460 ELSE IF( NTYPES.LT.0 ) THEN
461 INFO = -3
462 ELSE IF( THRESH.LT.ZERO ) THEN
463 INFO = -6
464 ELSE IF( LDA.LE.1 .OR. LDA.LT.NMAX ) THEN
465 INFO = -10
466 ELSE IF( LDU.LE.1 .OR. LDU.LT.NMAX ) THEN
467 INFO = -19
468 ELSE IF( LWKOPT.GT.LWORK ) THEN
469 INFO = -30
470 END IF
471 *
472 IF( INFO.NE.0 ) THEN
473 CALL XERBLA( 'ZCHKGG', -INFO )
474 RETURN
475 END IF
476 *
477 * Quick return if possible
478 *
479 IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 )
480 $ RETURN
481 *
482 SAFMIN = DLAMCH( 'Safe minimum' )
483 ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' )
484 SAFMIN = SAFMIN / ULP
485 SAFMAX = ONE / SAFMIN
486 CALL DLABAD( SAFMIN, SAFMAX )
487 ULPINV = ONE / ULP
488 *
489 * The values RMAGN(2:3) depend on N, see below.
490 *
491 RMAGN( 0 ) = ZERO
492 RMAGN( 1 ) = ONE
493 *
494 * Loop over sizes, types
495 *
496 NTESTT = 0
497 NERRS = 0
498 NMATS = 0
499 *
500 DO 240 JSIZE = 1, NSIZES
501 N = NN( JSIZE )
502 N1 = MAX( 1, N )
503 RMAGN( 2 ) = SAFMAX*ULP / DBLE( N1 )
504 RMAGN( 3 ) = SAFMIN*ULPINV*N1
505 *
506 IF( NSIZES.NE.1 ) THEN
507 MTYPES = MIN( MAXTYP, NTYPES )
508 ELSE
509 MTYPES = MIN( MAXTYP+1, NTYPES )
510 END IF
511 *
512 DO 230 JTYPE = 1, MTYPES
513 IF( .NOT.DOTYPE( JTYPE ) )
514 $ GO TO 230
515 NMATS = NMATS + 1
516 NTEST = 0
517 *
518 * Save ISEED in case of an error.
519 *
520 DO 20 J = 1, 4
521 IOLDSD( J ) = ISEED( J )
522 20 CONTINUE
523 *
524 * Initialize RESULT
525 *
526 DO 30 J = 1, 15
527 RESULT( J ) = ZERO
528 30 CONTINUE
529 *
530 * Compute A and B
531 *
532 * Description of control parameters:
533 *
534 * KZLASS: =1 means w/o rotation, =2 means w/ rotation,
535 * =3 means random.
536 * KATYPE: the "type" to be passed to ZLATM4 for computing A.
537 * KAZERO: the pattern of zeros on the diagonal for A:
538 * =1: ( xxx ), =2: (0, xxx ) =3: ( 0, 0, xxx, 0 ),
539 * =4: ( 0, xxx, 0, 0 ), =5: ( 0, 0, 1, xxx, 0 ),
540 * =6: ( 0, 1, 0, xxx, 0 ). (xxx means a string of
541 * non-zero entries.)
542 * KAMAGN: the magnitude of the matrix: =0: zero, =1: O(1),
543 * =2: large, =3: small.
544 * LASIGN: .TRUE. if the diagonal elements of A are to be
545 * multiplied by a random magnitude 1 number.
546 * KBTYPE, KBZERO, KBMAGN, LBSIGN: the same, but for B.
547 * KTRIAN: =0: don't fill in the upper triangle, =1: do.
548 * KZ1, KZ2, KADD: used to implement KAZERO and KBZERO.
549 * RMAGN: used to implement KAMAGN and KBMAGN.
550 *
551 IF( MTYPES.GT.MAXTYP )
552 $ GO TO 110
553 IINFO = 0
554 IF( KCLASS( JTYPE ).LT.3 ) THEN
555 *
556 * Generate A (w/o rotation)
557 *
558 IF( ABS( KATYPE( JTYPE ) ).EQ.3 ) THEN
559 IN = 2*( ( N-1 ) / 2 ) + 1
560 IF( IN.NE.N )
561 $ CALL ZLASET( 'Full', N, N, CZERO, CZERO, A, LDA )
562 ELSE
563 IN = N
564 END IF
565 CALL ZLATM4( KATYPE( JTYPE ), IN, KZ1( KAZERO( JTYPE ) ),
566 $ KZ2( KAZERO( JTYPE ) ), LASIGN( JTYPE ),
567 $ RMAGN( KAMAGN( JTYPE ) ), ULP,
568 $ RMAGN( KTRIAN( JTYPE )*KAMAGN( JTYPE ) ), 4,
569 $ ISEED, A, LDA )
570 IADD = KADD( KAZERO( JTYPE ) )
571 IF( IADD.GT.0 .AND. IADD.LE.N )
572 $ A( IADD, IADD ) = RMAGN( KAMAGN( JTYPE ) )
573 *
574 * Generate B (w/o rotation)
575 *
576 IF( ABS( KBTYPE( JTYPE ) ).EQ.3 ) THEN
577 IN = 2*( ( N-1 ) / 2 ) + 1
578 IF( IN.NE.N )
579 $ CALL ZLASET( 'Full', N, N, CZERO, CZERO, B, LDA )
580 ELSE
581 IN = N
582 END IF
583 CALL ZLATM4( KBTYPE( JTYPE ), IN, KZ1( KBZERO( JTYPE ) ),
584 $ KZ2( KBZERO( JTYPE ) ), LBSIGN( JTYPE ),
585 $ RMAGN( KBMAGN( JTYPE ) ), ONE,
586 $ RMAGN( KTRIAN( JTYPE )*KBMAGN( JTYPE ) ), 4,
587 $ ISEED, B, LDA )
588 IADD = KADD( KBZERO( JTYPE ) )
589 IF( IADD.NE.0 )
590 $ B( IADD, IADD ) = RMAGN( KBMAGN( JTYPE ) )
591 *
592 IF( KCLASS( JTYPE ).EQ.2 .AND. N.GT.0 ) THEN
593 *
594 * Include rotations
595 *
596 * Generate U, V as Householder transformations times a
597 * diagonal matrix. (Note that ZLARFG makes U(j,j) and
598 * V(j,j) real.)
599 *
600 DO 50 JC = 1, N - 1
601 DO 40 JR = JC, N
602 U( JR, JC ) = ZLARND( 3, ISEED )
603 V( JR, JC ) = ZLARND( 3, ISEED )
604 40 CONTINUE
605 CALL ZLARFG( N+1-JC, U( JC, JC ), U( JC+1, JC ), 1,
606 $ WORK( JC ) )
607 WORK( 2*N+JC ) = SIGN( ONE, DBLE( U( JC, JC ) ) )
608 U( JC, JC ) = CONE
609 CALL ZLARFG( N+1-JC, V( JC, JC ), V( JC+1, JC ), 1,
610 $ WORK( N+JC ) )
611 WORK( 3*N+JC ) = SIGN( ONE, DBLE( V( JC, JC ) ) )
612 V( JC, JC ) = CONE
613 50 CONTINUE
614 CTEMP = ZLARND( 3, ISEED )
615 U( N, N ) = CONE
616 WORK( N ) = CZERO
617 WORK( 3*N ) = CTEMP / ABS( CTEMP )
618 CTEMP = ZLARND( 3, ISEED )
619 V( N, N ) = CONE
620 WORK( 2*N ) = CZERO
621 WORK( 4*N ) = CTEMP / ABS( CTEMP )
622 *
623 * Apply the diagonal matrices
624 *
625 DO 70 JC = 1, N
626 DO 60 JR = 1, N
627 A( JR, JC ) = WORK( 2*N+JR )*
628 $ DCONJG( WORK( 3*N+JC ) )*
629 $ A( JR, JC )
630 B( JR, JC ) = WORK( 2*N+JR )*
631 $ DCONJG( WORK( 3*N+JC ) )*
632 $ B( JR, JC )
633 60 CONTINUE
634 70 CONTINUE
635 CALL ZUNM2R( 'L', 'N', N, N, N-1, U, LDU, WORK, A,
636 $ LDA, WORK( 2*N+1 ), IINFO )
637 IF( IINFO.NE.0 )
638 $ GO TO 100
639 CALL ZUNM2R( 'R', 'C', N, N, N-1, V, LDU, WORK( N+1 ),
640 $ A, LDA, WORK( 2*N+1 ), IINFO )
641 IF( IINFO.NE.0 )
642 $ GO TO 100
643 CALL ZUNM2R( 'L', 'N', N, N, N-1, U, LDU, WORK, B,
644 $ LDA, WORK( 2*N+1 ), IINFO )
645 IF( IINFO.NE.0 )
646 $ GO TO 100
647 CALL ZUNM2R( 'R', 'C', N, N, N-1, V, LDU, WORK( N+1 ),
648 $ B, LDA, WORK( 2*N+1 ), IINFO )
649 IF( IINFO.NE.0 )
650 $ GO TO 100
651 END IF
652 ELSE
653 *
654 * Random matrices
655 *
656 DO 90 JC = 1, N
657 DO 80 JR = 1, N
658 A( JR, JC ) = RMAGN( KAMAGN( JTYPE ) )*
659 $ ZLARND( 4, ISEED )
660 B( JR, JC ) = RMAGN( KBMAGN( JTYPE ) )*
661 $ ZLARND( 4, ISEED )
662 80 CONTINUE
663 90 CONTINUE
664 END IF
665 *
666 ANORM = ZLANGE( '1', N, N, A, LDA, RWORK )
667 BNORM = ZLANGE( '1', N, N, B, LDA, RWORK )
668 *
669 100 CONTINUE
670 *
671 IF( IINFO.NE.0 ) THEN
672 WRITE( NOUNIT, FMT = 9999 )'Generator', IINFO, N, JTYPE,
673 $ IOLDSD
674 INFO = ABS( IINFO )
675 RETURN
676 END IF
677 *
678 110 CONTINUE
679 *
680 * Call ZGEQR2, ZUNM2R, and ZGGHRD to compute H, T, U, and V
681 *
682 CALL ZLACPY( ' ', N, N, A, LDA, H, LDA )
683 CALL ZLACPY( ' ', N, N, B, LDA, T, LDA )
684 NTEST = 1
685 RESULT( 1 ) = ULPINV
686 *
687 CALL ZGEQR2( N, N, T, LDA, WORK, WORK( N+1 ), IINFO )
688 IF( IINFO.NE.0 ) THEN
689 WRITE( NOUNIT, FMT = 9999 )'ZGEQR2', IINFO, N, JTYPE,
690 $ IOLDSD
691 INFO = ABS( IINFO )
692 GO TO 210
693 END IF
694 *
695 CALL ZUNM2R( 'L', 'C', N, N, N, T, LDA, WORK, H, LDA,
696 $ WORK( N+1 ), IINFO )
697 IF( IINFO.NE.0 ) THEN
698 WRITE( NOUNIT, FMT = 9999 )'ZUNM2R', IINFO, N, JTYPE,
699 $ IOLDSD
700 INFO = ABS( IINFO )
701 GO TO 210
702 END IF
703 *
704 CALL ZLASET( 'Full', N, N, CZERO, CONE, U, LDU )
705 CALL ZUNM2R( 'R', 'N', N, N, N, T, LDA, WORK, U, LDU,
706 $ WORK( N+1 ), IINFO )
707 IF( IINFO.NE.0 ) THEN
708 WRITE( NOUNIT, FMT = 9999 )'ZUNM2R', IINFO, N, JTYPE,
709 $ IOLDSD
710 INFO = ABS( IINFO )
711 GO TO 210
712 END IF
713 *
714 CALL ZGGHRD( 'V', 'I', N, 1, N, H, LDA, T, LDA, U, LDU, V,
715 $ LDU, IINFO )
716 IF( IINFO.NE.0 ) THEN
717 WRITE( NOUNIT, FMT = 9999 )'ZGGHRD', IINFO, N, JTYPE,
718 $ IOLDSD
719 INFO = ABS( IINFO )
720 GO TO 210
721 END IF
722 NTEST = 4
723 *
724 * Do tests 1--4
725 *
726 CALL ZGET51( 1, N, A, LDA, H, LDA, U, LDU, V, LDU, WORK,
727 $ RWORK, RESULT( 1 ) )
728 CALL ZGET51( 1, N, B, LDA, T, LDA, U, LDU, V, LDU, WORK,
729 $ RWORK, RESULT( 2 ) )
730 CALL ZGET51( 3, N, B, LDA, T, LDA, U, LDU, U, LDU, WORK,
731 $ RWORK, RESULT( 3 ) )
732 CALL ZGET51( 3, N, B, LDA, T, LDA, V, LDU, V, LDU, WORK,
733 $ RWORK, RESULT( 4 ) )
734 *
735 * Call ZHGEQZ to compute S1, P1, S2, P2, Q, and Z, do tests.
736 *
737 * Compute T1 and UZ
738 *
739 * Eigenvalues only
740 *
741 CALL ZLACPY( ' ', N, N, H, LDA, S2, LDA )
742 CALL ZLACPY( ' ', N, N, T, LDA, P2, LDA )
743 NTEST = 5
744 RESULT( 5 ) = ULPINV
745 *
746 CALL ZHGEQZ( 'E', 'N', 'N', N, 1, N, S2, LDA, P2, LDA,
747 $ ALPHA3, BETA3, Q, LDU, Z, LDU, WORK, LWORK,
748 $ RWORK, IINFO )
749 IF( IINFO.NE.0 ) THEN
750 WRITE( NOUNIT, FMT = 9999 )'ZHGEQZ(E)', IINFO, N, JTYPE,
751 $ IOLDSD
752 INFO = ABS( IINFO )
753 GO TO 210
754 END IF
755 *
756 * Eigenvalues and Full Schur Form
757 *
758 CALL ZLACPY( ' ', N, N, H, LDA, S2, LDA )
759 CALL ZLACPY( ' ', N, N, T, LDA, P2, LDA )
760 *
761 CALL ZHGEQZ( 'S', 'N', 'N', N, 1, N, S2, LDA, P2, LDA,
762 $ ALPHA1, BETA1, Q, LDU, Z, LDU, WORK, LWORK,
763 $ RWORK, IINFO )
764 IF( IINFO.NE.0 ) THEN
765 WRITE( NOUNIT, FMT = 9999 )'ZHGEQZ(S)', IINFO, N, JTYPE,
766 $ IOLDSD
767 INFO = ABS( IINFO )
768 GO TO 210
769 END IF
770 *
771 * Eigenvalues, Schur Form, and Schur Vectors
772 *
773 CALL ZLACPY( ' ', N, N, H, LDA, S1, LDA )
774 CALL ZLACPY( ' ', N, N, T, LDA, P1, LDA )
775 *
776 CALL ZHGEQZ( 'S', 'I', 'I', N, 1, N, S1, LDA, P1, LDA,
777 $ ALPHA1, BETA1, Q, LDU, Z, LDU, WORK, LWORK,
778 $ RWORK, IINFO )
779 IF( IINFO.NE.0 ) THEN
780 WRITE( NOUNIT, FMT = 9999 )'ZHGEQZ(V)', IINFO, N, JTYPE,
781 $ IOLDSD
782 INFO = ABS( IINFO )
783 GO TO 210
784 END IF
785 *
786 NTEST = 8
787 *
788 * Do Tests 5--8
789 *
790 CALL ZGET51( 1, N, H, LDA, S1, LDA, Q, LDU, Z, LDU, WORK,
791 $ RWORK, RESULT( 5 ) )
792 CALL ZGET51( 1, N, T, LDA, P1, LDA, Q, LDU, Z, LDU, WORK,
793 $ RWORK, RESULT( 6 ) )
794 CALL ZGET51( 3, N, T, LDA, P1, LDA, Q, LDU, Q, LDU, WORK,
795 $ RWORK, RESULT( 7 ) )
796 CALL ZGET51( 3, N, T, LDA, P1, LDA, Z, LDU, Z, LDU, WORK,
797 $ RWORK, RESULT( 8 ) )
798 *
799 * Compute the Left and Right Eigenvectors of (S1,P1)
800 *
801 * 9: Compute the left eigenvector Matrix without
802 * back transforming:
803 *
804 NTEST = 9
805 RESULT( 9 ) = ULPINV
806 *
807 * To test "SELECT" option, compute half of the eigenvectors
808 * in one call, and half in another
809 *
810 I1 = N / 2
811 DO 120 J = 1, I1
812 LLWORK( J ) = .TRUE.
813 120 CONTINUE
814 DO 130 J = I1 + 1, N
815 LLWORK( J ) = .FALSE.
816 130 CONTINUE
817 *
818 CALL ZTGEVC( 'L', 'S', LLWORK, N, S1, LDA, P1, LDA, EVECTL,
819 $ LDU, CDUMMA, LDU, N, IN, WORK, RWORK, IINFO )
820 IF( IINFO.NE.0 ) THEN
821 WRITE( NOUNIT, FMT = 9999 )'ZTGEVC(L,S1)', IINFO, N,
822 $ JTYPE, IOLDSD
823 INFO = ABS( IINFO )
824 GO TO 210
825 END IF
826 *
827 I1 = IN
828 DO 140 J = 1, I1
829 LLWORK( J ) = .FALSE.
830 140 CONTINUE
831 DO 150 J = I1 + 1, N
832 LLWORK( J ) = .TRUE.
833 150 CONTINUE
834 *
835 CALL ZTGEVC( 'L', 'S', LLWORK, N, S1, LDA, P1, LDA,
836 $ EVECTL( 1, I1+1 ), LDU, CDUMMA, LDU, N, IN,
837 $ WORK, RWORK, IINFO )
838 IF( IINFO.NE.0 ) THEN
839 WRITE( NOUNIT, FMT = 9999 )'ZTGEVC(L,S2)', IINFO, N,
840 $ JTYPE, IOLDSD
841 INFO = ABS( IINFO )
842 GO TO 210
843 END IF
844 *
845 CALL ZGET52( .TRUE., N, S1, LDA, P1, LDA, EVECTL, LDU,
846 $ ALPHA1, BETA1, WORK, RWORK, DUMMA( 1 ) )
847 RESULT( 9 ) = DUMMA( 1 )
848 IF( DUMMA( 2 ).GT.THRSHN ) THEN
849 WRITE( NOUNIT, FMT = 9998 )'Left', 'ZTGEVC(HOWMNY=S)',
850 $ DUMMA( 2 ), N, JTYPE, IOLDSD
851 END IF
852 *
853 * 10: Compute the left eigenvector Matrix with
854 * back transforming:
855 *
856 NTEST = 10
857 RESULT( 10 ) = ULPINV
858 CALL ZLACPY( 'F', N, N, Q, LDU, EVECTL, LDU )
859 CALL ZTGEVC( 'L', 'B', LLWORK, N, S1, LDA, P1, LDA, EVECTL,
860 $ LDU, CDUMMA, LDU, N, IN, WORK, RWORK, IINFO )
861 IF( IINFO.NE.0 ) THEN
862 WRITE( NOUNIT, FMT = 9999 )'ZTGEVC(L,B)', IINFO, N,
863 $ JTYPE, IOLDSD
864 INFO = ABS( IINFO )
865 GO TO 210
866 END IF
867 *
868 CALL ZGET52( .TRUE., N, H, LDA, T, LDA, EVECTL, LDU, ALPHA1,
869 $ BETA1, WORK, RWORK, DUMMA( 1 ) )
870 RESULT( 10 ) = DUMMA( 1 )
871 IF( DUMMA( 2 ).GT.THRSHN ) THEN
872 WRITE( NOUNIT, FMT = 9998 )'Left', 'ZTGEVC(HOWMNY=B)',
873 $ DUMMA( 2 ), N, JTYPE, IOLDSD
874 END IF
875 *
876 * 11: Compute the right eigenvector Matrix without
877 * back transforming:
878 *
879 NTEST = 11
880 RESULT( 11 ) = ULPINV
881 *
882 * To test "SELECT" option, compute half of the eigenvectors
883 * in one call, and half in another
884 *
885 I1 = N / 2
886 DO 160 J = 1, I1
887 LLWORK( J ) = .TRUE.
888 160 CONTINUE
889 DO 170 J = I1 + 1, N
890 LLWORK( J ) = .FALSE.
891 170 CONTINUE
892 *
893 CALL ZTGEVC( 'R', 'S', LLWORK, N, S1, LDA, P1, LDA, CDUMMA,
894 $ LDU, EVECTR, LDU, N, IN, WORK, RWORK, IINFO )
895 IF( IINFO.NE.0 ) THEN
896 WRITE( NOUNIT, FMT = 9999 )'ZTGEVC(R,S1)', IINFO, N,
897 $ JTYPE, IOLDSD
898 INFO = ABS( IINFO )
899 GO TO 210
900 END IF
901 *
902 I1 = IN
903 DO 180 J = 1, I1
904 LLWORK( J ) = .FALSE.
905 180 CONTINUE
906 DO 190 J = I1 + 1, N
907 LLWORK( J ) = .TRUE.
908 190 CONTINUE
909 *
910 CALL ZTGEVC( 'R', 'S', LLWORK, N, S1, LDA, P1, LDA, CDUMMA,
911 $ LDU, EVECTR( 1, I1+1 ), LDU, N, IN, WORK,
912 $ RWORK, IINFO )
913 IF( IINFO.NE.0 ) THEN
914 WRITE( NOUNIT, FMT = 9999 )'ZTGEVC(R,S2)', IINFO, N,
915 $ JTYPE, IOLDSD
916 INFO = ABS( IINFO )
917 GO TO 210
918 END IF
919 *
920 CALL ZGET52( .FALSE., N, S1, LDA, P1, LDA, EVECTR, LDU,
921 $ ALPHA1, BETA1, WORK, RWORK, DUMMA( 1 ) )
922 RESULT( 11 ) = DUMMA( 1 )
923 IF( DUMMA( 2 ).GT.THRESH ) THEN
924 WRITE( NOUNIT, FMT = 9998 )'Right', 'ZTGEVC(HOWMNY=S)',
925 $ DUMMA( 2 ), N, JTYPE, IOLDSD
926 END IF
927 *
928 * 12: Compute the right eigenvector Matrix with
929 * back transforming:
930 *
931 NTEST = 12
932 RESULT( 12 ) = ULPINV
933 CALL ZLACPY( 'F', N, N, Z, LDU, EVECTR, LDU )
934 CALL ZTGEVC( 'R', 'B', LLWORK, N, S1, LDA, P1, LDA, CDUMMA,
935 $ LDU, EVECTR, LDU, N, IN, WORK, RWORK, IINFO )
936 IF( IINFO.NE.0 ) THEN
937 WRITE( NOUNIT, FMT = 9999 )'ZTGEVC(R,B)', IINFO, N,
938 $ JTYPE, IOLDSD
939 INFO = ABS( IINFO )
940 GO TO 210
941 END IF
942 *
943 CALL ZGET52( .FALSE., N, H, LDA, T, LDA, EVECTR, LDU,
944 $ ALPHA1, BETA1, WORK, RWORK, DUMMA( 1 ) )
945 RESULT( 12 ) = DUMMA( 1 )
946 IF( DUMMA( 2 ).GT.THRESH ) THEN
947 WRITE( NOUNIT, FMT = 9998 )'Right', 'ZTGEVC(HOWMNY=B)',
948 $ DUMMA( 2 ), N, JTYPE, IOLDSD
949 END IF
950 *
951 * Tests 13--15 are done only on request
952 *
953 IF( TSTDIF ) THEN
954 *
955 * Do Tests 13--14
956 *
957 CALL ZGET51( 2, N, S1, LDA, S2, LDA, Q, LDU, Z, LDU,
958 $ WORK, RWORK, RESULT( 13 ) )
959 CALL ZGET51( 2, N, P1, LDA, P2, LDA, Q, LDU, Z, LDU,
960 $ WORK, RWORK, RESULT( 14 ) )
961 *
962 * Do Test 15
963 *
964 TEMP1 = ZERO
965 TEMP2 = ZERO
966 DO 200 J = 1, N
967 TEMP1 = MAX( TEMP1, ABS( ALPHA1( J )-ALPHA3( J ) ) )
968 TEMP2 = MAX( TEMP2, ABS( BETA1( J )-BETA3( J ) ) )
969 200 CONTINUE
970 *
971 TEMP1 = TEMP1 / MAX( SAFMIN, ULP*MAX( TEMP1, ANORM ) )
972 TEMP2 = TEMP2 / MAX( SAFMIN, ULP*MAX( TEMP2, BNORM ) )
973 RESULT( 15 ) = MAX( TEMP1, TEMP2 )
974 NTEST = 15
975 ELSE
976 RESULT( 13 ) = ZERO
977 RESULT( 14 ) = ZERO
978 RESULT( 15 ) = ZERO
979 NTEST = 12
980 END IF
981 *
982 * End of Loop -- Check for RESULT(j) > THRESH
983 *
984 210 CONTINUE
985 *
986 NTESTT = NTESTT + NTEST
987 *
988 * Print out tests which fail.
989 *
990 DO 220 JR = 1, NTEST
991 IF( RESULT( JR ).GE.THRESH ) THEN
992 *
993 * If this is the first test to fail,
994 * print a header to the data file.
995 *
996 IF( NERRS.EQ.0 ) THEN
997 WRITE( NOUNIT, FMT = 9997 )'ZGG'
998 *
999 * Matrix types
1000 *
1001 WRITE( NOUNIT, FMT = 9996 )
1002 WRITE( NOUNIT, FMT = 9995 )
1003 WRITE( NOUNIT, FMT = 9994 )'Unitary'
1004 *
1005 * Tests performed
1006 *
1007 WRITE( NOUNIT, FMT = 9993 )'unitary', '*',
1008 $ 'conjugate transpose', ( '*', J = 1, 10 )
1009 *
1010 END IF
1011 NERRS = NERRS + 1
1012 IF( RESULT( JR ).LT.10000.0D0 ) THEN
1013 WRITE( NOUNIT, FMT = 9992 )N, JTYPE, IOLDSD, JR,
1014 $ RESULT( JR )
1015 ELSE
1016 WRITE( NOUNIT, FMT = 9991 )N, JTYPE, IOLDSD, JR,
1017 $ RESULT( JR )
1018 END IF
1019 END IF
1020 220 CONTINUE
1021 *
1022 230 CONTINUE
1023 240 CONTINUE
1024 *
1025 * Summary
1026 *
1027 CALL DLASUM( 'ZGG', NOUNIT, NERRS, NTESTT )
1028 RETURN
1029 *
1030 9999 FORMAT( ' ZCHKGG: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
1031 $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' )
1032 *
1033 9998 FORMAT( ' ZCHKGG: ', A, ' Eigenvectors from ', A, ' incorrectly ',
1034 $ 'normalized.', / ' Bits of error=', 0P, G10.3, ',', 9X,
1035 $ 'N=', I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5,
1036 $ ')' )
1037 *
1038 9997 FORMAT( 1X, A3, ' -- Complex Generalized eigenvalue problem' )
1039 *
1040 9996 FORMAT( ' Matrix types (see ZCHKGG for details): ' )
1041 *
1042 9995 FORMAT( ' Special Matrices:', 23X,
1043 $ '(J''=transposed Jordan block)',
1044 $ / ' 1=(0,0) 2=(I,0) 3=(0,I) 4=(I,I) 5=(J'',J'') ',
1045 $ '6=(diag(J'',I), diag(I,J''))', / ' Diagonal Matrices: ( ',
1046 $ 'D=diag(0,1,2,...) )', / ' 7=(D,I) 9=(large*D, small*I',
1047 $ ') 11=(large*I, small*D) 13=(large*D, large*I)', /
1048 $ ' 8=(I,D) 10=(small*D, large*I) 12=(small*I, large*D) ',
1049 $ ' 14=(small*D, small*I)', / ' 15=(D, reversed D)' )
1050 9994 FORMAT( ' Matrices Rotated by Random ', A, ' Matrices U, V:',
1051 $ / ' 16=Transposed Jordan Blocks 19=geometric ',
1052 $ 'alpha, beta=0,1', / ' 17=arithm. alpha&beta ',
1053 $ ' 20=arithmetic alpha, beta=0,1', / ' 18=clustered ',
1054 $ 'alpha, beta=0,1 21=random alpha, beta=0,1',
1055 $ / ' Large & Small Matrices:', / ' 22=(large, small) ',
1056 $ '23=(small,large) 24=(small,small) 25=(large,large)',
1057 $ / ' 26=random O(1) matrices.' )
1058 *
1059 9993 FORMAT( / ' Tests performed: (H is Hessenberg, S is Schur, B, ',
1060 $ 'T, P are triangular,', / 20X, 'U, V, Q, and Z are ', A,
1061 $ ', l and r are the', / 20X,
1062 $ 'appropriate left and right eigenvectors, resp., a is',
1063 $ / 20X, 'alpha, b is beta, and ', A, ' means ', A, '.)',
1064 $ / ' 1 = | A - U H V', A,
1065 $ ' | / ( |A| n ulp ) 2 = | B - U T V', A,
1066 $ ' | / ( |B| n ulp )', / ' 3 = | I - UU', A,
1067 $ ' | / ( n ulp ) 4 = | I - VV', A,
1068 $ ' | / ( n ulp )', / ' 5 = | H - Q S Z', A,
1069 $ ' | / ( |H| n ulp )', 6X, '6 = | T - Q P Z', A,
1070 $ ' | / ( |T| n ulp )', / ' 7 = | I - QQ', A,
1071 $ ' | / ( n ulp ) 8 = | I - ZZ', A,
1072 $ ' | / ( n ulp )', / ' 9 = max | ( b S - a P )', A,
1073 $ ' l | / const. 10 = max | ( b H - a T )', A,
1074 $ ' l | / const.', /
1075 $ ' 11= max | ( b S - a P ) r | / const. 12 = max | ( b H',
1076 $ ' - a T ) r | / const.', / 1X )
1077 *
1078 9992 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=',
1079 $ 4( I4, ',' ), ' result ', I2, ' is', 0P, F8.2 )
1080 9991 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=',
1081 $ 4( I4, ',' ), ' result ', I2, ' is', 1P, D10.3 )
1082 *
1083 * End of ZCHKGG
1084 *
1085 END