1 SUBROUTINE ZDRGEV( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
2 $ NOUNIT, A, LDA, B, S, T, Q, LDQ, Z, QE, LDQE,
3 $ ALPHA, BETA, ALPHA1, BETA1, WORK, LWORK, RWORK,
4 $ RESULT, INFO )
5 *
6 * -- LAPACK test routine (version 3.1.1) --
7 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
8 * February 2007
9 *
10 * .. Scalar Arguments ..
11 INTEGER INFO, LDA, LDQ, LDQE, LWORK, NOUNIT, NSIZES,
12 $ NTYPES
13 DOUBLE PRECISION THRESH
14 * ..
15 * .. Array Arguments ..
16 LOGICAL DOTYPE( * )
17 INTEGER ISEED( 4 ), NN( * )
18 DOUBLE PRECISION RESULT( * ), RWORK( * )
19 COMPLEX*16 A( LDA, * ), ALPHA( * ), ALPHA1( * ),
20 $ B( LDA, * ), BETA( * ), BETA1( * ),
21 $ Q( LDQ, * ), QE( LDQE, * ), S( LDA, * ),
22 $ T( LDA, * ), WORK( * ), Z( LDQ, * )
23 * ..
24 *
25 * Purpose
26 * =======
27 *
28 * ZDRGEV checks the nonsymmetric generalized eigenvalue problem driver
29 * routine ZGGEV.
30 *
31 * ZGGEV computes for a pair of n-by-n nonsymmetric matrices (A,B) the
32 * generalized eigenvalues and, optionally, the left and right
33 * eigenvectors.
34 *
35 * A generalized eigenvalue for a pair of matrices (A,B) is a scalar w
36 * or a ratio alpha/beta = w, such that A - w*B is singular. It is
37 * usually represented as the pair (alpha,beta), as there is reasonalbe
38 * interpretation for beta=0, and even for both being zero.
39 *
40 * A right generalized eigenvector corresponding to a generalized
41 * eigenvalue w for a pair of matrices (A,B) is a vector r such that
42 * (A - wB) * r = 0. A left generalized eigenvector is a vector l such
43 * that l**H * (A - wB) = 0, where l**H is the conjugate-transpose of l.
44 *
45 * When ZDRGEV is called, a number of matrix "sizes" ("n's") and a
46 * number of matrix "types" are specified. For each size ("n")
47 * and each type of matrix, a pair of matrices (A, B) will be generated
48 * and used for testing. For each matrix pair, the following tests
49 * will be performed and compared with the threshhold THRESH.
50 *
51 * Results from ZGGEV:
52 *
53 * (1) max over all left eigenvalue/-vector pairs (alpha/beta,l) of
54 *
55 * | VL**H * (beta A - alpha B) |/( ulp max(|beta A|, |alpha B|) )
56 *
57 * where VL**H is the conjugate-transpose of VL.
58 *
59 * (2) | |VL(i)| - 1 | / ulp and whether largest component real
60 *
61 * VL(i) denotes the i-th column of VL.
62 *
63 * (3) max over all left eigenvalue/-vector pairs (alpha/beta,r) of
64 *
65 * | (beta A - alpha B) * VR | / ( ulp max(|beta A|, |alpha B|) )
66 *
67 * (4) | |VR(i)| - 1 | / ulp and whether largest component real
68 *
69 * VR(i) denotes the i-th column of VR.
70 *
71 * (5) W(full) = W(partial)
72 * W(full) denotes the eigenvalues computed when both l and r
73 * are also computed, and W(partial) denotes the eigenvalues
74 * computed when only W, only W and r, or only W and l are
75 * computed.
76 *
77 * (6) VL(full) = VL(partial)
78 * VL(full) denotes the left eigenvectors computed when both l
79 * and r are computed, and VL(partial) denotes the result
80 * when only l is computed.
81 *
82 * (7) VR(full) = VR(partial)
83 * VR(full) denotes the right eigenvectors computed when both l
84 * and r are also computed, and VR(partial) denotes the result
85 * when only l is computed.
86 *
87 *
88 * Test Matrices
89 * ---- --------
90 *
91 * The sizes of the test matrices are specified by an array
92 * NN(1:NSIZES); the value of each element NN(j) specifies one size.
93 * The "types" are specified by a logical array DOTYPE( 1:NTYPES ); if
94 * DOTYPE(j) is .TRUE., then matrix type "j" will be generated.
95 * Currently, the list of possible types is:
96 *
97 * (1) ( 0, 0 ) (a pair of zero matrices)
98 *
99 * (2) ( I, 0 ) (an identity and a zero matrix)
100 *
101 * (3) ( 0, I ) (an identity and a zero matrix)
102 *
103 * (4) ( I, I ) (a pair of identity matrices)
104 *
105 * t t
106 * (5) ( J , J ) (a pair of transposed Jordan blocks)
107 *
108 * t ( I 0 )
109 * (6) ( X, Y ) where X = ( J 0 ) and Y = ( t )
110 * ( 0 I ) ( 0 J )
111 * and I is a k x k identity and J a (k+1)x(k+1)
112 * Jordan block; k=(N-1)/2
113 *
114 * (7) ( D, I ) where D is diag( 0, 1,..., N-1 ) (a diagonal
115 * matrix with those diagonal entries.)
116 * (8) ( I, D )
117 *
118 * (9) ( big*D, small*I ) where "big" is near overflow and small=1/big
119 *
120 * (10) ( small*D, big*I )
121 *
122 * (11) ( big*I, small*D )
123 *
124 * (12) ( small*I, big*D )
125 *
126 * (13) ( big*D, big*I )
127 *
128 * (14) ( small*D, small*I )
129 *
130 * (15) ( D1, D2 ) where D1 is diag( 0, 0, 1, ..., N-3, 0 ) and
131 * D2 is diag( 0, N-3, N-4,..., 1, 0, 0 )
132 * t t
133 * (16) Q ( J , J ) Z where Q and Z are random orthogonal matrices.
134 *
135 * (17) Q ( T1, T2 ) Z where T1 and T2 are upper triangular matrices
136 * with random O(1) entries above the diagonal
137 * and diagonal entries diag(T1) =
138 * ( 0, 0, 1, ..., N-3, 0 ) and diag(T2) =
139 * ( 0, N-3, N-4,..., 1, 0, 0 )
140 *
141 * (18) Q ( T1, T2 ) Z diag(T1) = ( 0, 0, 1, 1, s, ..., s, 0 )
142 * diag(T2) = ( 0, 1, 0, 1,..., 1, 0 )
143 * s = machine precision.
144 *
145 * (19) Q ( T1, T2 ) Z diag(T1)=( 0,0,1,1, 1-d, ..., 1-(N-5)*d=s, 0 )
146 * diag(T2) = ( 0, 1, 0, 1, ..., 1, 0 )
147 *
148 * N-5
149 * (20) Q ( T1, T2 ) Z diag(T1)=( 0, 0, 1, 1, a, ..., a =s, 0 )
150 * diag(T2) = ( 0, 1, 0, 1, ..., 1, 0, 0 )
151 *
152 * (21) Q ( T1, T2 ) Z diag(T1)=( 0, 0, 1, r1, r2, ..., r(N-4), 0 )
153 * diag(T2) = ( 0, 1, 0, 1, ..., 1, 0, 0 )
154 * where r1,..., r(N-4) are random.
155 *
156 * (22) Q ( big*T1, small*T2 ) Z diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
157 * diag(T2) = ( 0, 1, ..., 1, 0, 0 )
158 *
159 * (23) Q ( small*T1, big*T2 ) Z diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
160 * diag(T2) = ( 0, 1, ..., 1, 0, 0 )
161 *
162 * (24) Q ( small*T1, small*T2 ) Z diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
163 * diag(T2) = ( 0, 1, ..., 1, 0, 0 )
164 *
165 * (25) Q ( big*T1, big*T2 ) Z diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
166 * diag(T2) = ( 0, 1, ..., 1, 0, 0 )
167 *
168 * (26) Q ( T1, T2 ) Z where T1 and T2 are random upper-triangular
169 * matrices.
170 *
171 *
172 * Arguments
173 * =========
174 *
175 * NSIZES (input) INTEGER
176 * The number of sizes of matrices to use. If it is zero,
177 * ZDRGES does nothing. NSIZES >= 0.
178 *
179 * NN (input) INTEGER array, dimension (NSIZES)
180 * An array containing the sizes to be used for the matrices.
181 * Zero values will be skipped. NN >= 0.
182 *
183 * NTYPES (input) INTEGER
184 * The number of elements in DOTYPE. If it is zero, ZDRGEV
185 * does nothing. It must be at least zero. If it is MAXTYP+1
186 * and NSIZES is 1, then an additional type, MAXTYP+1 is
187 * defined, which is to use whatever matrix is in A. This
188 * is only useful if DOTYPE(1:MAXTYP) is .FALSE. and
189 * DOTYPE(MAXTYP+1) is .TRUE. .
190 *
191 * DOTYPE (input) LOGICAL array, dimension (NTYPES)
192 * If DOTYPE(j) is .TRUE., then for each size in NN a
193 * matrix of that size and of type j will be generated.
194 * If NTYPES is smaller than the maximum number of types
195 * defined (PARAMETER MAXTYP), then types NTYPES+1 through
196 * MAXTYP will not be generated. If NTYPES is larger
197 * than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES)
198 * will be ignored.
199 *
200 * ISEED (input/output) INTEGER array, dimension (4)
201 * On entry ISEED specifies the seed of the random number
202 * generator. The array elements should be between 0 and 4095;
203 * if not they will be reduced mod 4096. Also, ISEED(4) must
204 * be odd. The random number generator uses a linear
205 * congruential sequence limited to small integers, and so
206 * should produce machine independent random numbers. The
207 * values of ISEED are changed on exit, and can be used in the
208 * next call to ZDRGES to continue the same random number
209 * sequence.
210 *
211 * THRESH (input) DOUBLE PRECISION
212 * A test will count as "failed" if the "error", computed as
213 * described above, exceeds THRESH. Note that the error is
214 * scaled to be O(1), so THRESH should be a reasonably small
215 * multiple of 1, e.g., 10 or 100. In particular, it should
216 * not depend on the precision (single vs. double) or the size
217 * of the matrix. It must be at least zero.
218 *
219 * NOUNIT (input) INTEGER
220 * The FORTRAN unit number for printing out error messages
221 * (e.g., if a routine returns IERR not equal to 0.)
222 *
223 * A (input/workspace) COMPLEX*16 array, dimension(LDA, max(NN))
224 * Used to hold the original A matrix. Used as input only
225 * if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and
226 * DOTYPE(MAXTYP+1)=.TRUE.
227 *
228 * LDA (input) INTEGER
229 * The leading dimension of A, B, S, and T.
230 * It must be at least 1 and at least max( NN ).
231 *
232 * B (input/workspace) COMPLEX*16 array, dimension(LDA, max(NN))
233 * Used to hold the original B matrix. Used as input only
234 * if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and
235 * DOTYPE(MAXTYP+1)=.TRUE.
236 *
237 * S (workspace) COMPLEX*16 array, dimension (LDA, max(NN))
238 * The Schur form matrix computed from A by ZGGEV. On exit, S
239 * contains the Schur form matrix corresponding to the matrix
240 * in A.
241 *
242 * T (workspace) COMPLEX*16 array, dimension (LDA, max(NN))
243 * The upper triangular matrix computed from B by ZGGEV.
244 *
245 * Q (workspace) COMPLEX*16 array, dimension (LDQ, max(NN))
246 * The (left) eigenvectors matrix computed by ZGGEV.
247 *
248 * LDQ (input) INTEGER
249 * The leading dimension of Q and Z. It must
250 * be at least 1 and at least max( NN ).
251 *
252 * Z (workspace) COMPLEX*16 array, dimension( LDQ, max(NN) )
253 * The (right) orthogonal matrix computed by ZGGEV.
254 *
255 * QE (workspace) COMPLEX*16 array, dimension( LDQ, max(NN) )
256 * QE holds the computed right or left eigenvectors.
257 *
258 * LDQE (input) INTEGER
259 * The leading dimension of QE. LDQE >= max(1,max(NN)).
260 *
261 * ALPHA (workspace) COMPLEX*16 array, dimension (max(NN))
262 * BETA (workspace) COMPLEX*16 array, dimension (max(NN))
263 * The generalized eigenvalues of (A,B) computed by ZGGEV.
264 * ( ALPHAR(k)+ALPHAI(k)*i ) / BETA(k) is the k-th
265 * generalized eigenvalue of A and B.
266 *
267 * ALPHA1 (workspace) COMPLEX*16 array, dimension (max(NN))
268 * BETA1 (workspace) COMPLEX*16 array, dimension (max(NN))
269 * Like ALPHAR, ALPHAI, BETA, these arrays contain the
270 * eigenvalues of A and B, but those computed when ZGGEV only
271 * computes a partial eigendecomposition, i.e. not the
272 * eigenvalues and left and right eigenvectors.
273 *
274 * WORK (workspace) COMPLEX*16 array, dimension (LWORK)
275 *
276 * LWORK (input) INTEGER
277 * The number of entries in WORK. LWORK >= N*(N+1)
278 *
279 * RWORK (workspace) DOUBLE PRECISION array, dimension (8*N)
280 * Real workspace.
281 *
282 * RESULT (output) DOUBLE PRECISION array, dimension (2)
283 * The values computed by the tests described above.
284 * The values are currently limited to 1/ulp, to avoid overflow.
285 *
286 * INFO (output) INTEGER
287 * = 0: successful exit
288 * < 0: if INFO = -i, the i-th argument had an illegal value.
289 * > 0: A routine returned an error code. INFO is the
290 * absolute value of the INFO value returned.
291 *
292 * =====================================================================
293 *
294 * .. Parameters ..
295 DOUBLE PRECISION ZERO, ONE
296 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
297 COMPLEX*16 CZERO, CONE
298 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
299 $ CONE = ( 1.0D+0, 0.0D+0 ) )
300 INTEGER MAXTYP
301 PARAMETER ( MAXTYP = 26 )
302 * ..
303 * .. Local Scalars ..
304 LOGICAL BADNN
305 INTEGER I, IADD, IERR, IN, J, JC, JR, JSIZE, JTYPE,
306 $ MAXWRK, MINWRK, MTYPES, N, N1, NB, NERRS,
307 $ NMATS, NMAX, NTESTT
308 DOUBLE PRECISION SAFMAX, SAFMIN, ULP, ULPINV
309 COMPLEX*16 CTEMP
310 * ..
311 * .. Local Arrays ..
312 LOGICAL LASIGN( MAXTYP ), LBSIGN( MAXTYP )
313 INTEGER IOLDSD( 4 ), KADD( 6 ), KAMAGN( MAXTYP ),
314 $ KATYPE( MAXTYP ), KAZERO( MAXTYP ),
315 $ KBMAGN( MAXTYP ), KBTYPE( MAXTYP ),
316 $ KBZERO( MAXTYP ), KCLASS( MAXTYP ),
317 $ KTRIAN( MAXTYP ), KZ1( 6 ), KZ2( 6 )
318 DOUBLE PRECISION RMAGN( 0: 3 )
319 * ..
320 * .. External Functions ..
321 INTEGER ILAENV
322 DOUBLE PRECISION DLAMCH
323 COMPLEX*16 ZLARND
324 EXTERNAL ILAENV, DLAMCH, ZLARND
325 * ..
326 * .. External Subroutines ..
327 EXTERNAL ALASVM, DLABAD, XERBLA, ZGET52, ZGGEV, ZLACPY,
328 $ ZLARFG, ZLASET, ZLATM4, ZUNM2R
329 * ..
330 * .. Intrinsic Functions ..
331 INTRINSIC ABS, DBLE, DCONJG, MAX, MIN, SIGN
332 * ..
333 * .. Data statements ..
334 DATA KCLASS / 15*1, 10*2, 1*3 /
335 DATA KZ1 / 0, 1, 2, 1, 3, 3 /
336 DATA KZ2 / 0, 0, 1, 2, 1, 1 /
337 DATA KADD / 0, 0, 0, 0, 3, 2 /
338 DATA KATYPE / 0, 1, 0, 1, 2, 3, 4, 1, 4, 4, 1, 1, 4,
339 $ 4, 4, 2, 4, 5, 8, 7, 9, 4*4, 0 /
340 DATA KBTYPE / 0, 0, 1, 1, 2, -3, 1, 4, 1, 1, 4, 4,
341 $ 1, 1, -4, 2, -4, 8*8, 0 /
342 DATA KAZERO / 6*1, 2, 1, 2*2, 2*1, 2*2, 3, 1, 3,
343 $ 4*5, 4*3, 1 /
344 DATA KBZERO / 6*1, 1, 2, 2*1, 2*2, 2*1, 4, 1, 4,
345 $ 4*6, 4*4, 1 /
346 DATA KAMAGN / 8*1, 2, 3, 2, 3, 2, 3, 7*1, 2, 3, 3,
347 $ 2, 1 /
348 DATA KBMAGN / 8*1, 3, 2, 3, 2, 2, 3, 7*1, 3, 2, 3,
349 $ 2, 1 /
350 DATA KTRIAN / 16*0, 10*1 /
351 DATA LASIGN / 6*.FALSE., .TRUE., .FALSE., 2*.TRUE.,
352 $ 2*.FALSE., 3*.TRUE., .FALSE., .TRUE.,
353 $ 3*.FALSE., 5*.TRUE., .FALSE. /
354 DATA LBSIGN / 7*.FALSE., .TRUE., 2*.FALSE.,
355 $ 2*.TRUE., 2*.FALSE., .TRUE., .FALSE., .TRUE.,
356 $ 9*.FALSE. /
357 * ..
358 * .. Executable Statements ..
359 *
360 * Check for errors
361 *
362 INFO = 0
363 *
364 BADNN = .FALSE.
365 NMAX = 1
366 DO 10 J = 1, NSIZES
367 NMAX = MAX( NMAX, NN( J ) )
368 IF( NN( J ).LT.0 )
369 $ BADNN = .TRUE.
370 10 CONTINUE
371 *
372 IF( NSIZES.LT.0 ) THEN
373 INFO = -1
374 ELSE IF( BADNN ) THEN
375 INFO = -2
376 ELSE IF( NTYPES.LT.0 ) THEN
377 INFO = -3
378 ELSE IF( THRESH.LT.ZERO ) THEN
379 INFO = -6
380 ELSE IF( LDA.LE.1 .OR. LDA.LT.NMAX ) THEN
381 INFO = -9
382 ELSE IF( LDQ.LE.1 .OR. LDQ.LT.NMAX ) THEN
383 INFO = -14
384 ELSE IF( LDQE.LE.1 .OR. LDQE.LT.NMAX ) THEN
385 INFO = -17
386 END IF
387 *
388 * Compute workspace
389 * (Note: Comments in the code beginning "Workspace:" describe the
390 * minimal amount of workspace needed at that point in the code,
391 * as well as the preferred amount for good performance.
392 * NB refers to the optimal block size for the immediately
393 * following subroutine, as returned by ILAENV.
394 *
395 MINWRK = 1
396 IF( INFO.EQ.0 .AND. LWORK.GE.1 ) THEN
397 MINWRK = NMAX*( NMAX+1 )
398 NB = MAX( 1, ILAENV( 1, 'ZGEQRF', ' ', NMAX, NMAX, -1, -1 ),
399 $ ILAENV( 1, 'ZUNMQR', 'LC', NMAX, NMAX, NMAX, -1 ),
400 $ ILAENV( 1, 'ZUNGQR', ' ', NMAX, NMAX, NMAX, -1 ) )
401 MAXWRK = MAX( 2*NMAX, NMAX*( NB+1 ), NMAX*( NMAX+1 ) )
402 WORK( 1 ) = MAXWRK
403 END IF
404 *
405 IF( LWORK.LT.MINWRK )
406 $ INFO = -23
407 *
408 IF( INFO.NE.0 ) THEN
409 CALL XERBLA( 'ZDRGEV', -INFO )
410 RETURN
411 END IF
412 *
413 * Quick return if possible
414 *
415 IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 )
416 $ RETURN
417 *
418 ULP = DLAMCH( 'Precision' )
419 SAFMIN = DLAMCH( 'Safe minimum' )
420 SAFMIN = SAFMIN / ULP
421 SAFMAX = ONE / SAFMIN
422 CALL DLABAD( SAFMIN, SAFMAX )
423 ULPINV = ONE / ULP
424 *
425 * The values RMAGN(2:3) depend on N, see below.
426 *
427 RMAGN( 0 ) = ZERO
428 RMAGN( 1 ) = ONE
429 *
430 * Loop over sizes, types
431 *
432 NTESTT = 0
433 NERRS = 0
434 NMATS = 0
435 *
436 DO 220 JSIZE = 1, NSIZES
437 N = NN( JSIZE )
438 N1 = MAX( 1, N )
439 RMAGN( 2 ) = SAFMAX*ULP / DBLE( N1 )
440 RMAGN( 3 ) = SAFMIN*ULPINV*N1
441 *
442 IF( NSIZES.NE.1 ) THEN
443 MTYPES = MIN( MAXTYP, NTYPES )
444 ELSE
445 MTYPES = MIN( MAXTYP+1, NTYPES )
446 END IF
447 *
448 DO 210 JTYPE = 1, MTYPES
449 IF( .NOT.DOTYPE( JTYPE ) )
450 $ GO TO 210
451 NMATS = NMATS + 1
452 *
453 * Save ISEED in case of an error.
454 *
455 DO 20 J = 1, 4
456 IOLDSD( J ) = ISEED( J )
457 20 CONTINUE
458 *
459 * Generate test matrices A and B
460 *
461 * Description of control parameters:
462 *
463 * KZLASS: =1 means w/o rotation, =2 means w/ rotation,
464 * =3 means random.
465 * KATYPE: the "type" to be passed to ZLATM4 for computing A.
466 * KAZERO: the pattern of zeros on the diagonal for A:
467 * =1: ( xxx ), =2: (0, xxx ) =3: ( 0, 0, xxx, 0 ),
468 * =4: ( 0, xxx, 0, 0 ), =5: ( 0, 0, 1, xxx, 0 ),
469 * =6: ( 0, 1, 0, xxx, 0 ). (xxx means a string of
470 * non-zero entries.)
471 * KAMAGN: the magnitude of the matrix: =0: zero, =1: O(1),
472 * =2: large, =3: small.
473 * LASIGN: .TRUE. if the diagonal elements of A are to be
474 * multiplied by a random magnitude 1 number.
475 * KBTYPE, KBZERO, KBMAGN, LBSIGN: the same, but for B.
476 * KTRIAN: =0: don't fill in the upper triangle, =1: do.
477 * KZ1, KZ2, KADD: used to implement KAZERO and KBZERO.
478 * RMAGN: used to implement KAMAGN and KBMAGN.
479 *
480 IF( MTYPES.GT.MAXTYP )
481 $ GO TO 100
482 IERR = 0
483 IF( KCLASS( JTYPE ).LT.3 ) THEN
484 *
485 * Generate A (w/o rotation)
486 *
487 IF( ABS( KATYPE( JTYPE ) ).EQ.3 ) THEN
488 IN = 2*( ( N-1 ) / 2 ) + 1
489 IF( IN.NE.N )
490 $ CALL ZLASET( 'Full', N, N, CZERO, CZERO, A, LDA )
491 ELSE
492 IN = N
493 END IF
494 CALL ZLATM4( KATYPE( JTYPE ), IN, KZ1( KAZERO( JTYPE ) ),
495 $ KZ2( KAZERO( JTYPE ) ), LASIGN( JTYPE ),
496 $ RMAGN( KAMAGN( JTYPE ) ), ULP,
497 $ RMAGN( KTRIAN( JTYPE )*KAMAGN( JTYPE ) ), 2,
498 $ ISEED, A, LDA )
499 IADD = KADD( KAZERO( JTYPE ) )
500 IF( IADD.GT.0 .AND. IADD.LE.N )
501 $ A( IADD, IADD ) = RMAGN( KAMAGN( JTYPE ) )
502 *
503 * Generate B (w/o rotation)
504 *
505 IF( ABS( KBTYPE( JTYPE ) ).EQ.3 ) THEN
506 IN = 2*( ( N-1 ) / 2 ) + 1
507 IF( IN.NE.N )
508 $ CALL ZLASET( 'Full', N, N, CZERO, CZERO, B, LDA )
509 ELSE
510 IN = N
511 END IF
512 CALL ZLATM4( KBTYPE( JTYPE ), IN, KZ1( KBZERO( JTYPE ) ),
513 $ KZ2( KBZERO( JTYPE ) ), LBSIGN( JTYPE ),
514 $ RMAGN( KBMAGN( JTYPE ) ), ONE,
515 $ RMAGN( KTRIAN( JTYPE )*KBMAGN( JTYPE ) ), 2,
516 $ ISEED, B, LDA )
517 IADD = KADD( KBZERO( JTYPE ) )
518 IF( IADD.NE.0 .AND. IADD.LE.N )
519 $ B( IADD, IADD ) = RMAGN( KBMAGN( JTYPE ) )
520 *
521 IF( KCLASS( JTYPE ).EQ.2 .AND. N.GT.0 ) THEN
522 *
523 * Include rotations
524 *
525 * Generate Q, Z as Householder transformations times
526 * a diagonal matrix.
527 *
528 DO 40 JC = 1, N - 1
529 DO 30 JR = JC, N
530 Q( JR, JC ) = ZLARND( 3, ISEED )
531 Z( JR, JC ) = ZLARND( 3, ISEED )
532 30 CONTINUE
533 CALL ZLARFG( N+1-JC, Q( JC, JC ), Q( JC+1, JC ), 1,
534 $ WORK( JC ) )
535 WORK( 2*N+JC ) = SIGN( ONE, DBLE( Q( JC, JC ) ) )
536 Q( JC, JC ) = CONE
537 CALL ZLARFG( N+1-JC, Z( JC, JC ), Z( JC+1, JC ), 1,
538 $ WORK( N+JC ) )
539 WORK( 3*N+JC ) = SIGN( ONE, DBLE( Z( JC, JC ) ) )
540 Z( JC, JC ) = CONE
541 40 CONTINUE
542 CTEMP = ZLARND( 3, ISEED )
543 Q( N, N ) = CONE
544 WORK( N ) = CZERO
545 WORK( 3*N ) = CTEMP / ABS( CTEMP )
546 CTEMP = ZLARND( 3, ISEED )
547 Z( N, N ) = CONE
548 WORK( 2*N ) = CZERO
549 WORK( 4*N ) = CTEMP / ABS( CTEMP )
550 *
551 * Apply the diagonal matrices
552 *
553 DO 60 JC = 1, N
554 DO 50 JR = 1, N
555 A( JR, JC ) = WORK( 2*N+JR )*
556 $ DCONJG( WORK( 3*N+JC ) )*
557 $ A( JR, JC )
558 B( JR, JC ) = WORK( 2*N+JR )*
559 $ DCONJG( WORK( 3*N+JC ) )*
560 $ B( JR, JC )
561 50 CONTINUE
562 60 CONTINUE
563 CALL ZUNM2R( 'L', 'N', N, N, N-1, Q, LDQ, WORK, A,
564 $ LDA, WORK( 2*N+1 ), IERR )
565 IF( IERR.NE.0 )
566 $ GO TO 90
567 CALL ZUNM2R( 'R', 'C', N, N, N-1, Z, LDQ, WORK( N+1 ),
568 $ A, LDA, WORK( 2*N+1 ), IERR )
569 IF( IERR.NE.0 )
570 $ GO TO 90
571 CALL ZUNM2R( 'L', 'N', N, N, N-1, Q, LDQ, WORK, B,
572 $ LDA, WORK( 2*N+1 ), IERR )
573 IF( IERR.NE.0 )
574 $ GO TO 90
575 CALL ZUNM2R( 'R', 'C', N, N, N-1, Z, LDQ, WORK( N+1 ),
576 $ B, LDA, WORK( 2*N+1 ), IERR )
577 IF( IERR.NE.0 )
578 $ GO TO 90
579 END IF
580 ELSE
581 *
582 * Random matrices
583 *
584 DO 80 JC = 1, N
585 DO 70 JR = 1, N
586 A( JR, JC ) = RMAGN( KAMAGN( JTYPE ) )*
587 $ ZLARND( 4, ISEED )
588 B( JR, JC ) = RMAGN( KBMAGN( JTYPE ) )*
589 $ ZLARND( 4, ISEED )
590 70 CONTINUE
591 80 CONTINUE
592 END IF
593 *
594 90 CONTINUE
595 *
596 IF( IERR.NE.0 ) THEN
597 WRITE( NOUNIT, FMT = 9999 )'Generator', IERR, N, JTYPE,
598 $ IOLDSD
599 INFO = ABS( IERR )
600 RETURN
601 END IF
602 *
603 100 CONTINUE
604 *
605 DO 110 I = 1, 7
606 RESULT( I ) = -ONE
607 110 CONTINUE
608 *
609 * Call ZGGEV to compute eigenvalues and eigenvectors.
610 *
611 CALL ZLACPY( ' ', N, N, A, LDA, S, LDA )
612 CALL ZLACPY( ' ', N, N, B, LDA, T, LDA )
613 CALL ZGGEV( 'V', 'V', N, S, LDA, T, LDA, ALPHA, BETA, Q,
614 $ LDQ, Z, LDQ, WORK, LWORK, RWORK, IERR )
615 IF( IERR.NE.0 .AND. IERR.NE.N+1 ) THEN
616 RESULT( 1 ) = ULPINV
617 WRITE( NOUNIT, FMT = 9999 )'ZGGEV1', IERR, N, JTYPE,
618 $ IOLDSD
619 INFO = ABS( IERR )
620 GO TO 190
621 END IF
622 *
623 * Do the tests (1) and (2)
624 *
625 CALL ZGET52( .TRUE., N, A, LDA, B, LDA, Q, LDQ, ALPHA, BETA,
626 $ WORK, RWORK, RESULT( 1 ) )
627 IF( RESULT( 2 ).GT.THRESH ) THEN
628 WRITE( NOUNIT, FMT = 9998 )'Left', 'ZGGEV1',
629 $ RESULT( 2 ), N, JTYPE, IOLDSD
630 END IF
631 *
632 * Do the tests (3) and (4)
633 *
634 CALL ZGET52( .FALSE., N, A, LDA, B, LDA, Z, LDQ, ALPHA,
635 $ BETA, WORK, RWORK, RESULT( 3 ) )
636 IF( RESULT( 4 ).GT.THRESH ) THEN
637 WRITE( NOUNIT, FMT = 9998 )'Right', 'ZGGEV1',
638 $ RESULT( 4 ), N, JTYPE, IOLDSD
639 END IF
640 *
641 * Do test (5)
642 *
643 CALL ZLACPY( ' ', N, N, A, LDA, S, LDA )
644 CALL ZLACPY( ' ', N, N, B, LDA, T, LDA )
645 CALL ZGGEV( 'N', 'N', N, S, LDA, T, LDA, ALPHA1, BETA1, Q,
646 $ LDQ, Z, LDQ, WORK, LWORK, RWORK, IERR )
647 IF( IERR.NE.0 .AND. IERR.NE.N+1 ) THEN
648 RESULT( 1 ) = ULPINV
649 WRITE( NOUNIT, FMT = 9999 )'ZGGEV2', IERR, N, JTYPE,
650 $ IOLDSD
651 INFO = ABS( IERR )
652 GO TO 190
653 END IF
654 *
655 DO 120 J = 1, N
656 IF( ALPHA( J ).NE.ALPHA1( J ) .OR. BETA( J ).NE.
657 $ BETA1( J ) )RESULT( 5 ) = ULPINV
658 120 CONTINUE
659 *
660 * Do test (6): Compute eigenvalues and left eigenvectors,
661 * and test them
662 *
663 CALL ZLACPY( ' ', N, N, A, LDA, S, LDA )
664 CALL ZLACPY( ' ', N, N, B, LDA, T, LDA )
665 CALL ZGGEV( 'V', 'N', N, S, LDA, T, LDA, ALPHA1, BETA1, QE,
666 $ LDQE, Z, LDQ, WORK, LWORK, RWORK, IERR )
667 IF( IERR.NE.0 .AND. IERR.NE.N+1 ) THEN
668 RESULT( 1 ) = ULPINV
669 WRITE( NOUNIT, FMT = 9999 )'ZGGEV3', IERR, N, JTYPE,
670 $ IOLDSD
671 INFO = ABS( IERR )
672 GO TO 190
673 END IF
674 *
675 DO 130 J = 1, N
676 IF( ALPHA( J ).NE.ALPHA1( J ) .OR. BETA( J ).NE.
677 $ BETA1( J ) )RESULT( 6 ) = ULPINV
678 130 CONTINUE
679 *
680 DO 150 J = 1, N
681 DO 140 JC = 1, N
682 IF( Q( J, JC ).NE.QE( J, JC ) )
683 $ RESULT( 6 ) = ULPINV
684 140 CONTINUE
685 150 CONTINUE
686 *
687 * Do test (7): Compute eigenvalues and right eigenvectors,
688 * and test them
689 *
690 CALL ZLACPY( ' ', N, N, A, LDA, S, LDA )
691 CALL ZLACPY( ' ', N, N, B, LDA, T, LDA )
692 CALL ZGGEV( 'N', 'V', N, S, LDA, T, LDA, ALPHA1, BETA1, Q,
693 $ LDQ, QE, LDQE, WORK, LWORK, RWORK, IERR )
694 IF( IERR.NE.0 .AND. IERR.NE.N+1 ) THEN
695 RESULT( 1 ) = ULPINV
696 WRITE( NOUNIT, FMT = 9999 )'ZGGEV4', IERR, N, JTYPE,
697 $ IOLDSD
698 INFO = ABS( IERR )
699 GO TO 190
700 END IF
701 *
702 DO 160 J = 1, N
703 IF( ALPHA( J ).NE.ALPHA1( J ) .OR. BETA( J ).NE.
704 $ BETA1( J ) )RESULT( 7 ) = ULPINV
705 160 CONTINUE
706 *
707 DO 180 J = 1, N
708 DO 170 JC = 1, N
709 IF( Z( J, JC ).NE.QE( J, JC ) )
710 $ RESULT( 7 ) = ULPINV
711 170 CONTINUE
712 180 CONTINUE
713 *
714 * End of Loop -- Check for RESULT(j) > THRESH
715 *
716 190 CONTINUE
717 *
718 NTESTT = NTESTT + 7
719 *
720 * Print out tests which fail.
721 *
722 DO 200 JR = 1, 7
723 IF( RESULT( JR ).GE.THRESH ) THEN
724 *
725 * If this is the first test to fail,
726 * print a header to the data file.
727 *
728 IF( NERRS.EQ.0 ) THEN
729 WRITE( NOUNIT, FMT = 9997 )'ZGV'
730 *
731 * Matrix types
732 *
733 WRITE( NOUNIT, FMT = 9996 )
734 WRITE( NOUNIT, FMT = 9995 )
735 WRITE( NOUNIT, FMT = 9994 )'Orthogonal'
736 *
737 * Tests performed
738 *
739 WRITE( NOUNIT, FMT = 9993 )
740 *
741 END IF
742 NERRS = NERRS + 1
743 IF( RESULT( JR ).LT.10000.0D0 ) THEN
744 WRITE( NOUNIT, FMT = 9992 )N, JTYPE, IOLDSD, JR,
745 $ RESULT( JR )
746 ELSE
747 WRITE( NOUNIT, FMT = 9991 )N, JTYPE, IOLDSD, JR,
748 $ RESULT( JR )
749 END IF
750 END IF
751 200 CONTINUE
752 *
753 210 CONTINUE
754 220 CONTINUE
755 *
756 * Summary
757 *
758 CALL ALASVM( 'ZGV', NOUNIT, NERRS, NTESTT, 0 )
759 *
760 WORK( 1 ) = MAXWRK
761 *
762 RETURN
763 *
764 9999 FORMAT( ' ZDRGEV: ', A, ' returned INFO=', I6, '.', / 3X, 'N=',
765 $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' )
766 *
767 9998 FORMAT( ' ZDRGEV: ', A, ' Eigenvectors from ', A, ' incorrectly ',
768 $ 'normalized.', / ' Bits of error=', 0P, G10.3, ',', 3X,
769 $ 'N=', I4, ', JTYPE=', I3, ', ISEED=(', 3( I4, ',' ), I5,
770 $ ')' )
771 *
772 9997 FORMAT( / 1X, A3, ' -- Complex Generalized eigenvalue problem ',
773 $ 'driver' )
774 *
775 9996 FORMAT( ' Matrix types (see ZDRGEV for details): ' )
776 *
777 9995 FORMAT( ' Special Matrices:', 23X,
778 $ '(J''=transposed Jordan block)',
779 $ / ' 1=(0,0) 2=(I,0) 3=(0,I) 4=(I,I) 5=(J'',J'') ',
780 $ '6=(diag(J'',I), diag(I,J''))', / ' Diagonal Matrices: ( ',
781 $ 'D=diag(0,1,2,...) )', / ' 7=(D,I) 9=(large*D, small*I',
782 $ ') 11=(large*I, small*D) 13=(large*D, large*I)', /
783 $ ' 8=(I,D) 10=(small*D, large*I) 12=(small*I, large*D) ',
784 $ ' 14=(small*D, small*I)', / ' 15=(D, reversed D)' )
785 9994 FORMAT( ' Matrices Rotated by Random ', A, ' Matrices U, V:',
786 $ / ' 16=Transposed Jordan Blocks 19=geometric ',
787 $ 'alpha, beta=0,1', / ' 17=arithm. alpha&beta ',
788 $ ' 20=arithmetic alpha, beta=0,1', / ' 18=clustered ',
789 $ 'alpha, beta=0,1 21=random alpha, beta=0,1',
790 $ / ' Large & Small Matrices:', / ' 22=(large, small) ',
791 $ '23=(small,large) 24=(small,small) 25=(large,large)',
792 $ / ' 26=random O(1) matrices.' )
793 *
794 9993 FORMAT( / ' Tests performed: ',
795 $ / ' 1 = max | ( b A - a B )''*l | / const.,',
796 $ / ' 2 = | |VR(i)| - 1 | / ulp,',
797 $ / ' 3 = max | ( b A - a B )*r | / const.',
798 $ / ' 4 = | |VL(i)| - 1 | / ulp,',
799 $ / ' 5 = 0 if W same no matter if r or l computed,',
800 $ / ' 6 = 0 if l same no matter if l computed,',
801 $ / ' 7 = 0 if r same no matter if r computed,', / 1X )
802 9992 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=',
803 $ 4( I4, ',' ), ' result ', I2, ' is', 0P, F8.2 )
804 9991 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=',
805 $ 4( I4, ',' ), ' result ', I2, ' is', 1P, D10.3 )
806 *
807 * End of ZDRGEV
808 *
809 END
2 $ NOUNIT, A, LDA, B, S, T, Q, LDQ, Z, QE, LDQE,
3 $ ALPHA, BETA, ALPHA1, BETA1, WORK, LWORK, RWORK,
4 $ RESULT, INFO )
5 *
6 * -- LAPACK test routine (version 3.1.1) --
7 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
8 * February 2007
9 *
10 * .. Scalar Arguments ..
11 INTEGER INFO, LDA, LDQ, LDQE, LWORK, NOUNIT, NSIZES,
12 $ NTYPES
13 DOUBLE PRECISION THRESH
14 * ..
15 * .. Array Arguments ..
16 LOGICAL DOTYPE( * )
17 INTEGER ISEED( 4 ), NN( * )
18 DOUBLE PRECISION RESULT( * ), RWORK( * )
19 COMPLEX*16 A( LDA, * ), ALPHA( * ), ALPHA1( * ),
20 $ B( LDA, * ), BETA( * ), BETA1( * ),
21 $ Q( LDQ, * ), QE( LDQE, * ), S( LDA, * ),
22 $ T( LDA, * ), WORK( * ), Z( LDQ, * )
23 * ..
24 *
25 * Purpose
26 * =======
27 *
28 * ZDRGEV checks the nonsymmetric generalized eigenvalue problem driver
29 * routine ZGGEV.
30 *
31 * ZGGEV computes for a pair of n-by-n nonsymmetric matrices (A,B) the
32 * generalized eigenvalues and, optionally, the left and right
33 * eigenvectors.
34 *
35 * A generalized eigenvalue for a pair of matrices (A,B) is a scalar w
36 * or a ratio alpha/beta = w, such that A - w*B is singular. It is
37 * usually represented as the pair (alpha,beta), as there is reasonalbe
38 * interpretation for beta=0, and even for both being zero.
39 *
40 * A right generalized eigenvector corresponding to a generalized
41 * eigenvalue w for a pair of matrices (A,B) is a vector r such that
42 * (A - wB) * r = 0. A left generalized eigenvector is a vector l such
43 * that l**H * (A - wB) = 0, where l**H is the conjugate-transpose of l.
44 *
45 * When ZDRGEV is called, a number of matrix "sizes" ("n's") and a
46 * number of matrix "types" are specified. For each size ("n")
47 * and each type of matrix, a pair of matrices (A, B) will be generated
48 * and used for testing. For each matrix pair, the following tests
49 * will be performed and compared with the threshhold THRESH.
50 *
51 * Results from ZGGEV:
52 *
53 * (1) max over all left eigenvalue/-vector pairs (alpha/beta,l) of
54 *
55 * | VL**H * (beta A - alpha B) |/( ulp max(|beta A|, |alpha B|) )
56 *
57 * where VL**H is the conjugate-transpose of VL.
58 *
59 * (2) | |VL(i)| - 1 | / ulp and whether largest component real
60 *
61 * VL(i) denotes the i-th column of VL.
62 *
63 * (3) max over all left eigenvalue/-vector pairs (alpha/beta,r) of
64 *
65 * | (beta A - alpha B) * VR | / ( ulp max(|beta A|, |alpha B|) )
66 *
67 * (4) | |VR(i)| - 1 | / ulp and whether largest component real
68 *
69 * VR(i) denotes the i-th column of VR.
70 *
71 * (5) W(full) = W(partial)
72 * W(full) denotes the eigenvalues computed when both l and r
73 * are also computed, and W(partial) denotes the eigenvalues
74 * computed when only W, only W and r, or only W and l are
75 * computed.
76 *
77 * (6) VL(full) = VL(partial)
78 * VL(full) denotes the left eigenvectors computed when both l
79 * and r are computed, and VL(partial) denotes the result
80 * when only l is computed.
81 *
82 * (7) VR(full) = VR(partial)
83 * VR(full) denotes the right eigenvectors computed when both l
84 * and r are also computed, and VR(partial) denotes the result
85 * when only l is computed.
86 *
87 *
88 * Test Matrices
89 * ---- --------
90 *
91 * The sizes of the test matrices are specified by an array
92 * NN(1:NSIZES); the value of each element NN(j) specifies one size.
93 * The "types" are specified by a logical array DOTYPE( 1:NTYPES ); if
94 * DOTYPE(j) is .TRUE., then matrix type "j" will be generated.
95 * Currently, the list of possible types is:
96 *
97 * (1) ( 0, 0 ) (a pair of zero matrices)
98 *
99 * (2) ( I, 0 ) (an identity and a zero matrix)
100 *
101 * (3) ( 0, I ) (an identity and a zero matrix)
102 *
103 * (4) ( I, I ) (a pair of identity matrices)
104 *
105 * t t
106 * (5) ( J , J ) (a pair of transposed Jordan blocks)
107 *
108 * t ( I 0 )
109 * (6) ( X, Y ) where X = ( J 0 ) and Y = ( t )
110 * ( 0 I ) ( 0 J )
111 * and I is a k x k identity and J a (k+1)x(k+1)
112 * Jordan block; k=(N-1)/2
113 *
114 * (7) ( D, I ) where D is diag( 0, 1,..., N-1 ) (a diagonal
115 * matrix with those diagonal entries.)
116 * (8) ( I, D )
117 *
118 * (9) ( big*D, small*I ) where "big" is near overflow and small=1/big
119 *
120 * (10) ( small*D, big*I )
121 *
122 * (11) ( big*I, small*D )
123 *
124 * (12) ( small*I, big*D )
125 *
126 * (13) ( big*D, big*I )
127 *
128 * (14) ( small*D, small*I )
129 *
130 * (15) ( D1, D2 ) where D1 is diag( 0, 0, 1, ..., N-3, 0 ) and
131 * D2 is diag( 0, N-3, N-4,..., 1, 0, 0 )
132 * t t
133 * (16) Q ( J , J ) Z where Q and Z are random orthogonal matrices.
134 *
135 * (17) Q ( T1, T2 ) Z where T1 and T2 are upper triangular matrices
136 * with random O(1) entries above the diagonal
137 * and diagonal entries diag(T1) =
138 * ( 0, 0, 1, ..., N-3, 0 ) and diag(T2) =
139 * ( 0, N-3, N-4,..., 1, 0, 0 )
140 *
141 * (18) Q ( T1, T2 ) Z diag(T1) = ( 0, 0, 1, 1, s, ..., s, 0 )
142 * diag(T2) = ( 0, 1, 0, 1,..., 1, 0 )
143 * s = machine precision.
144 *
145 * (19) Q ( T1, T2 ) Z diag(T1)=( 0,0,1,1, 1-d, ..., 1-(N-5)*d=s, 0 )
146 * diag(T2) = ( 0, 1, 0, 1, ..., 1, 0 )
147 *
148 * N-5
149 * (20) Q ( T1, T2 ) Z diag(T1)=( 0, 0, 1, 1, a, ..., a =s, 0 )
150 * diag(T2) = ( 0, 1, 0, 1, ..., 1, 0, 0 )
151 *
152 * (21) Q ( T1, T2 ) Z diag(T1)=( 0, 0, 1, r1, r2, ..., r(N-4), 0 )
153 * diag(T2) = ( 0, 1, 0, 1, ..., 1, 0, 0 )
154 * where r1,..., r(N-4) are random.
155 *
156 * (22) Q ( big*T1, small*T2 ) Z diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
157 * diag(T2) = ( 0, 1, ..., 1, 0, 0 )
158 *
159 * (23) Q ( small*T1, big*T2 ) Z diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
160 * diag(T2) = ( 0, 1, ..., 1, 0, 0 )
161 *
162 * (24) Q ( small*T1, small*T2 ) Z diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
163 * diag(T2) = ( 0, 1, ..., 1, 0, 0 )
164 *
165 * (25) Q ( big*T1, big*T2 ) Z diag(T1) = ( 0, 0, 1, ..., N-3, 0 )
166 * diag(T2) = ( 0, 1, ..., 1, 0, 0 )
167 *
168 * (26) Q ( T1, T2 ) Z where T1 and T2 are random upper-triangular
169 * matrices.
170 *
171 *
172 * Arguments
173 * =========
174 *
175 * NSIZES (input) INTEGER
176 * The number of sizes of matrices to use. If it is zero,
177 * ZDRGES does nothing. NSIZES >= 0.
178 *
179 * NN (input) INTEGER array, dimension (NSIZES)
180 * An array containing the sizes to be used for the matrices.
181 * Zero values will be skipped. NN >= 0.
182 *
183 * NTYPES (input) INTEGER
184 * The number of elements in DOTYPE. If it is zero, ZDRGEV
185 * does nothing. It must be at least zero. If it is MAXTYP+1
186 * and NSIZES is 1, then an additional type, MAXTYP+1 is
187 * defined, which is to use whatever matrix is in A. This
188 * is only useful if DOTYPE(1:MAXTYP) is .FALSE. and
189 * DOTYPE(MAXTYP+1) is .TRUE. .
190 *
191 * DOTYPE (input) LOGICAL array, dimension (NTYPES)
192 * If DOTYPE(j) is .TRUE., then for each size in NN a
193 * matrix of that size and of type j will be generated.
194 * If NTYPES is smaller than the maximum number of types
195 * defined (PARAMETER MAXTYP), then types NTYPES+1 through
196 * MAXTYP will not be generated. If NTYPES is larger
197 * than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES)
198 * will be ignored.
199 *
200 * ISEED (input/output) INTEGER array, dimension (4)
201 * On entry ISEED specifies the seed of the random number
202 * generator. The array elements should be between 0 and 4095;
203 * if not they will be reduced mod 4096. Also, ISEED(4) must
204 * be odd. The random number generator uses a linear
205 * congruential sequence limited to small integers, and so
206 * should produce machine independent random numbers. The
207 * values of ISEED are changed on exit, and can be used in the
208 * next call to ZDRGES to continue the same random number
209 * sequence.
210 *
211 * THRESH (input) DOUBLE PRECISION
212 * A test will count as "failed" if the "error", computed as
213 * described above, exceeds THRESH. Note that the error is
214 * scaled to be O(1), so THRESH should be a reasonably small
215 * multiple of 1, e.g., 10 or 100. In particular, it should
216 * not depend on the precision (single vs. double) or the size
217 * of the matrix. It must be at least zero.
218 *
219 * NOUNIT (input) INTEGER
220 * The FORTRAN unit number for printing out error messages
221 * (e.g., if a routine returns IERR not equal to 0.)
222 *
223 * A (input/workspace) COMPLEX*16 array, dimension(LDA, max(NN))
224 * Used to hold the original A matrix. Used as input only
225 * if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and
226 * DOTYPE(MAXTYP+1)=.TRUE.
227 *
228 * LDA (input) INTEGER
229 * The leading dimension of A, B, S, and T.
230 * It must be at least 1 and at least max( NN ).
231 *
232 * B (input/workspace) COMPLEX*16 array, dimension(LDA, max(NN))
233 * Used to hold the original B matrix. Used as input only
234 * if NTYPES=MAXTYP+1, DOTYPE(1:MAXTYP)=.FALSE., and
235 * DOTYPE(MAXTYP+1)=.TRUE.
236 *
237 * S (workspace) COMPLEX*16 array, dimension (LDA, max(NN))
238 * The Schur form matrix computed from A by ZGGEV. On exit, S
239 * contains the Schur form matrix corresponding to the matrix
240 * in A.
241 *
242 * T (workspace) COMPLEX*16 array, dimension (LDA, max(NN))
243 * The upper triangular matrix computed from B by ZGGEV.
244 *
245 * Q (workspace) COMPLEX*16 array, dimension (LDQ, max(NN))
246 * The (left) eigenvectors matrix computed by ZGGEV.
247 *
248 * LDQ (input) INTEGER
249 * The leading dimension of Q and Z. It must
250 * be at least 1 and at least max( NN ).
251 *
252 * Z (workspace) COMPLEX*16 array, dimension( LDQ, max(NN) )
253 * The (right) orthogonal matrix computed by ZGGEV.
254 *
255 * QE (workspace) COMPLEX*16 array, dimension( LDQ, max(NN) )
256 * QE holds the computed right or left eigenvectors.
257 *
258 * LDQE (input) INTEGER
259 * The leading dimension of QE. LDQE >= max(1,max(NN)).
260 *
261 * ALPHA (workspace) COMPLEX*16 array, dimension (max(NN))
262 * BETA (workspace) COMPLEX*16 array, dimension (max(NN))
263 * The generalized eigenvalues of (A,B) computed by ZGGEV.
264 * ( ALPHAR(k)+ALPHAI(k)*i ) / BETA(k) is the k-th
265 * generalized eigenvalue of A and B.
266 *
267 * ALPHA1 (workspace) COMPLEX*16 array, dimension (max(NN))
268 * BETA1 (workspace) COMPLEX*16 array, dimension (max(NN))
269 * Like ALPHAR, ALPHAI, BETA, these arrays contain the
270 * eigenvalues of A and B, but those computed when ZGGEV only
271 * computes a partial eigendecomposition, i.e. not the
272 * eigenvalues and left and right eigenvectors.
273 *
274 * WORK (workspace) COMPLEX*16 array, dimension (LWORK)
275 *
276 * LWORK (input) INTEGER
277 * The number of entries in WORK. LWORK >= N*(N+1)
278 *
279 * RWORK (workspace) DOUBLE PRECISION array, dimension (8*N)
280 * Real workspace.
281 *
282 * RESULT (output) DOUBLE PRECISION array, dimension (2)
283 * The values computed by the tests described above.
284 * The values are currently limited to 1/ulp, to avoid overflow.
285 *
286 * INFO (output) INTEGER
287 * = 0: successful exit
288 * < 0: if INFO = -i, the i-th argument had an illegal value.
289 * > 0: A routine returned an error code. INFO is the
290 * absolute value of the INFO value returned.
291 *
292 * =====================================================================
293 *
294 * .. Parameters ..
295 DOUBLE PRECISION ZERO, ONE
296 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
297 COMPLEX*16 CZERO, CONE
298 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ),
299 $ CONE = ( 1.0D+0, 0.0D+0 ) )
300 INTEGER MAXTYP
301 PARAMETER ( MAXTYP = 26 )
302 * ..
303 * .. Local Scalars ..
304 LOGICAL BADNN
305 INTEGER I, IADD, IERR, IN, J, JC, JR, JSIZE, JTYPE,
306 $ MAXWRK, MINWRK, MTYPES, N, N1, NB, NERRS,
307 $ NMATS, NMAX, NTESTT
308 DOUBLE PRECISION SAFMAX, SAFMIN, ULP, ULPINV
309 COMPLEX*16 CTEMP
310 * ..
311 * .. Local Arrays ..
312 LOGICAL LASIGN( MAXTYP ), LBSIGN( MAXTYP )
313 INTEGER IOLDSD( 4 ), KADD( 6 ), KAMAGN( MAXTYP ),
314 $ KATYPE( MAXTYP ), KAZERO( MAXTYP ),
315 $ KBMAGN( MAXTYP ), KBTYPE( MAXTYP ),
316 $ KBZERO( MAXTYP ), KCLASS( MAXTYP ),
317 $ KTRIAN( MAXTYP ), KZ1( 6 ), KZ2( 6 )
318 DOUBLE PRECISION RMAGN( 0: 3 )
319 * ..
320 * .. External Functions ..
321 INTEGER ILAENV
322 DOUBLE PRECISION DLAMCH
323 COMPLEX*16 ZLARND
324 EXTERNAL ILAENV, DLAMCH, ZLARND
325 * ..
326 * .. External Subroutines ..
327 EXTERNAL ALASVM, DLABAD, XERBLA, ZGET52, ZGGEV, ZLACPY,
328 $ ZLARFG, ZLASET, ZLATM4, ZUNM2R
329 * ..
330 * .. Intrinsic Functions ..
331 INTRINSIC ABS, DBLE, DCONJG, MAX, MIN, SIGN
332 * ..
333 * .. Data statements ..
334 DATA KCLASS / 15*1, 10*2, 1*3 /
335 DATA KZ1 / 0, 1, 2, 1, 3, 3 /
336 DATA KZ2 / 0, 0, 1, 2, 1, 1 /
337 DATA KADD / 0, 0, 0, 0, 3, 2 /
338 DATA KATYPE / 0, 1, 0, 1, 2, 3, 4, 1, 4, 4, 1, 1, 4,
339 $ 4, 4, 2, 4, 5, 8, 7, 9, 4*4, 0 /
340 DATA KBTYPE / 0, 0, 1, 1, 2, -3, 1, 4, 1, 1, 4, 4,
341 $ 1, 1, -4, 2, -4, 8*8, 0 /
342 DATA KAZERO / 6*1, 2, 1, 2*2, 2*1, 2*2, 3, 1, 3,
343 $ 4*5, 4*3, 1 /
344 DATA KBZERO / 6*1, 1, 2, 2*1, 2*2, 2*1, 4, 1, 4,
345 $ 4*6, 4*4, 1 /
346 DATA KAMAGN / 8*1, 2, 3, 2, 3, 2, 3, 7*1, 2, 3, 3,
347 $ 2, 1 /
348 DATA KBMAGN / 8*1, 3, 2, 3, 2, 2, 3, 7*1, 3, 2, 3,
349 $ 2, 1 /
350 DATA KTRIAN / 16*0, 10*1 /
351 DATA LASIGN / 6*.FALSE., .TRUE., .FALSE., 2*.TRUE.,
352 $ 2*.FALSE., 3*.TRUE., .FALSE., .TRUE.,
353 $ 3*.FALSE., 5*.TRUE., .FALSE. /
354 DATA LBSIGN / 7*.FALSE., .TRUE., 2*.FALSE.,
355 $ 2*.TRUE., 2*.FALSE., .TRUE., .FALSE., .TRUE.,
356 $ 9*.FALSE. /
357 * ..
358 * .. Executable Statements ..
359 *
360 * Check for errors
361 *
362 INFO = 0
363 *
364 BADNN = .FALSE.
365 NMAX = 1
366 DO 10 J = 1, NSIZES
367 NMAX = MAX( NMAX, NN( J ) )
368 IF( NN( J ).LT.0 )
369 $ BADNN = .TRUE.
370 10 CONTINUE
371 *
372 IF( NSIZES.LT.0 ) THEN
373 INFO = -1
374 ELSE IF( BADNN ) THEN
375 INFO = -2
376 ELSE IF( NTYPES.LT.0 ) THEN
377 INFO = -3
378 ELSE IF( THRESH.LT.ZERO ) THEN
379 INFO = -6
380 ELSE IF( LDA.LE.1 .OR. LDA.LT.NMAX ) THEN
381 INFO = -9
382 ELSE IF( LDQ.LE.1 .OR. LDQ.LT.NMAX ) THEN
383 INFO = -14
384 ELSE IF( LDQE.LE.1 .OR. LDQE.LT.NMAX ) THEN
385 INFO = -17
386 END IF
387 *
388 * Compute workspace
389 * (Note: Comments in the code beginning "Workspace:" describe the
390 * minimal amount of workspace needed at that point in the code,
391 * as well as the preferred amount for good performance.
392 * NB refers to the optimal block size for the immediately
393 * following subroutine, as returned by ILAENV.
394 *
395 MINWRK = 1
396 IF( INFO.EQ.0 .AND. LWORK.GE.1 ) THEN
397 MINWRK = NMAX*( NMAX+1 )
398 NB = MAX( 1, ILAENV( 1, 'ZGEQRF', ' ', NMAX, NMAX, -1, -1 ),
399 $ ILAENV( 1, 'ZUNMQR', 'LC', NMAX, NMAX, NMAX, -1 ),
400 $ ILAENV( 1, 'ZUNGQR', ' ', NMAX, NMAX, NMAX, -1 ) )
401 MAXWRK = MAX( 2*NMAX, NMAX*( NB+1 ), NMAX*( NMAX+1 ) )
402 WORK( 1 ) = MAXWRK
403 END IF
404 *
405 IF( LWORK.LT.MINWRK )
406 $ INFO = -23
407 *
408 IF( INFO.NE.0 ) THEN
409 CALL XERBLA( 'ZDRGEV', -INFO )
410 RETURN
411 END IF
412 *
413 * Quick return if possible
414 *
415 IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 )
416 $ RETURN
417 *
418 ULP = DLAMCH( 'Precision' )
419 SAFMIN = DLAMCH( 'Safe minimum' )
420 SAFMIN = SAFMIN / ULP
421 SAFMAX = ONE / SAFMIN
422 CALL DLABAD( SAFMIN, SAFMAX )
423 ULPINV = ONE / ULP
424 *
425 * The values RMAGN(2:3) depend on N, see below.
426 *
427 RMAGN( 0 ) = ZERO
428 RMAGN( 1 ) = ONE
429 *
430 * Loop over sizes, types
431 *
432 NTESTT = 0
433 NERRS = 0
434 NMATS = 0
435 *
436 DO 220 JSIZE = 1, NSIZES
437 N = NN( JSIZE )
438 N1 = MAX( 1, N )
439 RMAGN( 2 ) = SAFMAX*ULP / DBLE( N1 )
440 RMAGN( 3 ) = SAFMIN*ULPINV*N1
441 *
442 IF( NSIZES.NE.1 ) THEN
443 MTYPES = MIN( MAXTYP, NTYPES )
444 ELSE
445 MTYPES = MIN( MAXTYP+1, NTYPES )
446 END IF
447 *
448 DO 210 JTYPE = 1, MTYPES
449 IF( .NOT.DOTYPE( JTYPE ) )
450 $ GO TO 210
451 NMATS = NMATS + 1
452 *
453 * Save ISEED in case of an error.
454 *
455 DO 20 J = 1, 4
456 IOLDSD( J ) = ISEED( J )
457 20 CONTINUE
458 *
459 * Generate test matrices A and B
460 *
461 * Description of control parameters:
462 *
463 * KZLASS: =1 means w/o rotation, =2 means w/ rotation,
464 * =3 means random.
465 * KATYPE: the "type" to be passed to ZLATM4 for computing A.
466 * KAZERO: the pattern of zeros on the diagonal for A:
467 * =1: ( xxx ), =2: (0, xxx ) =3: ( 0, 0, xxx, 0 ),
468 * =4: ( 0, xxx, 0, 0 ), =5: ( 0, 0, 1, xxx, 0 ),
469 * =6: ( 0, 1, 0, xxx, 0 ). (xxx means a string of
470 * non-zero entries.)
471 * KAMAGN: the magnitude of the matrix: =0: zero, =1: O(1),
472 * =2: large, =3: small.
473 * LASIGN: .TRUE. if the diagonal elements of A are to be
474 * multiplied by a random magnitude 1 number.
475 * KBTYPE, KBZERO, KBMAGN, LBSIGN: the same, but for B.
476 * KTRIAN: =0: don't fill in the upper triangle, =1: do.
477 * KZ1, KZ2, KADD: used to implement KAZERO and KBZERO.
478 * RMAGN: used to implement KAMAGN and KBMAGN.
479 *
480 IF( MTYPES.GT.MAXTYP )
481 $ GO TO 100
482 IERR = 0
483 IF( KCLASS( JTYPE ).LT.3 ) THEN
484 *
485 * Generate A (w/o rotation)
486 *
487 IF( ABS( KATYPE( JTYPE ) ).EQ.3 ) THEN
488 IN = 2*( ( N-1 ) / 2 ) + 1
489 IF( IN.NE.N )
490 $ CALL ZLASET( 'Full', N, N, CZERO, CZERO, A, LDA )
491 ELSE
492 IN = N
493 END IF
494 CALL ZLATM4( KATYPE( JTYPE ), IN, KZ1( KAZERO( JTYPE ) ),
495 $ KZ2( KAZERO( JTYPE ) ), LASIGN( JTYPE ),
496 $ RMAGN( KAMAGN( JTYPE ) ), ULP,
497 $ RMAGN( KTRIAN( JTYPE )*KAMAGN( JTYPE ) ), 2,
498 $ ISEED, A, LDA )
499 IADD = KADD( KAZERO( JTYPE ) )
500 IF( IADD.GT.0 .AND. IADD.LE.N )
501 $ A( IADD, IADD ) = RMAGN( KAMAGN( JTYPE ) )
502 *
503 * Generate B (w/o rotation)
504 *
505 IF( ABS( KBTYPE( JTYPE ) ).EQ.3 ) THEN
506 IN = 2*( ( N-1 ) / 2 ) + 1
507 IF( IN.NE.N )
508 $ CALL ZLASET( 'Full', N, N, CZERO, CZERO, B, LDA )
509 ELSE
510 IN = N
511 END IF
512 CALL ZLATM4( KBTYPE( JTYPE ), IN, KZ1( KBZERO( JTYPE ) ),
513 $ KZ2( KBZERO( JTYPE ) ), LBSIGN( JTYPE ),
514 $ RMAGN( KBMAGN( JTYPE ) ), ONE,
515 $ RMAGN( KTRIAN( JTYPE )*KBMAGN( JTYPE ) ), 2,
516 $ ISEED, B, LDA )
517 IADD = KADD( KBZERO( JTYPE ) )
518 IF( IADD.NE.0 .AND. IADD.LE.N )
519 $ B( IADD, IADD ) = RMAGN( KBMAGN( JTYPE ) )
520 *
521 IF( KCLASS( JTYPE ).EQ.2 .AND. N.GT.0 ) THEN
522 *
523 * Include rotations
524 *
525 * Generate Q, Z as Householder transformations times
526 * a diagonal matrix.
527 *
528 DO 40 JC = 1, N - 1
529 DO 30 JR = JC, N
530 Q( JR, JC ) = ZLARND( 3, ISEED )
531 Z( JR, JC ) = ZLARND( 3, ISEED )
532 30 CONTINUE
533 CALL ZLARFG( N+1-JC, Q( JC, JC ), Q( JC+1, JC ), 1,
534 $ WORK( JC ) )
535 WORK( 2*N+JC ) = SIGN( ONE, DBLE( Q( JC, JC ) ) )
536 Q( JC, JC ) = CONE
537 CALL ZLARFG( N+1-JC, Z( JC, JC ), Z( JC+1, JC ), 1,
538 $ WORK( N+JC ) )
539 WORK( 3*N+JC ) = SIGN( ONE, DBLE( Z( JC, JC ) ) )
540 Z( JC, JC ) = CONE
541 40 CONTINUE
542 CTEMP = ZLARND( 3, ISEED )
543 Q( N, N ) = CONE
544 WORK( N ) = CZERO
545 WORK( 3*N ) = CTEMP / ABS( CTEMP )
546 CTEMP = ZLARND( 3, ISEED )
547 Z( N, N ) = CONE
548 WORK( 2*N ) = CZERO
549 WORK( 4*N ) = CTEMP / ABS( CTEMP )
550 *
551 * Apply the diagonal matrices
552 *
553 DO 60 JC = 1, N
554 DO 50 JR = 1, N
555 A( JR, JC ) = WORK( 2*N+JR )*
556 $ DCONJG( WORK( 3*N+JC ) )*
557 $ A( JR, JC )
558 B( JR, JC ) = WORK( 2*N+JR )*
559 $ DCONJG( WORK( 3*N+JC ) )*
560 $ B( JR, JC )
561 50 CONTINUE
562 60 CONTINUE
563 CALL ZUNM2R( 'L', 'N', N, N, N-1, Q, LDQ, WORK, A,
564 $ LDA, WORK( 2*N+1 ), IERR )
565 IF( IERR.NE.0 )
566 $ GO TO 90
567 CALL ZUNM2R( 'R', 'C', N, N, N-1, Z, LDQ, WORK( N+1 ),
568 $ A, LDA, WORK( 2*N+1 ), IERR )
569 IF( IERR.NE.0 )
570 $ GO TO 90
571 CALL ZUNM2R( 'L', 'N', N, N, N-1, Q, LDQ, WORK, B,
572 $ LDA, WORK( 2*N+1 ), IERR )
573 IF( IERR.NE.0 )
574 $ GO TO 90
575 CALL ZUNM2R( 'R', 'C', N, N, N-1, Z, LDQ, WORK( N+1 ),
576 $ B, LDA, WORK( 2*N+1 ), IERR )
577 IF( IERR.NE.0 )
578 $ GO TO 90
579 END IF
580 ELSE
581 *
582 * Random matrices
583 *
584 DO 80 JC = 1, N
585 DO 70 JR = 1, N
586 A( JR, JC ) = RMAGN( KAMAGN( JTYPE ) )*
587 $ ZLARND( 4, ISEED )
588 B( JR, JC ) = RMAGN( KBMAGN( JTYPE ) )*
589 $ ZLARND( 4, ISEED )
590 70 CONTINUE
591 80 CONTINUE
592 END IF
593 *
594 90 CONTINUE
595 *
596 IF( IERR.NE.0 ) THEN
597 WRITE( NOUNIT, FMT = 9999 )'Generator', IERR, N, JTYPE,
598 $ IOLDSD
599 INFO = ABS( IERR )
600 RETURN
601 END IF
602 *
603 100 CONTINUE
604 *
605 DO 110 I = 1, 7
606 RESULT( I ) = -ONE
607 110 CONTINUE
608 *
609 * Call ZGGEV to compute eigenvalues and eigenvectors.
610 *
611 CALL ZLACPY( ' ', N, N, A, LDA, S, LDA )
612 CALL ZLACPY( ' ', N, N, B, LDA, T, LDA )
613 CALL ZGGEV( 'V', 'V', N, S, LDA, T, LDA, ALPHA, BETA, Q,
614 $ LDQ, Z, LDQ, WORK, LWORK, RWORK, IERR )
615 IF( IERR.NE.0 .AND. IERR.NE.N+1 ) THEN
616 RESULT( 1 ) = ULPINV
617 WRITE( NOUNIT, FMT = 9999 )'ZGGEV1', IERR, N, JTYPE,
618 $ IOLDSD
619 INFO = ABS( IERR )
620 GO TO 190
621 END IF
622 *
623 * Do the tests (1) and (2)
624 *
625 CALL ZGET52( .TRUE., N, A, LDA, B, LDA, Q, LDQ, ALPHA, BETA,
626 $ WORK, RWORK, RESULT( 1 ) )
627 IF( RESULT( 2 ).GT.THRESH ) THEN
628 WRITE( NOUNIT, FMT = 9998 )'Left', 'ZGGEV1',
629 $ RESULT( 2 ), N, JTYPE, IOLDSD
630 END IF
631 *
632 * Do the tests (3) and (4)
633 *
634 CALL ZGET52( .FALSE., N, A, LDA, B, LDA, Z, LDQ, ALPHA,
635 $ BETA, WORK, RWORK, RESULT( 3 ) )
636 IF( RESULT( 4 ).GT.THRESH ) THEN
637 WRITE( NOUNIT, FMT = 9998 )'Right', 'ZGGEV1',
638 $ RESULT( 4 ), N, JTYPE, IOLDSD
639 END IF
640 *
641 * Do test (5)
642 *
643 CALL ZLACPY( ' ', N, N, A, LDA, S, LDA )
644 CALL ZLACPY( ' ', N, N, B, LDA, T, LDA )
645 CALL ZGGEV( 'N', 'N', N, S, LDA, T, LDA, ALPHA1, BETA1, Q,
646 $ LDQ, Z, LDQ, WORK, LWORK, RWORK, IERR )
647 IF( IERR.NE.0 .AND. IERR.NE.N+1 ) THEN
648 RESULT( 1 ) = ULPINV
649 WRITE( NOUNIT, FMT = 9999 )'ZGGEV2', IERR, N, JTYPE,
650 $ IOLDSD
651 INFO = ABS( IERR )
652 GO TO 190
653 END IF
654 *
655 DO 120 J = 1, N
656 IF( ALPHA( J ).NE.ALPHA1( J ) .OR. BETA( J ).NE.
657 $ BETA1( J ) )RESULT( 5 ) = ULPINV
658 120 CONTINUE
659 *
660 * Do test (6): Compute eigenvalues and left eigenvectors,
661 * and test them
662 *
663 CALL ZLACPY( ' ', N, N, A, LDA, S, LDA )
664 CALL ZLACPY( ' ', N, N, B, LDA, T, LDA )
665 CALL ZGGEV( 'V', 'N', N, S, LDA, T, LDA, ALPHA1, BETA1, QE,
666 $ LDQE, Z, LDQ, WORK, LWORK, RWORK, IERR )
667 IF( IERR.NE.0 .AND. IERR.NE.N+1 ) THEN
668 RESULT( 1 ) = ULPINV
669 WRITE( NOUNIT, FMT = 9999 )'ZGGEV3', IERR, N, JTYPE,
670 $ IOLDSD
671 INFO = ABS( IERR )
672 GO TO 190
673 END IF
674 *
675 DO 130 J = 1, N
676 IF( ALPHA( J ).NE.ALPHA1( J ) .OR. BETA( J ).NE.
677 $ BETA1( J ) )RESULT( 6 ) = ULPINV
678 130 CONTINUE
679 *
680 DO 150 J = 1, N
681 DO 140 JC = 1, N
682 IF( Q( J, JC ).NE.QE( J, JC ) )
683 $ RESULT( 6 ) = ULPINV
684 140 CONTINUE
685 150 CONTINUE
686 *
687 * Do test (7): Compute eigenvalues and right eigenvectors,
688 * and test them
689 *
690 CALL ZLACPY( ' ', N, N, A, LDA, S, LDA )
691 CALL ZLACPY( ' ', N, N, B, LDA, T, LDA )
692 CALL ZGGEV( 'N', 'V', N, S, LDA, T, LDA, ALPHA1, BETA1, Q,
693 $ LDQ, QE, LDQE, WORK, LWORK, RWORK, IERR )
694 IF( IERR.NE.0 .AND. IERR.NE.N+1 ) THEN
695 RESULT( 1 ) = ULPINV
696 WRITE( NOUNIT, FMT = 9999 )'ZGGEV4', IERR, N, JTYPE,
697 $ IOLDSD
698 INFO = ABS( IERR )
699 GO TO 190
700 END IF
701 *
702 DO 160 J = 1, N
703 IF( ALPHA( J ).NE.ALPHA1( J ) .OR. BETA( J ).NE.
704 $ BETA1( J ) )RESULT( 7 ) = ULPINV
705 160 CONTINUE
706 *
707 DO 180 J = 1, N
708 DO 170 JC = 1, N
709 IF( Z( J, JC ).NE.QE( J, JC ) )
710 $ RESULT( 7 ) = ULPINV
711 170 CONTINUE
712 180 CONTINUE
713 *
714 * End of Loop -- Check for RESULT(j) > THRESH
715 *
716 190 CONTINUE
717 *
718 NTESTT = NTESTT + 7
719 *
720 * Print out tests which fail.
721 *
722 DO 200 JR = 1, 7
723 IF( RESULT( JR ).GE.THRESH ) THEN
724 *
725 * If this is the first test to fail,
726 * print a header to the data file.
727 *
728 IF( NERRS.EQ.0 ) THEN
729 WRITE( NOUNIT, FMT = 9997 )'ZGV'
730 *
731 * Matrix types
732 *
733 WRITE( NOUNIT, FMT = 9996 )
734 WRITE( NOUNIT, FMT = 9995 )
735 WRITE( NOUNIT, FMT = 9994 )'Orthogonal'
736 *
737 * Tests performed
738 *
739 WRITE( NOUNIT, FMT = 9993 )
740 *
741 END IF
742 NERRS = NERRS + 1
743 IF( RESULT( JR ).LT.10000.0D0 ) THEN
744 WRITE( NOUNIT, FMT = 9992 )N, JTYPE, IOLDSD, JR,
745 $ RESULT( JR )
746 ELSE
747 WRITE( NOUNIT, FMT = 9991 )N, JTYPE, IOLDSD, JR,
748 $ RESULT( JR )
749 END IF
750 END IF
751 200 CONTINUE
752 *
753 210 CONTINUE
754 220 CONTINUE
755 *
756 * Summary
757 *
758 CALL ALASVM( 'ZGV', NOUNIT, NERRS, NTESTT, 0 )
759 *
760 WORK( 1 ) = MAXWRK
761 *
762 RETURN
763 *
764 9999 FORMAT( ' ZDRGEV: ', A, ' returned INFO=', I6, '.', / 3X, 'N=',
765 $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' )
766 *
767 9998 FORMAT( ' ZDRGEV: ', A, ' Eigenvectors from ', A, ' incorrectly ',
768 $ 'normalized.', / ' Bits of error=', 0P, G10.3, ',', 3X,
769 $ 'N=', I4, ', JTYPE=', I3, ', ISEED=(', 3( I4, ',' ), I5,
770 $ ')' )
771 *
772 9997 FORMAT( / 1X, A3, ' -- Complex Generalized eigenvalue problem ',
773 $ 'driver' )
774 *
775 9996 FORMAT( ' Matrix types (see ZDRGEV for details): ' )
776 *
777 9995 FORMAT( ' Special Matrices:', 23X,
778 $ '(J''=transposed Jordan block)',
779 $ / ' 1=(0,0) 2=(I,0) 3=(0,I) 4=(I,I) 5=(J'',J'') ',
780 $ '6=(diag(J'',I), diag(I,J''))', / ' Diagonal Matrices: ( ',
781 $ 'D=diag(0,1,2,...) )', / ' 7=(D,I) 9=(large*D, small*I',
782 $ ') 11=(large*I, small*D) 13=(large*D, large*I)', /
783 $ ' 8=(I,D) 10=(small*D, large*I) 12=(small*I, large*D) ',
784 $ ' 14=(small*D, small*I)', / ' 15=(D, reversed D)' )
785 9994 FORMAT( ' Matrices Rotated by Random ', A, ' Matrices U, V:',
786 $ / ' 16=Transposed Jordan Blocks 19=geometric ',
787 $ 'alpha, beta=0,1', / ' 17=arithm. alpha&beta ',
788 $ ' 20=arithmetic alpha, beta=0,1', / ' 18=clustered ',
789 $ 'alpha, beta=0,1 21=random alpha, beta=0,1',
790 $ / ' Large & Small Matrices:', / ' 22=(large, small) ',
791 $ '23=(small,large) 24=(small,small) 25=(large,large)',
792 $ / ' 26=random O(1) matrices.' )
793 *
794 9993 FORMAT( / ' Tests performed: ',
795 $ / ' 1 = max | ( b A - a B )''*l | / const.,',
796 $ / ' 2 = | |VR(i)| - 1 | / ulp,',
797 $ / ' 3 = max | ( b A - a B )*r | / const.',
798 $ / ' 4 = | |VL(i)| - 1 | / ulp,',
799 $ / ' 5 = 0 if W same no matter if r or l computed,',
800 $ / ' 6 = 0 if l same no matter if l computed,',
801 $ / ' 7 = 0 if r same no matter if r computed,', / 1X )
802 9992 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=',
803 $ 4( I4, ',' ), ' result ', I2, ' is', 0P, F8.2 )
804 9991 FORMAT( ' Matrix order=', I5, ', type=', I2, ', seed=',
805 $ 4( I4, ',' ), ' result ', I2, ' is', 1P, D10.3 )
806 *
807 * End of ZDRGEV
808 *
809 END