1 SUBROUTINE ZDRVVX( NSIZES, NN, NTYPES, DOTYPE, ISEED, THRESH,
2 $ NIUNIT, NOUNIT, A, LDA, H, W, W1, VL, LDVL, VR,
3 $ LDVR, LRE, LDLRE, RCONDV, RCNDV1, RCDVIN,
4 $ RCONDE, RCNDE1, RCDEIN, SCALE, SCALE1, RESULT,
5 $ WORK, NWORK, RWORK, 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 INTEGER INFO, LDA, LDLRE, LDVL, LDVR, NIUNIT, NOUNIT,
13 $ NSIZES, NTYPES, NWORK
14 DOUBLE PRECISION THRESH
15 * ..
16 * .. Array Arguments ..
17 LOGICAL DOTYPE( * )
18 INTEGER ISEED( 4 ), NN( * )
19 DOUBLE PRECISION RCDEIN( * ), RCDVIN( * ), RCNDE1( * ),
20 $ RCNDV1( * ), RCONDE( * ), RCONDV( * ),
21 $ RESULT( 11 ), RWORK( * ), SCALE( * ),
22 $ SCALE1( * )
23 COMPLEX*16 A( LDA, * ), H( LDA, * ), LRE( LDLRE, * ),
24 $ VL( LDVL, * ), VR( LDVR, * ), W( * ), W1( * ),
25 $ WORK( * )
26 * ..
27 *
28 * Purpose
29 * =======
30 *
31 * ZDRVVX checks the nonsymmetric eigenvalue problem expert driver
32 * ZGEEVX.
33 *
34 * ZDRVVX uses both test matrices generated randomly depending on
35 * data supplied in the calling sequence, as well as on data
36 * read from an input file and including precomputed condition
37 * numbers to which it compares the ones it computes.
38 *
39 * When ZDRVVX is called, a number of matrix "sizes" ("n's") and a
40 * number of matrix "types" are specified in the calling sequence.
41 * For each size ("n") and each type of matrix, one matrix will be
42 * generated and used to test the nonsymmetric eigenroutines. For
43 * each matrix, 9 tests will be performed:
44 *
45 * (1) | A * VR - VR * W | / ( n |A| ulp )
46 *
47 * Here VR is the matrix of unit right eigenvectors.
48 * W is a diagonal matrix with diagonal entries W(j).
49 *
50 * (2) | A**H * VL - VL * W**H | / ( n |A| ulp )
51 *
52 * Here VL is the matrix of unit left eigenvectors, A**H is the
53 * conjugate transpose of A, and W is as above.
54 *
55 * (3) | |VR(i)| - 1 | / ulp and largest component real
56 *
57 * VR(i) denotes the i-th column of VR.
58 *
59 * (4) | |VL(i)| - 1 | / ulp and largest component real
60 *
61 * VL(i) denotes the i-th column of VL.
62 *
63 * (5) W(full) = W(partial)
64 *
65 * W(full) denotes the eigenvalues computed when VR, VL, RCONDV
66 * and RCONDE are also computed, and W(partial) denotes the
67 * eigenvalues computed when only some of VR, VL, RCONDV, and
68 * RCONDE are computed.
69 *
70 * (6) VR(full) = VR(partial)
71 *
72 * VR(full) denotes the right eigenvectors computed when VL, RCONDV
73 * and RCONDE are computed, and VR(partial) denotes the result
74 * when only some of VL and RCONDV are computed.
75 *
76 * (7) VL(full) = VL(partial)
77 *
78 * VL(full) denotes the left eigenvectors computed when VR, RCONDV
79 * and RCONDE are computed, and VL(partial) denotes the result
80 * when only some of VR and RCONDV are computed.
81 *
82 * (8) 0 if SCALE, ILO, IHI, ABNRM (full) =
83 * SCALE, ILO, IHI, ABNRM (partial)
84 * 1/ulp otherwise
85 *
86 * SCALE, ILO, IHI and ABNRM describe how the matrix is balanced.
87 * (full) is when VR, VL, RCONDE and RCONDV are also computed, and
88 * (partial) is when some are not computed.
89 *
90 * (9) RCONDV(full) = RCONDV(partial)
91 *
92 * RCONDV(full) denotes the reciprocal condition numbers of the
93 * right eigenvectors computed when VR, VL and RCONDE are also
94 * computed. RCONDV(partial) denotes the reciprocal condition
95 * numbers when only some of VR, VL and RCONDE are computed.
96 *
97 * The "sizes" are specified by an array NN(1:NSIZES); the value of
98 * each element NN(j) specifies one size.
99 * The "types" are specified by a logical array DOTYPE( 1:NTYPES );
100 * if DOTYPE(j) is .TRUE., then matrix type "j" will be generated.
101 * Currently, the list of possible types is:
102 *
103 * (1) The zero matrix.
104 * (2) The identity matrix.
105 * (3) A (transposed) Jordan block, with 1's on the diagonal.
106 *
107 * (4) A diagonal matrix with evenly spaced entries
108 * 1, ..., ULP and random complex angles.
109 * (ULP = (first number larger than 1) - 1 )
110 * (5) A diagonal matrix with geometrically spaced entries
111 * 1, ..., ULP and random complex angles.
112 * (6) A diagonal matrix with "clustered" entries 1, ULP, ..., ULP
113 * and random complex angles.
114 *
115 * (7) Same as (4), but multiplied by a constant near
116 * the overflow threshold
117 * (8) Same as (4), but multiplied by a constant near
118 * the underflow threshold
119 *
120 * (9) A matrix of the form U' T U, where U is unitary and
121 * T has evenly spaced entries 1, ..., ULP with random complex
122 * angles on the diagonal and random O(1) entries in the upper
123 * triangle.
124 *
125 * (10) A matrix of the form U' T U, where U is unitary and
126 * T has geometrically spaced entries 1, ..., ULP with random
127 * complex angles on the diagonal and random O(1) entries in
128 * the upper triangle.
129 *
130 * (11) A matrix of the form U' T U, where U is unitary and
131 * T has "clustered" entries 1, ULP,..., ULP with random
132 * complex angles on the diagonal and random O(1) entries in
133 * the upper triangle.
134 *
135 * (12) A matrix of the form U' T U, where U is unitary and
136 * T has complex eigenvalues randomly chosen from
137 * ULP < |z| < 1 and random O(1) entries in the upper
138 * triangle.
139 *
140 * (13) A matrix of the form X' T X, where X has condition
141 * SQRT( ULP ) and T has evenly spaced entries 1, ..., ULP
142 * with random complex angles on the diagonal and random O(1)
143 * entries in the upper triangle.
144 *
145 * (14) A matrix of the form X' T X, where X has condition
146 * SQRT( ULP ) and T has geometrically spaced entries
147 * 1, ..., ULP with random complex angles on the diagonal
148 * and random O(1) entries in the upper triangle.
149 *
150 * (15) A matrix of the form X' T X, where X has condition
151 * SQRT( ULP ) and T has "clustered" entries 1, ULP,..., ULP
152 * with random complex angles on the diagonal and random O(1)
153 * entries in the upper triangle.
154 *
155 * (16) A matrix of the form X' T X, where X has condition
156 * SQRT( ULP ) and T has complex eigenvalues randomly chosen
157 * from ULP < |z| < 1 and random O(1) entries in the upper
158 * triangle.
159 *
160 * (17) Same as (16), but multiplied by a constant
161 * near the overflow threshold
162 * (18) Same as (16), but multiplied by a constant
163 * near the underflow threshold
164 *
165 * (19) Nonsymmetric matrix with random entries chosen from |z| < 1
166 * If N is at least 4, all entries in first two rows and last
167 * row, and first column and last two columns are zero.
168 * (20) Same as (19), but multiplied by a constant
169 * near the overflow threshold
170 * (21) Same as (19), but multiplied by a constant
171 * near the underflow threshold
172 *
173 * In addition, an input file will be read from logical unit number
174 * NIUNIT. The file contains matrices along with precomputed
175 * eigenvalues and reciprocal condition numbers for the eigenvalues
176 * and right eigenvectors. For these matrices, in addition to tests
177 * (1) to (9) we will compute the following two tests:
178 *
179 * (10) |RCONDV - RCDVIN| / cond(RCONDV)
180 *
181 * RCONDV is the reciprocal right eigenvector condition number
182 * computed by ZGEEVX and RCDVIN (the precomputed true value)
183 * is supplied as input. cond(RCONDV) is the condition number of
184 * RCONDV, and takes errors in computing RCONDV into account, so
185 * that the resulting quantity should be O(ULP). cond(RCONDV) is
186 * essentially given by norm(A)/RCONDE.
187 *
188 * (11) |RCONDE - RCDEIN| / cond(RCONDE)
189 *
190 * RCONDE is the reciprocal eigenvalue condition number
191 * computed by ZGEEVX and RCDEIN (the precomputed true value)
192 * is supplied as input. cond(RCONDE) is the condition number
193 * of RCONDE, and takes errors in computing RCONDE into account,
194 * so that the resulting quantity should be O(ULP). cond(RCONDE)
195 * is essentially given by norm(A)/RCONDV.
196 *
197 * Arguments
198 * ==========
199 *
200 * NSIZES (input) INTEGER
201 * The number of sizes of matrices to use. NSIZES must be at
202 * least zero. If it is zero, no randomly generated matrices
203 * are tested, but any test matrices read from NIUNIT will be
204 * tested.
205 *
206 * NN (input) INTEGER array, dimension (NSIZES)
207 * An array containing the sizes to be used for the matrices.
208 * Zero values will be skipped. The values must be at least
209 * zero.
210 *
211 * NTYPES (input) INTEGER
212 * The number of elements in DOTYPE. NTYPES must be at least
213 * zero. If it is zero, no randomly generated test matrices
214 * are tested, but and test matrices read from NIUNIT will be
215 * tested. If it is MAXTYP+1 and NSIZES is 1, then an
216 * additional type, MAXTYP+1 is defined, which is to use
217 * whatever matrix is in A. This is only useful if
218 * DOTYPE(1:MAXTYP) is .FALSE. and DOTYPE(MAXTYP+1) is .TRUE. .
219 *
220 * DOTYPE (input) LOGICAL array, dimension (NTYPES)
221 * If DOTYPE(j) is .TRUE., then for each size in NN a
222 * matrix of that size and of type j will be generated.
223 * If NTYPES is smaller than the maximum number of types
224 * defined (PARAMETER MAXTYP), then types NTYPES+1 through
225 * MAXTYP will not be generated. If NTYPES is larger
226 * than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES)
227 * will be ignored.
228 *
229 * ISEED (input/output) INTEGER array, dimension (4)
230 * On entry ISEED specifies the seed of the random number
231 * generator. The array elements should be between 0 and 4095;
232 * if not they will be reduced mod 4096. Also, ISEED(4) must
233 * be odd. The random number generator uses a linear
234 * congruential sequence limited to small integers, and so
235 * should produce machine independent random numbers. The
236 * values of ISEED are changed on exit, and can be used in the
237 * next call to ZDRVVX to continue the same random number
238 * sequence.
239 *
240 * THRESH (input) DOUBLE PRECISION
241 * A test will count as "failed" if the "error", computed as
242 * described above, exceeds THRESH. Note that the error
243 * is scaled to be O(1), so THRESH should be a reasonably
244 * small multiple of 1, e.g., 10 or 100. In particular,
245 * it should not depend on the precision (single vs. double)
246 * or the size of the matrix. It must be at least zero.
247 *
248 * NIUNIT (input) INTEGER
249 * The FORTRAN unit number for reading in the data file of
250 * problems to solve.
251 *
252 * NOUNIT (input) INTEGER
253 * The FORTRAN unit number for printing out error messages
254 * (e.g., if a routine returns INFO not equal to 0.)
255 *
256 * A (workspace) COMPLEX*16 array, dimension (LDA, max(NN,12))
257 * Used to hold the matrix whose eigenvalues are to be
258 * computed. On exit, A contains the last matrix actually used.
259 *
260 * LDA (input) INTEGER
261 * The leading dimension of A, and H. LDA must be at
262 * least 1 and at least max( NN, 12 ). (12 is the
263 * dimension of the largest matrix on the precomputed
264 * input file.)
265 *
266 * H (workspace) COMPLEX*16 array, dimension (LDA, max(NN,12))
267 * Another copy of the test matrix A, modified by ZGEEVX.
268 *
269 * W (workspace) COMPLEX*16 array, dimension (max(NN,12))
270 * Contains the eigenvalues of A.
271 *
272 * W1 (workspace) COMPLEX*16 array, dimension (max(NN,12))
273 * Like W, this array contains the eigenvalues of A,
274 * but those computed when ZGEEVX only computes a partial
275 * eigendecomposition, i.e. not the eigenvalues and left
276 * and right eigenvectors.
277 *
278 * VL (workspace) COMPLEX*16 array, dimension (LDVL, max(NN,12))
279 * VL holds the computed left eigenvectors.
280 *
281 * LDVL (input) INTEGER
282 * Leading dimension of VL. Must be at least max(1,max(NN,12)).
283 *
284 * VR (workspace) COMPLEX*16 array, dimension (LDVR, max(NN,12))
285 * VR holds the computed right eigenvectors.
286 *
287 * LDVR (input) INTEGER
288 * Leading dimension of VR. Must be at least max(1,max(NN,12)).
289 *
290 * LRE (workspace) COMPLEX*16 array, dimension (LDLRE, max(NN,12))
291 * LRE holds the computed right or left eigenvectors.
292 *
293 * LDLRE (input) INTEGER
294 * Leading dimension of LRE. Must be at least max(1,max(NN,12))
295 *
296 * RESULT (output) DOUBLE PRECISION array, dimension (11)
297 * The values computed by the seven tests described above.
298 * The values are currently limited to 1/ulp, to avoid
299 * overflow.
300 *
301 * WORK (workspace) COMPLEX*16 array, dimension (NWORK)
302 *
303 * NWORK (input) INTEGER
304 * The number of entries in WORK. This must be at least
305 * max(6*12+2*12**2,6*NN(j)+2*NN(j)**2) =
306 * max( 360 ,6*NN(j)+2*NN(j)**2) for all j.
307 *
308 * RWORK (workspace) DOUBLE PRECISION array, dimension (2*max(NN,12))
309 *
310 * INFO (output) INTEGER
311 * If 0, then successful exit.
312 * If <0, then input paramter -INFO is incorrect.
313 * If >0, ZLATMR, CLATMS, CLATME or ZGET23 returned an error
314 * code, and INFO is its absolute value.
315 *
316 *-----------------------------------------------------------------------
317 *
318 * Some Local Variables and Parameters:
319 * ---- ----- --------- --- ----------
320 *
321 * ZERO, ONE Real 0 and 1.
322 * MAXTYP The number of types defined.
323 * NMAX Largest value in NN or 12.
324 * NERRS The number of tests which have exceeded THRESH
325 * COND, CONDS,
326 * IMODE Values to be passed to the matrix generators.
327 * ANORM Norm of A; passed to matrix generators.
328 *
329 * OVFL, UNFL Overflow and underflow thresholds.
330 * ULP, ULPINV Finest relative precision and its inverse.
331 * RTULP, RTULPI Square roots of the previous 4 values.
332 *
333 * The following four arrays decode JTYPE:
334 * KTYPE(j) The general type (1-10) for type "j".
335 * KMODE(j) The MODE value to be passed to the matrix
336 * generator for type "j".
337 * KMAGN(j) The order of magnitude ( O(1),
338 * O(overflow^(1/2) ), O(underflow^(1/2) )
339 * KCONDS(j) Selectw whether CONDS is to be 1 or
340 * 1/sqrt(ulp). (0 means irrelevant.)
341 *
342 * =====================================================================
343 *
344 * .. Parameters ..
345 COMPLEX*16 CZERO
346 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ) )
347 COMPLEX*16 CONE
348 PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) )
349 DOUBLE PRECISION ZERO, ONE
350 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
351 INTEGER MAXTYP
352 PARAMETER ( MAXTYP = 21 )
353 * ..
354 * .. Local Scalars ..
355 LOGICAL BADNN
356 CHARACTER BALANC
357 CHARACTER*3 PATH
358 INTEGER I, IBAL, IINFO, IMODE, ISRT, ITYPE, IWK, J,
359 $ JCOL, JSIZE, JTYPE, MTYPES, N, NERRS, NFAIL,
360 $ NMAX, NNWORK, NTEST, NTESTF, NTESTT
361 DOUBLE PRECISION ANORM, COND, CONDS, OVFL, RTULP, RTULPI, ULP,
362 $ ULPINV, UNFL, WI, WR
363 * ..
364 * .. Local Arrays ..
365 CHARACTER BAL( 4 )
366 INTEGER IDUMMA( 1 ), IOLDSD( 4 ), KCONDS( MAXTYP ),
367 $ KMAGN( MAXTYP ), KMODE( MAXTYP ),
368 $ KTYPE( MAXTYP )
369 * ..
370 * .. External Functions ..
371 DOUBLE PRECISION DLAMCH
372 EXTERNAL DLAMCH
373 * ..
374 * .. External Subroutines ..
375 EXTERNAL DLABAD, DLASUM, XERBLA, ZGET23, ZLASET, ZLATME,
376 $ ZLATMR, ZLATMS
377 * ..
378 * .. Intrinsic Functions ..
379 INTRINSIC ABS, DCMPLX, MAX, MIN, SQRT
380 * ..
381 * .. Data statements ..
382 DATA KTYPE / 1, 2, 3, 5*4, 4*6, 6*6, 3*9 /
383 DATA KMAGN / 3*1, 1, 1, 1, 2, 3, 4*1, 1, 1, 1, 1, 2,
384 $ 3, 1, 2, 3 /
385 DATA KMODE / 3*0, 4, 3, 1, 4, 4, 4, 3, 1, 5, 4, 3,
386 $ 1, 5, 5, 5, 4, 3, 1 /
387 DATA KCONDS / 3*0, 5*0, 4*1, 6*2, 3*0 /
388 DATA BAL / 'N', 'P', 'S', 'B' /
389 * ..
390 * .. Executable Statements ..
391 *
392 PATH( 1: 1 ) = 'Zomplex precision'
393 PATH( 2: 3 ) = 'VX'
394 *
395 * Check for errors
396 *
397 NTESTT = 0
398 NTESTF = 0
399 INFO = 0
400 *
401 * Important constants
402 *
403 BADNN = .FALSE.
404 *
405 * 7 is the largest dimension in the input file of precomputed
406 * problems
407 *
408 NMAX = 7
409 DO 10 J = 1, NSIZES
410 NMAX = MAX( NMAX, NN( J ) )
411 IF( NN( J ).LT.0 )
412 $ BADNN = .TRUE.
413 10 CONTINUE
414 *
415 * Check for errors
416 *
417 IF( NSIZES.LT.0 ) THEN
418 INFO = -1
419 ELSE IF( BADNN ) THEN
420 INFO = -2
421 ELSE IF( NTYPES.LT.0 ) THEN
422 INFO = -3
423 ELSE IF( THRESH.LT.ZERO ) THEN
424 INFO = -6
425 ELSE IF( LDA.LT.1 .OR. LDA.LT.NMAX ) THEN
426 INFO = -10
427 ELSE IF( LDVL.LT.1 .OR. LDVL.LT.NMAX ) THEN
428 INFO = -15
429 ELSE IF( LDVR.LT.1 .OR. LDVR.LT.NMAX ) THEN
430 INFO = -17
431 ELSE IF( LDLRE.LT.1 .OR. LDLRE.LT.NMAX ) THEN
432 INFO = -19
433 ELSE IF( 6*NMAX+2*NMAX**2.GT.NWORK ) THEN
434 INFO = -30
435 END IF
436 *
437 IF( INFO.NE.0 ) THEN
438 CALL XERBLA( 'ZDRVVX', -INFO )
439 RETURN
440 END IF
441 *
442 * If nothing to do check on NIUNIT
443 *
444 IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 )
445 $ GO TO 160
446 *
447 * More Important constants
448 *
449 UNFL = DLAMCH( 'Safe minimum' )
450 OVFL = ONE / UNFL
451 CALL DLABAD( UNFL, OVFL )
452 ULP = DLAMCH( 'Precision' )
453 ULPINV = ONE / ULP
454 RTULP = SQRT( ULP )
455 RTULPI = ONE / RTULP
456 *
457 * Loop over sizes, types
458 *
459 NERRS = 0
460 *
461 DO 150 JSIZE = 1, NSIZES
462 N = NN( JSIZE )
463 IF( NSIZES.NE.1 ) THEN
464 MTYPES = MIN( MAXTYP, NTYPES )
465 ELSE
466 MTYPES = MIN( MAXTYP+1, NTYPES )
467 END IF
468 *
469 DO 140 JTYPE = 1, MTYPES
470 IF( .NOT.DOTYPE( JTYPE ) )
471 $ GO TO 140
472 *
473 * Save ISEED in case of an error.
474 *
475 DO 20 J = 1, 4
476 IOLDSD( J ) = ISEED( J )
477 20 CONTINUE
478 *
479 * Compute "A"
480 *
481 * Control parameters:
482 *
483 * KMAGN KCONDS KMODE KTYPE
484 * =1 O(1) 1 clustered 1 zero
485 * =2 large large clustered 2 identity
486 * =3 small exponential Jordan
487 * =4 arithmetic diagonal, (w/ eigenvalues)
488 * =5 random log symmetric, w/ eigenvalues
489 * =6 random general, w/ eigenvalues
490 * =7 random diagonal
491 * =8 random symmetric
492 * =9 random general
493 * =10 random triangular
494 *
495 IF( MTYPES.GT.MAXTYP )
496 $ GO TO 90
497 *
498 ITYPE = KTYPE( JTYPE )
499 IMODE = KMODE( JTYPE )
500 *
501 * Compute norm
502 *
503 GO TO ( 30, 40, 50 )KMAGN( JTYPE )
504 *
505 30 CONTINUE
506 ANORM = ONE
507 GO TO 60
508 *
509 40 CONTINUE
510 ANORM = OVFL*ULP
511 GO TO 60
512 *
513 50 CONTINUE
514 ANORM = UNFL*ULPINV
515 GO TO 60
516 *
517 60 CONTINUE
518 *
519 CALL ZLASET( 'Full', LDA, N, CZERO, CZERO, A, LDA )
520 IINFO = 0
521 COND = ULPINV
522 *
523 * Special Matrices -- Identity & Jordan block
524 *
525 * Zero
526 *
527 IF( ITYPE.EQ.1 ) THEN
528 IINFO = 0
529 *
530 ELSE IF( ITYPE.EQ.2 ) THEN
531 *
532 * Identity
533 *
534 DO 70 JCOL = 1, N
535 A( JCOL, JCOL ) = ANORM
536 70 CONTINUE
537 *
538 ELSE IF( ITYPE.EQ.3 ) THEN
539 *
540 * Jordan Block
541 *
542 DO 80 JCOL = 1, N
543 A( JCOL, JCOL ) = ANORM
544 IF( JCOL.GT.1 )
545 $ A( JCOL, JCOL-1 ) = ONE
546 80 CONTINUE
547 *
548 ELSE IF( ITYPE.EQ.4 ) THEN
549 *
550 * Diagonal Matrix, [Eigen]values Specified
551 *
552 CALL ZLATMS( N, N, 'S', ISEED, 'H', RWORK, IMODE, COND,
553 $ ANORM, 0, 0, 'N', A, LDA, WORK( N+1 ),
554 $ IINFO )
555 *
556 ELSE IF( ITYPE.EQ.5 ) THEN
557 *
558 * Symmetric, eigenvalues specified
559 *
560 CALL ZLATMS( N, N, 'S', ISEED, 'H', RWORK, IMODE, COND,
561 $ ANORM, N, N, 'N', A, LDA, WORK( N+1 ),
562 $ IINFO )
563 *
564 ELSE IF( ITYPE.EQ.6 ) THEN
565 *
566 * General, eigenvalues specified
567 *
568 IF( KCONDS( JTYPE ).EQ.1 ) THEN
569 CONDS = ONE
570 ELSE IF( KCONDS( JTYPE ).EQ.2 ) THEN
571 CONDS = RTULPI
572 ELSE
573 CONDS = ZERO
574 END IF
575 *
576 CALL ZLATME( N, 'D', ISEED, WORK, IMODE, COND, CONE, ' ',
577 $ 'T', 'T', 'T', RWORK, 4, CONDS, N, N, ANORM,
578 $ A, LDA, WORK( 2*N+1 ), IINFO )
579 *
580 ELSE IF( ITYPE.EQ.7 ) THEN
581 *
582 * Diagonal, random eigenvalues
583 *
584 CALL ZLATMR( N, N, 'D', ISEED, 'S', WORK, 6, ONE, CONE,
585 $ 'T', 'N', WORK( N+1 ), 1, ONE,
586 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, 0, 0,
587 $ ZERO, ANORM, 'NO', A, LDA, IDUMMA, IINFO )
588 *
589 ELSE IF( ITYPE.EQ.8 ) THEN
590 *
591 * Symmetric, random eigenvalues
592 *
593 CALL ZLATMR( N, N, 'D', ISEED, 'H', WORK, 6, ONE, CONE,
594 $ 'T', 'N', WORK( N+1 ), 1, ONE,
595 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N,
596 $ ZERO, ANORM, 'NO', A, LDA, IDUMMA, IINFO )
597 *
598 ELSE IF( ITYPE.EQ.9 ) THEN
599 *
600 * General, random eigenvalues
601 *
602 CALL ZLATMR( N, N, 'D', ISEED, 'N', WORK, 6, ONE, CONE,
603 $ 'T', 'N', WORK( N+1 ), 1, ONE,
604 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N,
605 $ ZERO, ANORM, 'NO', A, LDA, IDUMMA, IINFO )
606 IF( N.GE.4 ) THEN
607 CALL ZLASET( 'Full', 2, N, CZERO, CZERO, A, LDA )
608 CALL ZLASET( 'Full', N-3, 1, CZERO, CZERO, A( 3, 1 ),
609 $ LDA )
610 CALL ZLASET( 'Full', N-3, 2, CZERO, CZERO,
611 $ A( 3, N-1 ), LDA )
612 CALL ZLASET( 'Full', 1, N, CZERO, CZERO, A( N, 1 ),
613 $ LDA )
614 END IF
615 *
616 ELSE IF( ITYPE.EQ.10 ) THEN
617 *
618 * Triangular, random eigenvalues
619 *
620 CALL ZLATMR( N, N, 'D', ISEED, 'N', WORK, 6, ONE, CONE,
621 $ 'T', 'N', WORK( N+1 ), 1, ONE,
622 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, 0,
623 $ ZERO, ANORM, 'NO', A, LDA, IDUMMA, IINFO )
624 *
625 ELSE
626 *
627 IINFO = 1
628 END IF
629 *
630 IF( IINFO.NE.0 ) THEN
631 WRITE( NOUNIT, FMT = 9992 )'Generator', IINFO, N, JTYPE,
632 $ IOLDSD
633 INFO = ABS( IINFO )
634 RETURN
635 END IF
636 *
637 90 CONTINUE
638 *
639 * Test for minimal and generous workspace
640 *
641 DO 130 IWK = 1, 3
642 IF( IWK.EQ.1 ) THEN
643 NNWORK = 2*N
644 ELSE IF( IWK.EQ.2 ) THEN
645 NNWORK = 2*N + N**2
646 ELSE
647 NNWORK = 6*N + 2*N**2
648 END IF
649 NNWORK = MAX( NNWORK, 1 )
650 *
651 * Test for all balancing options
652 *
653 DO 120 IBAL = 1, 4
654 BALANC = BAL( IBAL )
655 *
656 * Perform tests
657 *
658 CALL ZGET23( .FALSE., 0, BALANC, JTYPE, THRESH,
659 $ IOLDSD, NOUNIT, N, A, LDA, H, W, W1, VL,
660 $ LDVL, VR, LDVR, LRE, LDLRE, RCONDV,
661 $ RCNDV1, RCDVIN, RCONDE, RCNDE1, RCDEIN,
662 $ SCALE, SCALE1, RESULT, WORK, NNWORK,
663 $ RWORK, INFO )
664 *
665 * Check for RESULT(j) > THRESH
666 *
667 NTEST = 0
668 NFAIL = 0
669 DO 100 J = 1, 9
670 IF( RESULT( J ).GE.ZERO )
671 $ NTEST = NTEST + 1
672 IF( RESULT( J ).GE.THRESH )
673 $ NFAIL = NFAIL + 1
674 100 CONTINUE
675 *
676 IF( NFAIL.GT.0 )
677 $ NTESTF = NTESTF + 1
678 IF( NTESTF.EQ.1 ) THEN
679 WRITE( NOUNIT, FMT = 9999 )PATH
680 WRITE( NOUNIT, FMT = 9998 )
681 WRITE( NOUNIT, FMT = 9997 )
682 WRITE( NOUNIT, FMT = 9996 )
683 WRITE( NOUNIT, FMT = 9995 )THRESH
684 NTESTF = 2
685 END IF
686 *
687 DO 110 J = 1, 9
688 IF( RESULT( J ).GE.THRESH ) THEN
689 WRITE( NOUNIT, FMT = 9994 )BALANC, N, IWK,
690 $ IOLDSD, JTYPE, J, RESULT( J )
691 END IF
692 110 CONTINUE
693 *
694 NERRS = NERRS + NFAIL
695 NTESTT = NTESTT + NTEST
696 *
697 120 CONTINUE
698 130 CONTINUE
699 140 CONTINUE
700 150 CONTINUE
701 *
702 160 CONTINUE
703 *
704 * Read in data from file to check accuracy of condition estimation.
705 * Assume input eigenvalues are sorted lexicographically (increasing
706 * by real part, then decreasing by imaginary part)
707 *
708 JTYPE = 0
709 170 CONTINUE
710 READ( NIUNIT, FMT = *, END = 220 )N, ISRT
711 *
712 * Read input data until N=0
713 *
714 IF( N.EQ.0 )
715 $ GO TO 220
716 JTYPE = JTYPE + 1
717 ISEED( 1 ) = JTYPE
718 DO 180 I = 1, N
719 READ( NIUNIT, FMT = * )( A( I, J ), J = 1, N )
720 180 CONTINUE
721 DO 190 I = 1, N
722 READ( NIUNIT, FMT = * )WR, WI, RCDEIN( I ), RCDVIN( I )
723 W1( I ) = DCMPLX( WR, WI )
724 190 CONTINUE
725 CALL ZGET23( .TRUE., ISRT, 'N', 22, THRESH, ISEED, NOUNIT, N, A,
726 $ LDA, H, W, W1, VL, LDVL, VR, LDVR, LRE, LDLRE,
727 $ RCONDV, RCNDV1, RCDVIN, RCONDE, RCNDE1, RCDEIN,
728 $ SCALE, SCALE1, RESULT, WORK, 6*N+2*N**2, RWORK,
729 $ INFO )
730 *
731 * Check for RESULT(j) > THRESH
732 *
733 NTEST = 0
734 NFAIL = 0
735 DO 200 J = 1, 11
736 IF( RESULT( J ).GE.ZERO )
737 $ NTEST = NTEST + 1
738 IF( RESULT( J ).GE.THRESH )
739 $ NFAIL = NFAIL + 1
740 200 CONTINUE
741 *
742 IF( NFAIL.GT.0 )
743 $ NTESTF = NTESTF + 1
744 IF( NTESTF.EQ.1 ) THEN
745 WRITE( NOUNIT, FMT = 9999 )PATH
746 WRITE( NOUNIT, FMT = 9998 )
747 WRITE( NOUNIT, FMT = 9997 )
748 WRITE( NOUNIT, FMT = 9996 )
749 WRITE( NOUNIT, FMT = 9995 )THRESH
750 NTESTF = 2
751 END IF
752 *
753 DO 210 J = 1, 11
754 IF( RESULT( J ).GE.THRESH ) THEN
755 WRITE( NOUNIT, FMT = 9993 )N, JTYPE, J, RESULT( J )
756 END IF
757 210 CONTINUE
758 *
759 NERRS = NERRS + NFAIL
760 NTESTT = NTESTT + NTEST
761 GO TO 170
762 220 CONTINUE
763 *
764 * Summary
765 *
766 CALL DLASUM( PATH, NOUNIT, NERRS, NTESTT )
767 *
768 9999 FORMAT( / 1X, A3, ' -- Complex Eigenvalue-Eigenvector ',
769 $ 'Decomposition Expert Driver',
770 $ / ' Matrix types (see ZDRVVX for details): ' )
771 *
772 9998 FORMAT( / ' Special Matrices:', / ' 1=Zero matrix. ',
773 $ ' ', ' 5=Diagonal: geometr. spaced entries.',
774 $ / ' 2=Identity matrix. ', ' 6=Diagona',
775 $ 'l: clustered entries.', / ' 3=Transposed Jordan block. ',
776 $ ' ', ' 7=Diagonal: large, evenly spaced.', / ' ',
777 $ '4=Diagonal: evenly spaced entries. ', ' 8=Diagonal: s',
778 $ 'mall, evenly spaced.' )
779 9997 FORMAT( ' Dense, Non-Symmetric Matrices:', / ' 9=Well-cond., ev',
780 $ 'enly spaced eigenvals.', ' 14=Ill-cond., geomet. spaced e',
781 $ 'igenals.', / ' 10=Well-cond., geom. spaced eigenvals. ',
782 $ ' 15=Ill-conditioned, clustered e.vals.', / ' 11=Well-cond',
783 $ 'itioned, clustered e.vals. ', ' 16=Ill-cond., random comp',
784 $ 'lex ', / ' 12=Well-cond., random complex ', ' ',
785 $ ' 17=Ill-cond., large rand. complx ', / ' 13=Ill-condi',
786 $ 'tioned, evenly spaced. ', ' 18=Ill-cond., small rand.',
787 $ ' complx ' )
788 9996 FORMAT( ' 19=Matrix with random O(1) entries. ', ' 21=Matrix ',
789 $ 'with small random entries.', / ' 20=Matrix with large ran',
790 $ 'dom entries. ', ' 22=Matrix read from input file', / )
791 9995 FORMAT( ' Tests performed with test threshold =', F8.2,
792 $ / / ' 1 = | A VR - VR W | / ( n |A| ulp ) ',
793 $ / ' 2 = | transpose(A) VL - VL W | / ( n |A| ulp ) ',
794 $ / ' 3 = | |VR(i)| - 1 | / ulp ',
795 $ / ' 4 = | |VL(i)| - 1 | / ulp ',
796 $ / ' 5 = 0 if W same no matter if VR or VL computed,',
797 $ ' 1/ulp otherwise', /
798 $ ' 6 = 0 if VR same no matter what else computed,',
799 $ ' 1/ulp otherwise', /
800 $ ' 7 = 0 if VL same no matter what else computed,',
801 $ ' 1/ulp otherwise', /
802 $ ' 8 = 0 if RCONDV same no matter what else computed,',
803 $ ' 1/ulp otherwise', /
804 $ ' 9 = 0 if SCALE, ILO, IHI, ABNRM same no matter what else',
805 $ ' computed, 1/ulp otherwise',
806 $ / ' 10 = | RCONDV - RCONDV(precomputed) | / cond(RCONDV),',
807 $ / ' 11 = | RCONDE - RCONDE(precomputed) | / cond(RCONDE),' )
808 9994 FORMAT( ' BALANC=''', A1, ''',N=', I4, ',IWK=', I1, ', seed=',
809 $ 4( I4, ',' ), ' type ', I2, ', test(', I2, ')=', G10.3 )
810 9993 FORMAT( ' N=', I5, ', input example =', I3, ', test(', I2, ')=',
811 $ G10.3 )
812 9992 FORMAT( ' ZDRVVX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
813 $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' )
814 *
815 RETURN
816 *
817 * End of ZDRVVX
818 *
819 END
2 $ NIUNIT, NOUNIT, A, LDA, H, W, W1, VL, LDVL, VR,
3 $ LDVR, LRE, LDLRE, RCONDV, RCNDV1, RCDVIN,
4 $ RCONDE, RCNDE1, RCDEIN, SCALE, SCALE1, RESULT,
5 $ WORK, NWORK, RWORK, 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 INTEGER INFO, LDA, LDLRE, LDVL, LDVR, NIUNIT, NOUNIT,
13 $ NSIZES, NTYPES, NWORK
14 DOUBLE PRECISION THRESH
15 * ..
16 * .. Array Arguments ..
17 LOGICAL DOTYPE( * )
18 INTEGER ISEED( 4 ), NN( * )
19 DOUBLE PRECISION RCDEIN( * ), RCDVIN( * ), RCNDE1( * ),
20 $ RCNDV1( * ), RCONDE( * ), RCONDV( * ),
21 $ RESULT( 11 ), RWORK( * ), SCALE( * ),
22 $ SCALE1( * )
23 COMPLEX*16 A( LDA, * ), H( LDA, * ), LRE( LDLRE, * ),
24 $ VL( LDVL, * ), VR( LDVR, * ), W( * ), W1( * ),
25 $ WORK( * )
26 * ..
27 *
28 * Purpose
29 * =======
30 *
31 * ZDRVVX checks the nonsymmetric eigenvalue problem expert driver
32 * ZGEEVX.
33 *
34 * ZDRVVX uses both test matrices generated randomly depending on
35 * data supplied in the calling sequence, as well as on data
36 * read from an input file and including precomputed condition
37 * numbers to which it compares the ones it computes.
38 *
39 * When ZDRVVX is called, a number of matrix "sizes" ("n's") and a
40 * number of matrix "types" are specified in the calling sequence.
41 * For each size ("n") and each type of matrix, one matrix will be
42 * generated and used to test the nonsymmetric eigenroutines. For
43 * each matrix, 9 tests will be performed:
44 *
45 * (1) | A * VR - VR * W | / ( n |A| ulp )
46 *
47 * Here VR is the matrix of unit right eigenvectors.
48 * W is a diagonal matrix with diagonal entries W(j).
49 *
50 * (2) | A**H * VL - VL * W**H | / ( n |A| ulp )
51 *
52 * Here VL is the matrix of unit left eigenvectors, A**H is the
53 * conjugate transpose of A, and W is as above.
54 *
55 * (3) | |VR(i)| - 1 | / ulp and largest component real
56 *
57 * VR(i) denotes the i-th column of VR.
58 *
59 * (4) | |VL(i)| - 1 | / ulp and largest component real
60 *
61 * VL(i) denotes the i-th column of VL.
62 *
63 * (5) W(full) = W(partial)
64 *
65 * W(full) denotes the eigenvalues computed when VR, VL, RCONDV
66 * and RCONDE are also computed, and W(partial) denotes the
67 * eigenvalues computed when only some of VR, VL, RCONDV, and
68 * RCONDE are computed.
69 *
70 * (6) VR(full) = VR(partial)
71 *
72 * VR(full) denotes the right eigenvectors computed when VL, RCONDV
73 * and RCONDE are computed, and VR(partial) denotes the result
74 * when only some of VL and RCONDV are computed.
75 *
76 * (7) VL(full) = VL(partial)
77 *
78 * VL(full) denotes the left eigenvectors computed when VR, RCONDV
79 * and RCONDE are computed, and VL(partial) denotes the result
80 * when only some of VR and RCONDV are computed.
81 *
82 * (8) 0 if SCALE, ILO, IHI, ABNRM (full) =
83 * SCALE, ILO, IHI, ABNRM (partial)
84 * 1/ulp otherwise
85 *
86 * SCALE, ILO, IHI and ABNRM describe how the matrix is balanced.
87 * (full) is when VR, VL, RCONDE and RCONDV are also computed, and
88 * (partial) is when some are not computed.
89 *
90 * (9) RCONDV(full) = RCONDV(partial)
91 *
92 * RCONDV(full) denotes the reciprocal condition numbers of the
93 * right eigenvectors computed when VR, VL and RCONDE are also
94 * computed. RCONDV(partial) denotes the reciprocal condition
95 * numbers when only some of VR, VL and RCONDE are computed.
96 *
97 * The "sizes" are specified by an array NN(1:NSIZES); the value of
98 * each element NN(j) specifies one size.
99 * The "types" are specified by a logical array DOTYPE( 1:NTYPES );
100 * if DOTYPE(j) is .TRUE., then matrix type "j" will be generated.
101 * Currently, the list of possible types is:
102 *
103 * (1) The zero matrix.
104 * (2) The identity matrix.
105 * (3) A (transposed) Jordan block, with 1's on the diagonal.
106 *
107 * (4) A diagonal matrix with evenly spaced entries
108 * 1, ..., ULP and random complex angles.
109 * (ULP = (first number larger than 1) - 1 )
110 * (5) A diagonal matrix with geometrically spaced entries
111 * 1, ..., ULP and random complex angles.
112 * (6) A diagonal matrix with "clustered" entries 1, ULP, ..., ULP
113 * and random complex angles.
114 *
115 * (7) Same as (4), but multiplied by a constant near
116 * the overflow threshold
117 * (8) Same as (4), but multiplied by a constant near
118 * the underflow threshold
119 *
120 * (9) A matrix of the form U' T U, where U is unitary and
121 * T has evenly spaced entries 1, ..., ULP with random complex
122 * angles on the diagonal and random O(1) entries in the upper
123 * triangle.
124 *
125 * (10) A matrix of the form U' T U, where U is unitary and
126 * T has geometrically spaced entries 1, ..., ULP with random
127 * complex angles on the diagonal and random O(1) entries in
128 * the upper triangle.
129 *
130 * (11) A matrix of the form U' T U, where U is unitary and
131 * T has "clustered" entries 1, ULP,..., ULP with random
132 * complex angles on the diagonal and random O(1) entries in
133 * the upper triangle.
134 *
135 * (12) A matrix of the form U' T U, where U is unitary and
136 * T has complex eigenvalues randomly chosen from
137 * ULP < |z| < 1 and random O(1) entries in the upper
138 * triangle.
139 *
140 * (13) A matrix of the form X' T X, where X has condition
141 * SQRT( ULP ) and T has evenly spaced entries 1, ..., ULP
142 * with random complex angles on the diagonal and random O(1)
143 * entries in the upper triangle.
144 *
145 * (14) A matrix of the form X' T X, where X has condition
146 * SQRT( ULP ) and T has geometrically spaced entries
147 * 1, ..., ULP with random complex angles on the diagonal
148 * and random O(1) entries in the upper triangle.
149 *
150 * (15) A matrix of the form X' T X, where X has condition
151 * SQRT( ULP ) and T has "clustered" entries 1, ULP,..., ULP
152 * with random complex angles on the diagonal and random O(1)
153 * entries in the upper triangle.
154 *
155 * (16) A matrix of the form X' T X, where X has condition
156 * SQRT( ULP ) and T has complex eigenvalues randomly chosen
157 * from ULP < |z| < 1 and random O(1) entries in the upper
158 * triangle.
159 *
160 * (17) Same as (16), but multiplied by a constant
161 * near the overflow threshold
162 * (18) Same as (16), but multiplied by a constant
163 * near the underflow threshold
164 *
165 * (19) Nonsymmetric matrix with random entries chosen from |z| < 1
166 * If N is at least 4, all entries in first two rows and last
167 * row, and first column and last two columns are zero.
168 * (20) Same as (19), but multiplied by a constant
169 * near the overflow threshold
170 * (21) Same as (19), but multiplied by a constant
171 * near the underflow threshold
172 *
173 * In addition, an input file will be read from logical unit number
174 * NIUNIT. The file contains matrices along with precomputed
175 * eigenvalues and reciprocal condition numbers for the eigenvalues
176 * and right eigenvectors. For these matrices, in addition to tests
177 * (1) to (9) we will compute the following two tests:
178 *
179 * (10) |RCONDV - RCDVIN| / cond(RCONDV)
180 *
181 * RCONDV is the reciprocal right eigenvector condition number
182 * computed by ZGEEVX and RCDVIN (the precomputed true value)
183 * is supplied as input. cond(RCONDV) is the condition number of
184 * RCONDV, and takes errors in computing RCONDV into account, so
185 * that the resulting quantity should be O(ULP). cond(RCONDV) is
186 * essentially given by norm(A)/RCONDE.
187 *
188 * (11) |RCONDE - RCDEIN| / cond(RCONDE)
189 *
190 * RCONDE is the reciprocal eigenvalue condition number
191 * computed by ZGEEVX and RCDEIN (the precomputed true value)
192 * is supplied as input. cond(RCONDE) is the condition number
193 * of RCONDE, and takes errors in computing RCONDE into account,
194 * so that the resulting quantity should be O(ULP). cond(RCONDE)
195 * is essentially given by norm(A)/RCONDV.
196 *
197 * Arguments
198 * ==========
199 *
200 * NSIZES (input) INTEGER
201 * The number of sizes of matrices to use. NSIZES must be at
202 * least zero. If it is zero, no randomly generated matrices
203 * are tested, but any test matrices read from NIUNIT will be
204 * tested.
205 *
206 * NN (input) INTEGER array, dimension (NSIZES)
207 * An array containing the sizes to be used for the matrices.
208 * Zero values will be skipped. The values must be at least
209 * zero.
210 *
211 * NTYPES (input) INTEGER
212 * The number of elements in DOTYPE. NTYPES must be at least
213 * zero. If it is zero, no randomly generated test matrices
214 * are tested, but and test matrices read from NIUNIT will be
215 * tested. If it is MAXTYP+1 and NSIZES is 1, then an
216 * additional type, MAXTYP+1 is defined, which is to use
217 * whatever matrix is in A. This is only useful if
218 * DOTYPE(1:MAXTYP) is .FALSE. and DOTYPE(MAXTYP+1) is .TRUE. .
219 *
220 * DOTYPE (input) LOGICAL array, dimension (NTYPES)
221 * If DOTYPE(j) is .TRUE., then for each size in NN a
222 * matrix of that size and of type j will be generated.
223 * If NTYPES is smaller than the maximum number of types
224 * defined (PARAMETER MAXTYP), then types NTYPES+1 through
225 * MAXTYP will not be generated. If NTYPES is larger
226 * than MAXTYP, DOTYPE(MAXTYP+1) through DOTYPE(NTYPES)
227 * will be ignored.
228 *
229 * ISEED (input/output) INTEGER array, dimension (4)
230 * On entry ISEED specifies the seed of the random number
231 * generator. The array elements should be between 0 and 4095;
232 * if not they will be reduced mod 4096. Also, ISEED(4) must
233 * be odd. The random number generator uses a linear
234 * congruential sequence limited to small integers, and so
235 * should produce machine independent random numbers. The
236 * values of ISEED are changed on exit, and can be used in the
237 * next call to ZDRVVX to continue the same random number
238 * sequence.
239 *
240 * THRESH (input) DOUBLE PRECISION
241 * A test will count as "failed" if the "error", computed as
242 * described above, exceeds THRESH. Note that the error
243 * is scaled to be O(1), so THRESH should be a reasonably
244 * small multiple of 1, e.g., 10 or 100. In particular,
245 * it should not depend on the precision (single vs. double)
246 * or the size of the matrix. It must be at least zero.
247 *
248 * NIUNIT (input) INTEGER
249 * The FORTRAN unit number for reading in the data file of
250 * problems to solve.
251 *
252 * NOUNIT (input) INTEGER
253 * The FORTRAN unit number for printing out error messages
254 * (e.g., if a routine returns INFO not equal to 0.)
255 *
256 * A (workspace) COMPLEX*16 array, dimension (LDA, max(NN,12))
257 * Used to hold the matrix whose eigenvalues are to be
258 * computed. On exit, A contains the last matrix actually used.
259 *
260 * LDA (input) INTEGER
261 * The leading dimension of A, and H. LDA must be at
262 * least 1 and at least max( NN, 12 ). (12 is the
263 * dimension of the largest matrix on the precomputed
264 * input file.)
265 *
266 * H (workspace) COMPLEX*16 array, dimension (LDA, max(NN,12))
267 * Another copy of the test matrix A, modified by ZGEEVX.
268 *
269 * W (workspace) COMPLEX*16 array, dimension (max(NN,12))
270 * Contains the eigenvalues of A.
271 *
272 * W1 (workspace) COMPLEX*16 array, dimension (max(NN,12))
273 * Like W, this array contains the eigenvalues of A,
274 * but those computed when ZGEEVX only computes a partial
275 * eigendecomposition, i.e. not the eigenvalues and left
276 * and right eigenvectors.
277 *
278 * VL (workspace) COMPLEX*16 array, dimension (LDVL, max(NN,12))
279 * VL holds the computed left eigenvectors.
280 *
281 * LDVL (input) INTEGER
282 * Leading dimension of VL. Must be at least max(1,max(NN,12)).
283 *
284 * VR (workspace) COMPLEX*16 array, dimension (LDVR, max(NN,12))
285 * VR holds the computed right eigenvectors.
286 *
287 * LDVR (input) INTEGER
288 * Leading dimension of VR. Must be at least max(1,max(NN,12)).
289 *
290 * LRE (workspace) COMPLEX*16 array, dimension (LDLRE, max(NN,12))
291 * LRE holds the computed right or left eigenvectors.
292 *
293 * LDLRE (input) INTEGER
294 * Leading dimension of LRE. Must be at least max(1,max(NN,12))
295 *
296 * RESULT (output) DOUBLE PRECISION array, dimension (11)
297 * The values computed by the seven tests described above.
298 * The values are currently limited to 1/ulp, to avoid
299 * overflow.
300 *
301 * WORK (workspace) COMPLEX*16 array, dimension (NWORK)
302 *
303 * NWORK (input) INTEGER
304 * The number of entries in WORK. This must be at least
305 * max(6*12+2*12**2,6*NN(j)+2*NN(j)**2) =
306 * max( 360 ,6*NN(j)+2*NN(j)**2) for all j.
307 *
308 * RWORK (workspace) DOUBLE PRECISION array, dimension (2*max(NN,12))
309 *
310 * INFO (output) INTEGER
311 * If 0, then successful exit.
312 * If <0, then input paramter -INFO is incorrect.
313 * If >0, ZLATMR, CLATMS, CLATME or ZGET23 returned an error
314 * code, and INFO is its absolute value.
315 *
316 *-----------------------------------------------------------------------
317 *
318 * Some Local Variables and Parameters:
319 * ---- ----- --------- --- ----------
320 *
321 * ZERO, ONE Real 0 and 1.
322 * MAXTYP The number of types defined.
323 * NMAX Largest value in NN or 12.
324 * NERRS The number of tests which have exceeded THRESH
325 * COND, CONDS,
326 * IMODE Values to be passed to the matrix generators.
327 * ANORM Norm of A; passed to matrix generators.
328 *
329 * OVFL, UNFL Overflow and underflow thresholds.
330 * ULP, ULPINV Finest relative precision and its inverse.
331 * RTULP, RTULPI Square roots of the previous 4 values.
332 *
333 * The following four arrays decode JTYPE:
334 * KTYPE(j) The general type (1-10) for type "j".
335 * KMODE(j) The MODE value to be passed to the matrix
336 * generator for type "j".
337 * KMAGN(j) The order of magnitude ( O(1),
338 * O(overflow^(1/2) ), O(underflow^(1/2) )
339 * KCONDS(j) Selectw whether CONDS is to be 1 or
340 * 1/sqrt(ulp). (0 means irrelevant.)
341 *
342 * =====================================================================
343 *
344 * .. Parameters ..
345 COMPLEX*16 CZERO
346 PARAMETER ( CZERO = ( 0.0D+0, 0.0D+0 ) )
347 COMPLEX*16 CONE
348 PARAMETER ( CONE = ( 1.0D+0, 0.0D+0 ) )
349 DOUBLE PRECISION ZERO, ONE
350 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
351 INTEGER MAXTYP
352 PARAMETER ( MAXTYP = 21 )
353 * ..
354 * .. Local Scalars ..
355 LOGICAL BADNN
356 CHARACTER BALANC
357 CHARACTER*3 PATH
358 INTEGER I, IBAL, IINFO, IMODE, ISRT, ITYPE, IWK, J,
359 $ JCOL, JSIZE, JTYPE, MTYPES, N, NERRS, NFAIL,
360 $ NMAX, NNWORK, NTEST, NTESTF, NTESTT
361 DOUBLE PRECISION ANORM, COND, CONDS, OVFL, RTULP, RTULPI, ULP,
362 $ ULPINV, UNFL, WI, WR
363 * ..
364 * .. Local Arrays ..
365 CHARACTER BAL( 4 )
366 INTEGER IDUMMA( 1 ), IOLDSD( 4 ), KCONDS( MAXTYP ),
367 $ KMAGN( MAXTYP ), KMODE( MAXTYP ),
368 $ KTYPE( MAXTYP )
369 * ..
370 * .. External Functions ..
371 DOUBLE PRECISION DLAMCH
372 EXTERNAL DLAMCH
373 * ..
374 * .. External Subroutines ..
375 EXTERNAL DLABAD, DLASUM, XERBLA, ZGET23, ZLASET, ZLATME,
376 $ ZLATMR, ZLATMS
377 * ..
378 * .. Intrinsic Functions ..
379 INTRINSIC ABS, DCMPLX, MAX, MIN, SQRT
380 * ..
381 * .. Data statements ..
382 DATA KTYPE / 1, 2, 3, 5*4, 4*6, 6*6, 3*9 /
383 DATA KMAGN / 3*1, 1, 1, 1, 2, 3, 4*1, 1, 1, 1, 1, 2,
384 $ 3, 1, 2, 3 /
385 DATA KMODE / 3*0, 4, 3, 1, 4, 4, 4, 3, 1, 5, 4, 3,
386 $ 1, 5, 5, 5, 4, 3, 1 /
387 DATA KCONDS / 3*0, 5*0, 4*1, 6*2, 3*0 /
388 DATA BAL / 'N', 'P', 'S', 'B' /
389 * ..
390 * .. Executable Statements ..
391 *
392 PATH( 1: 1 ) = 'Zomplex precision'
393 PATH( 2: 3 ) = 'VX'
394 *
395 * Check for errors
396 *
397 NTESTT = 0
398 NTESTF = 0
399 INFO = 0
400 *
401 * Important constants
402 *
403 BADNN = .FALSE.
404 *
405 * 7 is the largest dimension in the input file of precomputed
406 * problems
407 *
408 NMAX = 7
409 DO 10 J = 1, NSIZES
410 NMAX = MAX( NMAX, NN( J ) )
411 IF( NN( J ).LT.0 )
412 $ BADNN = .TRUE.
413 10 CONTINUE
414 *
415 * Check for errors
416 *
417 IF( NSIZES.LT.0 ) THEN
418 INFO = -1
419 ELSE IF( BADNN ) THEN
420 INFO = -2
421 ELSE IF( NTYPES.LT.0 ) THEN
422 INFO = -3
423 ELSE IF( THRESH.LT.ZERO ) THEN
424 INFO = -6
425 ELSE IF( LDA.LT.1 .OR. LDA.LT.NMAX ) THEN
426 INFO = -10
427 ELSE IF( LDVL.LT.1 .OR. LDVL.LT.NMAX ) THEN
428 INFO = -15
429 ELSE IF( LDVR.LT.1 .OR. LDVR.LT.NMAX ) THEN
430 INFO = -17
431 ELSE IF( LDLRE.LT.1 .OR. LDLRE.LT.NMAX ) THEN
432 INFO = -19
433 ELSE IF( 6*NMAX+2*NMAX**2.GT.NWORK ) THEN
434 INFO = -30
435 END IF
436 *
437 IF( INFO.NE.0 ) THEN
438 CALL XERBLA( 'ZDRVVX', -INFO )
439 RETURN
440 END IF
441 *
442 * If nothing to do check on NIUNIT
443 *
444 IF( NSIZES.EQ.0 .OR. NTYPES.EQ.0 )
445 $ GO TO 160
446 *
447 * More Important constants
448 *
449 UNFL = DLAMCH( 'Safe minimum' )
450 OVFL = ONE / UNFL
451 CALL DLABAD( UNFL, OVFL )
452 ULP = DLAMCH( 'Precision' )
453 ULPINV = ONE / ULP
454 RTULP = SQRT( ULP )
455 RTULPI = ONE / RTULP
456 *
457 * Loop over sizes, types
458 *
459 NERRS = 0
460 *
461 DO 150 JSIZE = 1, NSIZES
462 N = NN( JSIZE )
463 IF( NSIZES.NE.1 ) THEN
464 MTYPES = MIN( MAXTYP, NTYPES )
465 ELSE
466 MTYPES = MIN( MAXTYP+1, NTYPES )
467 END IF
468 *
469 DO 140 JTYPE = 1, MTYPES
470 IF( .NOT.DOTYPE( JTYPE ) )
471 $ GO TO 140
472 *
473 * Save ISEED in case of an error.
474 *
475 DO 20 J = 1, 4
476 IOLDSD( J ) = ISEED( J )
477 20 CONTINUE
478 *
479 * Compute "A"
480 *
481 * Control parameters:
482 *
483 * KMAGN KCONDS KMODE KTYPE
484 * =1 O(1) 1 clustered 1 zero
485 * =2 large large clustered 2 identity
486 * =3 small exponential Jordan
487 * =4 arithmetic diagonal, (w/ eigenvalues)
488 * =5 random log symmetric, w/ eigenvalues
489 * =6 random general, w/ eigenvalues
490 * =7 random diagonal
491 * =8 random symmetric
492 * =9 random general
493 * =10 random triangular
494 *
495 IF( MTYPES.GT.MAXTYP )
496 $ GO TO 90
497 *
498 ITYPE = KTYPE( JTYPE )
499 IMODE = KMODE( JTYPE )
500 *
501 * Compute norm
502 *
503 GO TO ( 30, 40, 50 )KMAGN( JTYPE )
504 *
505 30 CONTINUE
506 ANORM = ONE
507 GO TO 60
508 *
509 40 CONTINUE
510 ANORM = OVFL*ULP
511 GO TO 60
512 *
513 50 CONTINUE
514 ANORM = UNFL*ULPINV
515 GO TO 60
516 *
517 60 CONTINUE
518 *
519 CALL ZLASET( 'Full', LDA, N, CZERO, CZERO, A, LDA )
520 IINFO = 0
521 COND = ULPINV
522 *
523 * Special Matrices -- Identity & Jordan block
524 *
525 * Zero
526 *
527 IF( ITYPE.EQ.1 ) THEN
528 IINFO = 0
529 *
530 ELSE IF( ITYPE.EQ.2 ) THEN
531 *
532 * Identity
533 *
534 DO 70 JCOL = 1, N
535 A( JCOL, JCOL ) = ANORM
536 70 CONTINUE
537 *
538 ELSE IF( ITYPE.EQ.3 ) THEN
539 *
540 * Jordan Block
541 *
542 DO 80 JCOL = 1, N
543 A( JCOL, JCOL ) = ANORM
544 IF( JCOL.GT.1 )
545 $ A( JCOL, JCOL-1 ) = ONE
546 80 CONTINUE
547 *
548 ELSE IF( ITYPE.EQ.4 ) THEN
549 *
550 * Diagonal Matrix, [Eigen]values Specified
551 *
552 CALL ZLATMS( N, N, 'S', ISEED, 'H', RWORK, IMODE, COND,
553 $ ANORM, 0, 0, 'N', A, LDA, WORK( N+1 ),
554 $ IINFO )
555 *
556 ELSE IF( ITYPE.EQ.5 ) THEN
557 *
558 * Symmetric, eigenvalues specified
559 *
560 CALL ZLATMS( N, N, 'S', ISEED, 'H', RWORK, IMODE, COND,
561 $ ANORM, N, N, 'N', A, LDA, WORK( N+1 ),
562 $ IINFO )
563 *
564 ELSE IF( ITYPE.EQ.6 ) THEN
565 *
566 * General, eigenvalues specified
567 *
568 IF( KCONDS( JTYPE ).EQ.1 ) THEN
569 CONDS = ONE
570 ELSE IF( KCONDS( JTYPE ).EQ.2 ) THEN
571 CONDS = RTULPI
572 ELSE
573 CONDS = ZERO
574 END IF
575 *
576 CALL ZLATME( N, 'D', ISEED, WORK, IMODE, COND, CONE, ' ',
577 $ 'T', 'T', 'T', RWORK, 4, CONDS, N, N, ANORM,
578 $ A, LDA, WORK( 2*N+1 ), IINFO )
579 *
580 ELSE IF( ITYPE.EQ.7 ) THEN
581 *
582 * Diagonal, random eigenvalues
583 *
584 CALL ZLATMR( N, N, 'D', ISEED, 'S', WORK, 6, ONE, CONE,
585 $ 'T', 'N', WORK( N+1 ), 1, ONE,
586 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, 0, 0,
587 $ ZERO, ANORM, 'NO', A, LDA, IDUMMA, IINFO )
588 *
589 ELSE IF( ITYPE.EQ.8 ) THEN
590 *
591 * Symmetric, random eigenvalues
592 *
593 CALL ZLATMR( N, N, 'D', ISEED, 'H', WORK, 6, ONE, CONE,
594 $ 'T', 'N', WORK( N+1 ), 1, ONE,
595 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N,
596 $ ZERO, ANORM, 'NO', A, LDA, IDUMMA, IINFO )
597 *
598 ELSE IF( ITYPE.EQ.9 ) THEN
599 *
600 * General, random eigenvalues
601 *
602 CALL ZLATMR( N, N, 'D', ISEED, 'N', WORK, 6, ONE, CONE,
603 $ 'T', 'N', WORK( N+1 ), 1, ONE,
604 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, N,
605 $ ZERO, ANORM, 'NO', A, LDA, IDUMMA, IINFO )
606 IF( N.GE.4 ) THEN
607 CALL ZLASET( 'Full', 2, N, CZERO, CZERO, A, LDA )
608 CALL ZLASET( 'Full', N-3, 1, CZERO, CZERO, A( 3, 1 ),
609 $ LDA )
610 CALL ZLASET( 'Full', N-3, 2, CZERO, CZERO,
611 $ A( 3, N-1 ), LDA )
612 CALL ZLASET( 'Full', 1, N, CZERO, CZERO, A( N, 1 ),
613 $ LDA )
614 END IF
615 *
616 ELSE IF( ITYPE.EQ.10 ) THEN
617 *
618 * Triangular, random eigenvalues
619 *
620 CALL ZLATMR( N, N, 'D', ISEED, 'N', WORK, 6, ONE, CONE,
621 $ 'T', 'N', WORK( N+1 ), 1, ONE,
622 $ WORK( 2*N+1 ), 1, ONE, 'N', IDUMMA, N, 0,
623 $ ZERO, ANORM, 'NO', A, LDA, IDUMMA, IINFO )
624 *
625 ELSE
626 *
627 IINFO = 1
628 END IF
629 *
630 IF( IINFO.NE.0 ) THEN
631 WRITE( NOUNIT, FMT = 9992 )'Generator', IINFO, N, JTYPE,
632 $ IOLDSD
633 INFO = ABS( IINFO )
634 RETURN
635 END IF
636 *
637 90 CONTINUE
638 *
639 * Test for minimal and generous workspace
640 *
641 DO 130 IWK = 1, 3
642 IF( IWK.EQ.1 ) THEN
643 NNWORK = 2*N
644 ELSE IF( IWK.EQ.2 ) THEN
645 NNWORK = 2*N + N**2
646 ELSE
647 NNWORK = 6*N + 2*N**2
648 END IF
649 NNWORK = MAX( NNWORK, 1 )
650 *
651 * Test for all balancing options
652 *
653 DO 120 IBAL = 1, 4
654 BALANC = BAL( IBAL )
655 *
656 * Perform tests
657 *
658 CALL ZGET23( .FALSE., 0, BALANC, JTYPE, THRESH,
659 $ IOLDSD, NOUNIT, N, A, LDA, H, W, W1, VL,
660 $ LDVL, VR, LDVR, LRE, LDLRE, RCONDV,
661 $ RCNDV1, RCDVIN, RCONDE, RCNDE1, RCDEIN,
662 $ SCALE, SCALE1, RESULT, WORK, NNWORK,
663 $ RWORK, INFO )
664 *
665 * Check for RESULT(j) > THRESH
666 *
667 NTEST = 0
668 NFAIL = 0
669 DO 100 J = 1, 9
670 IF( RESULT( J ).GE.ZERO )
671 $ NTEST = NTEST + 1
672 IF( RESULT( J ).GE.THRESH )
673 $ NFAIL = NFAIL + 1
674 100 CONTINUE
675 *
676 IF( NFAIL.GT.0 )
677 $ NTESTF = NTESTF + 1
678 IF( NTESTF.EQ.1 ) THEN
679 WRITE( NOUNIT, FMT = 9999 )PATH
680 WRITE( NOUNIT, FMT = 9998 )
681 WRITE( NOUNIT, FMT = 9997 )
682 WRITE( NOUNIT, FMT = 9996 )
683 WRITE( NOUNIT, FMT = 9995 )THRESH
684 NTESTF = 2
685 END IF
686 *
687 DO 110 J = 1, 9
688 IF( RESULT( J ).GE.THRESH ) THEN
689 WRITE( NOUNIT, FMT = 9994 )BALANC, N, IWK,
690 $ IOLDSD, JTYPE, J, RESULT( J )
691 END IF
692 110 CONTINUE
693 *
694 NERRS = NERRS + NFAIL
695 NTESTT = NTESTT + NTEST
696 *
697 120 CONTINUE
698 130 CONTINUE
699 140 CONTINUE
700 150 CONTINUE
701 *
702 160 CONTINUE
703 *
704 * Read in data from file to check accuracy of condition estimation.
705 * Assume input eigenvalues are sorted lexicographically (increasing
706 * by real part, then decreasing by imaginary part)
707 *
708 JTYPE = 0
709 170 CONTINUE
710 READ( NIUNIT, FMT = *, END = 220 )N, ISRT
711 *
712 * Read input data until N=0
713 *
714 IF( N.EQ.0 )
715 $ GO TO 220
716 JTYPE = JTYPE + 1
717 ISEED( 1 ) = JTYPE
718 DO 180 I = 1, N
719 READ( NIUNIT, FMT = * )( A( I, J ), J = 1, N )
720 180 CONTINUE
721 DO 190 I = 1, N
722 READ( NIUNIT, FMT = * )WR, WI, RCDEIN( I ), RCDVIN( I )
723 W1( I ) = DCMPLX( WR, WI )
724 190 CONTINUE
725 CALL ZGET23( .TRUE., ISRT, 'N', 22, THRESH, ISEED, NOUNIT, N, A,
726 $ LDA, H, W, W1, VL, LDVL, VR, LDVR, LRE, LDLRE,
727 $ RCONDV, RCNDV1, RCDVIN, RCONDE, RCNDE1, RCDEIN,
728 $ SCALE, SCALE1, RESULT, WORK, 6*N+2*N**2, RWORK,
729 $ INFO )
730 *
731 * Check for RESULT(j) > THRESH
732 *
733 NTEST = 0
734 NFAIL = 0
735 DO 200 J = 1, 11
736 IF( RESULT( J ).GE.ZERO )
737 $ NTEST = NTEST + 1
738 IF( RESULT( J ).GE.THRESH )
739 $ NFAIL = NFAIL + 1
740 200 CONTINUE
741 *
742 IF( NFAIL.GT.0 )
743 $ NTESTF = NTESTF + 1
744 IF( NTESTF.EQ.1 ) THEN
745 WRITE( NOUNIT, FMT = 9999 )PATH
746 WRITE( NOUNIT, FMT = 9998 )
747 WRITE( NOUNIT, FMT = 9997 )
748 WRITE( NOUNIT, FMT = 9996 )
749 WRITE( NOUNIT, FMT = 9995 )THRESH
750 NTESTF = 2
751 END IF
752 *
753 DO 210 J = 1, 11
754 IF( RESULT( J ).GE.THRESH ) THEN
755 WRITE( NOUNIT, FMT = 9993 )N, JTYPE, J, RESULT( J )
756 END IF
757 210 CONTINUE
758 *
759 NERRS = NERRS + NFAIL
760 NTESTT = NTESTT + NTEST
761 GO TO 170
762 220 CONTINUE
763 *
764 * Summary
765 *
766 CALL DLASUM( PATH, NOUNIT, NERRS, NTESTT )
767 *
768 9999 FORMAT( / 1X, A3, ' -- Complex Eigenvalue-Eigenvector ',
769 $ 'Decomposition Expert Driver',
770 $ / ' Matrix types (see ZDRVVX for details): ' )
771 *
772 9998 FORMAT( / ' Special Matrices:', / ' 1=Zero matrix. ',
773 $ ' ', ' 5=Diagonal: geometr. spaced entries.',
774 $ / ' 2=Identity matrix. ', ' 6=Diagona',
775 $ 'l: clustered entries.', / ' 3=Transposed Jordan block. ',
776 $ ' ', ' 7=Diagonal: large, evenly spaced.', / ' ',
777 $ '4=Diagonal: evenly spaced entries. ', ' 8=Diagonal: s',
778 $ 'mall, evenly spaced.' )
779 9997 FORMAT( ' Dense, Non-Symmetric Matrices:', / ' 9=Well-cond., ev',
780 $ 'enly spaced eigenvals.', ' 14=Ill-cond., geomet. spaced e',
781 $ 'igenals.', / ' 10=Well-cond., geom. spaced eigenvals. ',
782 $ ' 15=Ill-conditioned, clustered e.vals.', / ' 11=Well-cond',
783 $ 'itioned, clustered e.vals. ', ' 16=Ill-cond., random comp',
784 $ 'lex ', / ' 12=Well-cond., random complex ', ' ',
785 $ ' 17=Ill-cond., large rand. complx ', / ' 13=Ill-condi',
786 $ 'tioned, evenly spaced. ', ' 18=Ill-cond., small rand.',
787 $ ' complx ' )
788 9996 FORMAT( ' 19=Matrix with random O(1) entries. ', ' 21=Matrix ',
789 $ 'with small random entries.', / ' 20=Matrix with large ran',
790 $ 'dom entries. ', ' 22=Matrix read from input file', / )
791 9995 FORMAT( ' Tests performed with test threshold =', F8.2,
792 $ / / ' 1 = | A VR - VR W | / ( n |A| ulp ) ',
793 $ / ' 2 = | transpose(A) VL - VL W | / ( n |A| ulp ) ',
794 $ / ' 3 = | |VR(i)| - 1 | / ulp ',
795 $ / ' 4 = | |VL(i)| - 1 | / ulp ',
796 $ / ' 5 = 0 if W same no matter if VR or VL computed,',
797 $ ' 1/ulp otherwise', /
798 $ ' 6 = 0 if VR same no matter what else computed,',
799 $ ' 1/ulp otherwise', /
800 $ ' 7 = 0 if VL same no matter what else computed,',
801 $ ' 1/ulp otherwise', /
802 $ ' 8 = 0 if RCONDV same no matter what else computed,',
803 $ ' 1/ulp otherwise', /
804 $ ' 9 = 0 if SCALE, ILO, IHI, ABNRM same no matter what else',
805 $ ' computed, 1/ulp otherwise',
806 $ / ' 10 = | RCONDV - RCONDV(precomputed) | / cond(RCONDV),',
807 $ / ' 11 = | RCONDE - RCONDE(precomputed) | / cond(RCONDE),' )
808 9994 FORMAT( ' BALANC=''', A1, ''',N=', I4, ',IWK=', I1, ', seed=',
809 $ 4( I4, ',' ), ' type ', I2, ', test(', I2, ')=', G10.3 )
810 9993 FORMAT( ' N=', I5, ', input example =', I3, ', test(', I2, ')=',
811 $ G10.3 )
812 9992 FORMAT( ' ZDRVVX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
813 $ I6, ', JTYPE=', I6, ', ISEED=(', 3( I5, ',' ), I5, ')' )
814 *
815 RETURN
816 *
817 * End of ZDRVVX
818 *
819 END