1 SUBROUTINE DGET23( COMP, BALANC, JTYPE, THRESH, ISEED, NOUNIT, N,
2 $ A, LDA, H, WR, WI, WR1, WI1, VL, LDVL, VR,
3 $ LDVR, LRE, LDLRE, RCONDV, RCNDV1, RCDVIN,
4 $ RCONDE, RCNDE1, RCDEIN, SCALE, SCALE1, RESULT,
5 $ WORK, LWORK, IWORK, INFO )
6 *
7 * -- LAPACK test routine (version 3.1) --
8 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
9 * November 2006
10 *
11 * .. Scalar Arguments ..
12 LOGICAL COMP
13 CHARACTER BALANC
14 INTEGER INFO, JTYPE, LDA, LDLRE, LDVL, LDVR, LWORK, N,
15 $ NOUNIT
16 DOUBLE PRECISION THRESH
17 * ..
18 * .. Array Arguments ..
19 INTEGER ISEED( 4 ), IWORK( * )
20 DOUBLE PRECISION A( LDA, * ), H( LDA, * ), LRE( LDLRE, * ),
21 $ RCDEIN( * ), RCDVIN( * ), RCNDE1( * ),
22 $ RCNDV1( * ), RCONDE( * ), RCONDV( * ),
23 $ RESULT( 11 ), SCALE( * ), SCALE1( * ),
24 $ VL( LDVL, * ), VR( LDVR, * ), WI( * ),
25 $ WI1( * ), WORK( * ), WR( * ), WR1( * )
26 * ..
27 *
28 * Purpose
29 * =======
30 *
31 * DGET23 checks the nonsymmetric eigenvalue problem driver SGEEVX.
32 * If COMP = .FALSE., the first 8 of the following tests will be
33 * performed on the input matrix A, and also test 9 if LWORK is
34 * sufficiently large.
35 * if COMP is .TRUE. all 11 tests will be performed.
36 *
37 * (1) | A * VR - VR * W | / ( n |A| ulp )
38 *
39 * Here VR is the matrix of unit right eigenvectors.
40 * W is a block diagonal matrix, with a 1x1 block for each
41 * real eigenvalue and a 2x2 block for each complex conjugate
42 * pair. If eigenvalues j and j+1 are a complex conjugate pair,
43 * so WR(j) = WR(j+1) = wr and WI(j) = - WI(j+1) = wi, then the
44 * 2 x 2 block corresponding to the pair will be:
45 *
46 * ( wr wi )
47 * ( -wi wr )
48 *
49 * Such a block multiplying an n x 2 matrix ( ur ui ) on the
50 * right will be the same as multiplying ur + i*ui by wr + i*wi.
51 *
52 * (2) | A**H * VL - VL * W**H | / ( n |A| ulp )
53 *
54 * Here VL is the matrix of unit left eigenvectors, A**H is the
55 * conjugate transpose of A, and W is as above.
56 *
57 * (3) | |VR(i)| - 1 | / ulp and largest component real
58 *
59 * VR(i) denotes the i-th column of VR.
60 *
61 * (4) | |VL(i)| - 1 | / ulp and largest component real
62 *
63 * VL(i) denotes the i-th column of VL.
64 *
65 * (5) 0 if W(full) = W(partial), 1/ulp otherwise
66 *
67 * W(full) denotes the eigenvalues computed when VR, VL, RCONDV
68 * and RCONDE are also computed, and W(partial) denotes the
69 * eigenvalues computed when only some of VR, VL, RCONDV, and
70 * RCONDE are computed.
71 *
72 * (6) 0 if VR(full) = VR(partial), 1/ulp otherwise
73 *
74 * VR(full) denotes the right eigenvectors computed when VL, RCONDV
75 * and RCONDE are computed, and VR(partial) denotes the result
76 * when only some of VL and RCONDV are computed.
77 *
78 * (7) 0 if VL(full) = VL(partial), 1/ulp otherwise
79 *
80 * VL(full) denotes the left eigenvectors computed when VR, RCONDV
81 * and RCONDE are computed, and VL(partial) denotes the result
82 * when only some of VR and RCONDV are computed.
83 *
84 * (8) 0 if SCALE, ILO, IHI, ABNRM (full) =
85 * SCALE, ILO, IHI, ABNRM (partial)
86 * 1/ulp otherwise
87 *
88 * SCALE, ILO, IHI and ABNRM describe how the matrix is balanced.
89 * (full) is when VR, VL, RCONDE and RCONDV are also computed, and
90 * (partial) is when some are not computed.
91 *
92 * (9) 0 if RCONDV(full) = RCONDV(partial), 1/ulp otherwise
93 *
94 * RCONDV(full) denotes the reciprocal condition numbers of the
95 * right eigenvectors computed when VR, VL and RCONDE are also
96 * computed. RCONDV(partial) denotes the reciprocal condition
97 * numbers when only some of VR, VL and RCONDE are computed.
98 *
99 * (10) |RCONDV - RCDVIN| / cond(RCONDV)
100 *
101 * RCONDV is the reciprocal right eigenvector condition number
102 * computed by DGEEVX and RCDVIN (the precomputed true value)
103 * is supplied as input. cond(RCONDV) is the condition number of
104 * RCONDV, and takes errors in computing RCONDV into account, so
105 * that the resulting quantity should be O(ULP). cond(RCONDV) is
106 * essentially given by norm(A)/RCONDE.
107 *
108 * (11) |RCONDE - RCDEIN| / cond(RCONDE)
109 *
110 * RCONDE is the reciprocal eigenvalue condition number
111 * computed by DGEEVX and RCDEIN (the precomputed true value)
112 * is supplied as input. cond(RCONDE) is the condition number
113 * of RCONDE, and takes errors in computing RCONDE into account,
114 * so that the resulting quantity should be O(ULP). cond(RCONDE)
115 * is essentially given by norm(A)/RCONDV.
116 *
117 * Arguments
118 * =========
119 *
120 * COMP (input) LOGICAL
121 * COMP describes which input tests to perform:
122 * = .FALSE. if the computed condition numbers are not to
123 * be tested against RCDVIN and RCDEIN
124 * = .TRUE. if they are to be compared
125 *
126 * BALANC (input) CHARACTER
127 * Describes the balancing option to be tested.
128 * = 'N' for no permuting or diagonal scaling
129 * = 'P' for permuting but no diagonal scaling
130 * = 'S' for no permuting but diagonal scaling
131 * = 'B' for permuting and diagonal scaling
132 *
133 * JTYPE (input) INTEGER
134 * Type of input matrix. Used to label output if error occurs.
135 *
136 * THRESH (input) DOUBLE PRECISION
137 * A test will count as "failed" if the "error", computed as
138 * described above, exceeds THRESH. Note that the error
139 * is scaled to be O(1), so THRESH should be a reasonably
140 * small multiple of 1, e.g., 10 or 100. In particular,
141 * it should not depend on the precision (single vs. double)
142 * or the size of the matrix. It must be at least zero.
143 *
144 * ISEED (input) INTEGER array, dimension (4)
145 * If COMP = .FALSE., the random number generator seed
146 * used to produce matrix.
147 * If COMP = .TRUE., ISEED(1) = the number of the example.
148 * Used to label output if error occurs.
149 *
150 * NOUNIT (input) INTEGER
151 * The FORTRAN unit number for printing out error messages
152 * (e.g., if a routine returns INFO not equal to 0.)
153 *
154 * N (input) INTEGER
155 * The dimension of A. N must be at least 0.
156 *
157 * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
158 * Used to hold the matrix whose eigenvalues are to be
159 * computed.
160 *
161 * LDA (input) INTEGER
162 * The leading dimension of A, and H. LDA must be at
163 * least 1 and at least N.
164 *
165 * H (workspace) DOUBLE PRECISION array, dimension (LDA,N)
166 * Another copy of the test matrix A, modified by DGEEVX.
167 *
168 * WR (workspace) DOUBLE PRECISION array, dimension (N)
169 * WI (workspace) DOUBLE PRECISION array, dimension (N)
170 * The real and imaginary parts of the eigenvalues of A.
171 * On exit, WR + WI*i are the eigenvalues of the matrix in A.
172 *
173 * WR1 (workspace) DOUBLE PRECISION array, dimension (N)
174 * WI1 (workspace) DOUBLE PRECISION array, dimension (N)
175 * Like WR, WI, these arrays contain the eigenvalues of A,
176 * but those computed when DGEEVX only computes a partial
177 * eigendecomposition, i.e. not the eigenvalues and left
178 * and right eigenvectors.
179 *
180 * VL (workspace) DOUBLE PRECISION array, dimension (LDVL,N)
181 * VL holds the computed left eigenvectors.
182 *
183 * LDVL (input) INTEGER
184 * Leading dimension of VL. Must be at least max(1,N).
185 *
186 * VR (workspace) DOUBLE PRECISION array, dimension (LDVR,N)
187 * VR holds the computed right eigenvectors.
188 *
189 * LDVR (input) INTEGER
190 * Leading dimension of VR. Must be at least max(1,N).
191 *
192 * LRE (workspace) DOUBLE PRECISION array, dimension (LDLRE,N)
193 * LRE holds the computed right or left eigenvectors.
194 *
195 * LDLRE (input) INTEGER
196 * Leading dimension of LRE. Must be at least max(1,N).
197 *
198 * RCONDV (workspace) DOUBLE PRECISION array, dimension (N)
199 * RCONDV holds the computed reciprocal condition numbers
200 * for eigenvectors.
201 *
202 * RCNDV1 (workspace) DOUBLE PRECISION array, dimension (N)
203 * RCNDV1 holds more computed reciprocal condition numbers
204 * for eigenvectors.
205 *
206 * RCDVIN (input) DOUBLE PRECISION array, dimension (N)
207 * When COMP = .TRUE. RCDVIN holds the precomputed reciprocal
208 * condition numbers for eigenvectors to be compared with
209 * RCONDV.
210 *
211 * RCONDE (workspace) DOUBLE PRECISION array, dimension (N)
212 * RCONDE holds the computed reciprocal condition numbers
213 * for eigenvalues.
214 *
215 * RCNDE1 (workspace) DOUBLE PRECISION array, dimension (N)
216 * RCNDE1 holds more computed reciprocal condition numbers
217 * for eigenvalues.
218 *
219 * RCDEIN (input) DOUBLE PRECISION array, dimension (N)
220 * When COMP = .TRUE. RCDEIN holds the precomputed reciprocal
221 * condition numbers for eigenvalues to be compared with
222 * RCONDE.
223 *
224 * SCALE (workspace) DOUBLE PRECISION array, dimension (N)
225 * Holds information describing balancing of matrix.
226 *
227 * SCALE1 (workspace) DOUBLE PRECISION array, dimension (N)
228 * Holds information describing balancing of matrix.
229 *
230 * RESULT (output) DOUBLE PRECISION array, dimension (11)
231 * The values computed by the 11 tests described above.
232 * The values are currently limited to 1/ulp, to avoid
233 * overflow.
234 *
235 * WORK (workspace) DOUBLE PRECISION array, dimension (LWORK)
236 *
237 * LWORK (input) INTEGER
238 * The number of entries in WORK. This must be at least
239 * 3*N, and 6*N+N**2 if tests 9, 10 or 11 are to be performed.
240 *
241 * IWORK (workspace) INTEGER array, dimension (2*N)
242 *
243 * INFO (output) INTEGER
244 * If 0, successful exit.
245 * If <0, input parameter -INFO had an incorrect value.
246 * If >0, DGEEVX returned an error code, the absolute
247 * value of which is returned.
248 *
249 * =====================================================================
250 *
251 *
252 * .. Parameters ..
253 DOUBLE PRECISION ZERO, ONE, TWO
254 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
255 DOUBLE PRECISION EPSIN
256 PARAMETER ( EPSIN = 5.9605D-8 )
257 * ..
258 * .. Local Scalars ..
259 LOGICAL BALOK, NOBAL
260 CHARACTER SENSE
261 INTEGER I, IHI, IHI1, IINFO, ILO, ILO1, ISENS, ISENSM,
262 $ J, JJ, KMIN
263 DOUBLE PRECISION ABNRM, ABNRM1, EPS, SMLNUM, TNRM, TOL, TOLIN,
264 $ ULP, ULPINV, V, VIMIN, VMAX, VMX, VRMIN, VRMX,
265 $ VTST
266 * ..
267 * .. Local Arrays ..
268 CHARACTER SENS( 2 )
269 DOUBLE PRECISION DUM( 1 ), RES( 2 )
270 * ..
271 * .. External Functions ..
272 LOGICAL LSAME
273 DOUBLE PRECISION DLAMCH, DLAPY2, DNRM2
274 EXTERNAL LSAME, DLAMCH, DLAPY2, DNRM2
275 * ..
276 * .. External Subroutines ..
277 EXTERNAL DGEEVX, DGET22, DLACPY, XERBLA
278 * ..
279 * .. Intrinsic Functions ..
280 INTRINSIC ABS, DBLE, MAX, MIN
281 * ..
282 * .. Data statements ..
283 DATA SENS / 'N', 'V' /
284 * ..
285 * .. Executable Statements ..
286 *
287 * Check for errors
288 *
289 NOBAL = LSAME( BALANC, 'N' )
290 BALOK = NOBAL .OR. LSAME( BALANC, 'P' ) .OR.
291 $ LSAME( BALANC, 'S' ) .OR. LSAME( BALANC, 'B' )
292 INFO = 0
293 IF( .NOT.BALOK ) THEN
294 INFO = -2
295 ELSE IF( THRESH.LT.ZERO ) THEN
296 INFO = -4
297 ELSE IF( NOUNIT.LE.0 ) THEN
298 INFO = -6
299 ELSE IF( N.LT.0 ) THEN
300 INFO = -7
301 ELSE IF( LDA.LT.1 .OR. LDA.LT.N ) THEN
302 INFO = -9
303 ELSE IF( LDVL.LT.1 .OR. LDVL.LT.N ) THEN
304 INFO = -16
305 ELSE IF( LDVR.LT.1 .OR. LDVR.LT.N ) THEN
306 INFO = -18
307 ELSE IF( LDLRE.LT.1 .OR. LDLRE.LT.N ) THEN
308 INFO = -20
309 ELSE IF( LWORK.LT.3*N .OR. ( COMP .AND. LWORK.LT.6*N+N*N ) ) THEN
310 INFO = -31
311 END IF
312 *
313 IF( INFO.NE.0 ) THEN
314 CALL XERBLA( 'DGET23', -INFO )
315 RETURN
316 END IF
317 *
318 * Quick return if nothing to do
319 *
320 DO 10 I = 1, 11
321 RESULT( I ) = -ONE
322 10 CONTINUE
323 *
324 IF( N.EQ.0 )
325 $ RETURN
326 *
327 * More Important constants
328 *
329 ULP = DLAMCH( 'Precision' )
330 SMLNUM = DLAMCH( 'S' )
331 ULPINV = ONE / ULP
332 *
333 * Compute eigenvalues and eigenvectors, and test them
334 *
335 IF( LWORK.GE.6*N+N*N ) THEN
336 SENSE = 'B'
337 ISENSM = 2
338 ELSE
339 SENSE = 'E'
340 ISENSM = 1
341 END IF
342 CALL DLACPY( 'F', N, N, A, LDA, H, LDA )
343 CALL DGEEVX( BALANC, 'V', 'V', SENSE, N, H, LDA, WR, WI, VL, LDVL,
344 $ VR, LDVR, ILO, IHI, SCALE, ABNRM, RCONDE, RCONDV,
345 $ WORK, LWORK, IWORK, IINFO )
346 IF( IINFO.NE.0 ) THEN
347 RESULT( 1 ) = ULPINV
348 IF( JTYPE.NE.22 ) THEN
349 WRITE( NOUNIT, FMT = 9998 )'DGEEVX1', IINFO, N, JTYPE,
350 $ BALANC, ISEED
351 ELSE
352 WRITE( NOUNIT, FMT = 9999 )'DGEEVX1', IINFO, N, ISEED( 1 )
353 END IF
354 INFO = ABS( IINFO )
355 RETURN
356 END IF
357 *
358 * Do Test (1)
359 *
360 CALL DGET22( 'N', 'N', 'N', N, A, LDA, VR, LDVR, WR, WI, WORK,
361 $ RES )
362 RESULT( 1 ) = RES( 1 )
363 *
364 * Do Test (2)
365 *
366 CALL DGET22( 'T', 'N', 'T', N, A, LDA, VL, LDVL, WR, WI, WORK,
367 $ RES )
368 RESULT( 2 ) = RES( 1 )
369 *
370 * Do Test (3)
371 *
372 DO 30 J = 1, N
373 TNRM = ONE
374 IF( WI( J ).EQ.ZERO ) THEN
375 TNRM = DNRM2( N, VR( 1, J ), 1 )
376 ELSE IF( WI( J ).GT.ZERO ) THEN
377 TNRM = DLAPY2( DNRM2( N, VR( 1, J ), 1 ),
378 $ DNRM2( N, VR( 1, J+1 ), 1 ) )
379 END IF
380 RESULT( 3 ) = MAX( RESULT( 3 ),
381 $ MIN( ULPINV, ABS( TNRM-ONE ) / ULP ) )
382 IF( WI( J ).GT.ZERO ) THEN
383 VMX = ZERO
384 VRMX = ZERO
385 DO 20 JJ = 1, N
386 VTST = DLAPY2( VR( JJ, J ), VR( JJ, J+1 ) )
387 IF( VTST.GT.VMX )
388 $ VMX = VTST
389 IF( VR( JJ, J+1 ).EQ.ZERO .AND. ABS( VR( JJ, J ) ).GT.
390 $ VRMX )VRMX = ABS( VR( JJ, J ) )
391 20 CONTINUE
392 IF( VRMX / VMX.LT.ONE-TWO*ULP )
393 $ RESULT( 3 ) = ULPINV
394 END IF
395 30 CONTINUE
396 *
397 * Do Test (4)
398 *
399 DO 50 J = 1, N
400 TNRM = ONE
401 IF( WI( J ).EQ.ZERO ) THEN
402 TNRM = DNRM2( N, VL( 1, J ), 1 )
403 ELSE IF( WI( J ).GT.ZERO ) THEN
404 TNRM = DLAPY2( DNRM2( N, VL( 1, J ), 1 ),
405 $ DNRM2( N, VL( 1, J+1 ), 1 ) )
406 END IF
407 RESULT( 4 ) = MAX( RESULT( 4 ),
408 $ MIN( ULPINV, ABS( TNRM-ONE ) / ULP ) )
409 IF( WI( J ).GT.ZERO ) THEN
410 VMX = ZERO
411 VRMX = ZERO
412 DO 40 JJ = 1, N
413 VTST = DLAPY2( VL( JJ, J ), VL( JJ, J+1 ) )
414 IF( VTST.GT.VMX )
415 $ VMX = VTST
416 IF( VL( JJ, J+1 ).EQ.ZERO .AND. ABS( VL( JJ, J ) ).GT.
417 $ VRMX )VRMX = ABS( VL( JJ, J ) )
418 40 CONTINUE
419 IF( VRMX / VMX.LT.ONE-TWO*ULP )
420 $ RESULT( 4 ) = ULPINV
421 END IF
422 50 CONTINUE
423 *
424 * Test for all options of computing condition numbers
425 *
426 DO 200 ISENS = 1, ISENSM
427 *
428 SENSE = SENS( ISENS )
429 *
430 * Compute eigenvalues only, and test them
431 *
432 CALL DLACPY( 'F', N, N, A, LDA, H, LDA )
433 CALL DGEEVX( BALANC, 'N', 'N', SENSE, N, H, LDA, WR1, WI1, DUM,
434 $ 1, DUM, 1, ILO1, IHI1, SCALE1, ABNRM1, RCNDE1,
435 $ RCNDV1, WORK, LWORK, IWORK, IINFO )
436 IF( IINFO.NE.0 ) THEN
437 RESULT( 1 ) = ULPINV
438 IF( JTYPE.NE.22 ) THEN
439 WRITE( NOUNIT, FMT = 9998 )'DGEEVX2', IINFO, N, JTYPE,
440 $ BALANC, ISEED
441 ELSE
442 WRITE( NOUNIT, FMT = 9999 )'DGEEVX2', IINFO, N,
443 $ ISEED( 1 )
444 END IF
445 INFO = ABS( IINFO )
446 GO TO 190
447 END IF
448 *
449 * Do Test (5)
450 *
451 DO 60 J = 1, N
452 IF( WR( J ).NE.WR1( J ) .OR. WI( J ).NE.WI1( J ) )
453 $ RESULT( 5 ) = ULPINV
454 60 CONTINUE
455 *
456 * Do Test (8)
457 *
458 IF( .NOT.NOBAL ) THEN
459 DO 70 J = 1, N
460 IF( SCALE( J ).NE.SCALE1( J ) )
461 $ RESULT( 8 ) = ULPINV
462 70 CONTINUE
463 IF( ILO.NE.ILO1 )
464 $ RESULT( 8 ) = ULPINV
465 IF( IHI.NE.IHI1 )
466 $ RESULT( 8 ) = ULPINV
467 IF( ABNRM.NE.ABNRM1 )
468 $ RESULT( 8 ) = ULPINV
469 END IF
470 *
471 * Do Test (9)
472 *
473 IF( ISENS.EQ.2 .AND. N.GT.1 ) THEN
474 DO 80 J = 1, N
475 IF( RCONDV( J ).NE.RCNDV1( J ) )
476 $ RESULT( 9 ) = ULPINV
477 80 CONTINUE
478 END IF
479 *
480 * Compute eigenvalues and right eigenvectors, and test them
481 *
482 CALL DLACPY( 'F', N, N, A, LDA, H, LDA )
483 CALL DGEEVX( BALANC, 'N', 'V', SENSE, N, H, LDA, WR1, WI1, DUM,
484 $ 1, LRE, LDLRE, ILO1, IHI1, SCALE1, ABNRM1, RCNDE1,
485 $ RCNDV1, WORK, LWORK, IWORK, IINFO )
486 IF( IINFO.NE.0 ) THEN
487 RESULT( 1 ) = ULPINV
488 IF( JTYPE.NE.22 ) THEN
489 WRITE( NOUNIT, FMT = 9998 )'DGEEVX3', IINFO, N, JTYPE,
490 $ BALANC, ISEED
491 ELSE
492 WRITE( NOUNIT, FMT = 9999 )'DGEEVX3', IINFO, N,
493 $ ISEED( 1 )
494 END IF
495 INFO = ABS( IINFO )
496 GO TO 190
497 END IF
498 *
499 * Do Test (5) again
500 *
501 DO 90 J = 1, N
502 IF( WR( J ).NE.WR1( J ) .OR. WI( J ).NE.WI1( J ) )
503 $ RESULT( 5 ) = ULPINV
504 90 CONTINUE
505 *
506 * Do Test (6)
507 *
508 DO 110 J = 1, N
509 DO 100 JJ = 1, N
510 IF( VR( J, JJ ).NE.LRE( J, JJ ) )
511 $ RESULT( 6 ) = ULPINV
512 100 CONTINUE
513 110 CONTINUE
514 *
515 * Do Test (8) again
516 *
517 IF( .NOT.NOBAL ) THEN
518 DO 120 J = 1, N
519 IF( SCALE( J ).NE.SCALE1( J ) )
520 $ RESULT( 8 ) = ULPINV
521 120 CONTINUE
522 IF( ILO.NE.ILO1 )
523 $ RESULT( 8 ) = ULPINV
524 IF( IHI.NE.IHI1 )
525 $ RESULT( 8 ) = ULPINV
526 IF( ABNRM.NE.ABNRM1 )
527 $ RESULT( 8 ) = ULPINV
528 END IF
529 *
530 * Do Test (9) again
531 *
532 IF( ISENS.EQ.2 .AND. N.GT.1 ) THEN
533 DO 130 J = 1, N
534 IF( RCONDV( J ).NE.RCNDV1( J ) )
535 $ RESULT( 9 ) = ULPINV
536 130 CONTINUE
537 END IF
538 *
539 * Compute eigenvalues and left eigenvectors, and test them
540 *
541 CALL DLACPY( 'F', N, N, A, LDA, H, LDA )
542 CALL DGEEVX( BALANC, 'V', 'N', SENSE, N, H, LDA, WR1, WI1, LRE,
543 $ LDLRE, DUM, 1, ILO1, IHI1, SCALE1, ABNRM1, RCNDE1,
544 $ RCNDV1, WORK, LWORK, IWORK, IINFO )
545 IF( IINFO.NE.0 ) THEN
546 RESULT( 1 ) = ULPINV
547 IF( JTYPE.NE.22 ) THEN
548 WRITE( NOUNIT, FMT = 9998 )'DGEEVX4', IINFO, N, JTYPE,
549 $ BALANC, ISEED
550 ELSE
551 WRITE( NOUNIT, FMT = 9999 )'DGEEVX4', IINFO, N,
552 $ ISEED( 1 )
553 END IF
554 INFO = ABS( IINFO )
555 GO TO 190
556 END IF
557 *
558 * Do Test (5) again
559 *
560 DO 140 J = 1, N
561 IF( WR( J ).NE.WR1( J ) .OR. WI( J ).NE.WI1( J ) )
562 $ RESULT( 5 ) = ULPINV
563 140 CONTINUE
564 *
565 * Do Test (7)
566 *
567 DO 160 J = 1, N
568 DO 150 JJ = 1, N
569 IF( VL( J, JJ ).NE.LRE( J, JJ ) )
570 $ RESULT( 7 ) = ULPINV
571 150 CONTINUE
572 160 CONTINUE
573 *
574 * Do Test (8) again
575 *
576 IF( .NOT.NOBAL ) THEN
577 DO 170 J = 1, N
578 IF( SCALE( J ).NE.SCALE1( J ) )
579 $ RESULT( 8 ) = ULPINV
580 170 CONTINUE
581 IF( ILO.NE.ILO1 )
582 $ RESULT( 8 ) = ULPINV
583 IF( IHI.NE.IHI1 )
584 $ RESULT( 8 ) = ULPINV
585 IF( ABNRM.NE.ABNRM1 )
586 $ RESULT( 8 ) = ULPINV
587 END IF
588 *
589 * Do Test (9) again
590 *
591 IF( ISENS.EQ.2 .AND. N.GT.1 ) THEN
592 DO 180 J = 1, N
593 IF( RCONDV( J ).NE.RCNDV1( J ) )
594 $ RESULT( 9 ) = ULPINV
595 180 CONTINUE
596 END IF
597 *
598 190 CONTINUE
599 *
600 200 CONTINUE
601 *
602 * If COMP, compare condition numbers to precomputed ones
603 *
604 IF( COMP ) THEN
605 CALL DLACPY( 'F', N, N, A, LDA, H, LDA )
606 CALL DGEEVX( 'N', 'V', 'V', 'B', N, H, LDA, WR, WI, VL, LDVL,
607 $ VR, LDVR, ILO, IHI, SCALE, ABNRM, RCONDE, RCONDV,
608 $ WORK, LWORK, IWORK, IINFO )
609 IF( IINFO.NE.0 ) THEN
610 RESULT( 1 ) = ULPINV
611 WRITE( NOUNIT, FMT = 9999 )'DGEEVX5', IINFO, N, ISEED( 1 )
612 INFO = ABS( IINFO )
613 GO TO 250
614 END IF
615 *
616 * Sort eigenvalues and condition numbers lexicographically
617 * to compare with inputs
618 *
619 DO 220 I = 1, N - 1
620 KMIN = I
621 VRMIN = WR( I )
622 VIMIN = WI( I )
623 DO 210 J = I + 1, N
624 IF( WR( J ).LT.VRMIN ) THEN
625 KMIN = J
626 VRMIN = WR( J )
627 VIMIN = WI( J )
628 END IF
629 210 CONTINUE
630 WR( KMIN ) = WR( I )
631 WI( KMIN ) = WI( I )
632 WR( I ) = VRMIN
633 WI( I ) = VIMIN
634 VRMIN = RCONDE( KMIN )
635 RCONDE( KMIN ) = RCONDE( I )
636 RCONDE( I ) = VRMIN
637 VRMIN = RCONDV( KMIN )
638 RCONDV( KMIN ) = RCONDV( I )
639 RCONDV( I ) = VRMIN
640 220 CONTINUE
641 *
642 * Compare condition numbers for eigenvectors
643 * taking their condition numbers into account
644 *
645 RESULT( 10 ) = ZERO
646 EPS = MAX( EPSIN, ULP )
647 V = MAX( DBLE( N )*EPS*ABNRM, SMLNUM )
648 IF( ABNRM.EQ.ZERO )
649 $ V = ONE
650 DO 230 I = 1, N
651 IF( V.GT.RCONDV( I )*RCONDE( I ) ) THEN
652 TOL = RCONDV( I )
653 ELSE
654 TOL = V / RCONDE( I )
655 END IF
656 IF( V.GT.RCDVIN( I )*RCDEIN( I ) ) THEN
657 TOLIN = RCDVIN( I )
658 ELSE
659 TOLIN = V / RCDEIN( I )
660 END IF
661 TOL = MAX( TOL, SMLNUM / EPS )
662 TOLIN = MAX( TOLIN, SMLNUM / EPS )
663 IF( EPS*( RCDVIN( I )-TOLIN ).GT.RCONDV( I )+TOL ) THEN
664 VMAX = ONE / EPS
665 ELSE IF( RCDVIN( I )-TOLIN.GT.RCONDV( I )+TOL ) THEN
666 VMAX = ( RCDVIN( I )-TOLIN ) / ( RCONDV( I )+TOL )
667 ELSE IF( RCDVIN( I )+TOLIN.LT.EPS*( RCONDV( I )-TOL ) ) THEN
668 VMAX = ONE / EPS
669 ELSE IF( RCDVIN( I )+TOLIN.LT.RCONDV( I )-TOL ) THEN
670 VMAX = ( RCONDV( I )-TOL ) / ( RCDVIN( I )+TOLIN )
671 ELSE
672 VMAX = ONE
673 END IF
674 RESULT( 10 ) = MAX( RESULT( 10 ), VMAX )
675 230 CONTINUE
676 *
677 * Compare condition numbers for eigenvalues
678 * taking their condition numbers into account
679 *
680 RESULT( 11 ) = ZERO
681 DO 240 I = 1, N
682 IF( V.GT.RCONDV( I ) ) THEN
683 TOL = ONE
684 ELSE
685 TOL = V / RCONDV( I )
686 END IF
687 IF( V.GT.RCDVIN( I ) ) THEN
688 TOLIN = ONE
689 ELSE
690 TOLIN = V / RCDVIN( I )
691 END IF
692 TOL = MAX( TOL, SMLNUM / EPS )
693 TOLIN = MAX( TOLIN, SMLNUM / EPS )
694 IF( EPS*( RCDEIN( I )-TOLIN ).GT.RCONDE( I )+TOL ) THEN
695 VMAX = ONE / EPS
696 ELSE IF( RCDEIN( I )-TOLIN.GT.RCONDE( I )+TOL ) THEN
697 VMAX = ( RCDEIN( I )-TOLIN ) / ( RCONDE( I )+TOL )
698 ELSE IF( RCDEIN( I )+TOLIN.LT.EPS*( RCONDE( I )-TOL ) ) THEN
699 VMAX = ONE / EPS
700 ELSE IF( RCDEIN( I )+TOLIN.LT.RCONDE( I )-TOL ) THEN
701 VMAX = ( RCONDE( I )-TOL ) / ( RCDEIN( I )+TOLIN )
702 ELSE
703 VMAX = ONE
704 END IF
705 RESULT( 11 ) = MAX( RESULT( 11 ), VMAX )
706 240 CONTINUE
707 250 CONTINUE
708 *
709 END IF
710 *
711 9999 FORMAT( ' DGET23: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
712 $ I6, ', INPUT EXAMPLE NUMBER = ', I4 )
713 9998 FORMAT( ' DGET23: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
714 $ I6, ', JTYPE=', I6, ', BALANC = ', A, ', ISEED=(',
715 $ 3( I5, ',' ), I5, ')' )
716 *
717 RETURN
718 *
719 * End of DGET23
720 *
721 END
2 $ A, LDA, H, WR, WI, WR1, WI1, VL, LDVL, VR,
3 $ LDVR, LRE, LDLRE, RCONDV, RCNDV1, RCDVIN,
4 $ RCONDE, RCNDE1, RCDEIN, SCALE, SCALE1, RESULT,
5 $ WORK, LWORK, IWORK, INFO )
6 *
7 * -- LAPACK test routine (version 3.1) --
8 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
9 * November 2006
10 *
11 * .. Scalar Arguments ..
12 LOGICAL COMP
13 CHARACTER BALANC
14 INTEGER INFO, JTYPE, LDA, LDLRE, LDVL, LDVR, LWORK, N,
15 $ NOUNIT
16 DOUBLE PRECISION THRESH
17 * ..
18 * .. Array Arguments ..
19 INTEGER ISEED( 4 ), IWORK( * )
20 DOUBLE PRECISION A( LDA, * ), H( LDA, * ), LRE( LDLRE, * ),
21 $ RCDEIN( * ), RCDVIN( * ), RCNDE1( * ),
22 $ RCNDV1( * ), RCONDE( * ), RCONDV( * ),
23 $ RESULT( 11 ), SCALE( * ), SCALE1( * ),
24 $ VL( LDVL, * ), VR( LDVR, * ), WI( * ),
25 $ WI1( * ), WORK( * ), WR( * ), WR1( * )
26 * ..
27 *
28 * Purpose
29 * =======
30 *
31 * DGET23 checks the nonsymmetric eigenvalue problem driver SGEEVX.
32 * If COMP = .FALSE., the first 8 of the following tests will be
33 * performed on the input matrix A, and also test 9 if LWORK is
34 * sufficiently large.
35 * if COMP is .TRUE. all 11 tests will be performed.
36 *
37 * (1) | A * VR - VR * W | / ( n |A| ulp )
38 *
39 * Here VR is the matrix of unit right eigenvectors.
40 * W is a block diagonal matrix, with a 1x1 block for each
41 * real eigenvalue and a 2x2 block for each complex conjugate
42 * pair. If eigenvalues j and j+1 are a complex conjugate pair,
43 * so WR(j) = WR(j+1) = wr and WI(j) = - WI(j+1) = wi, then the
44 * 2 x 2 block corresponding to the pair will be:
45 *
46 * ( wr wi )
47 * ( -wi wr )
48 *
49 * Such a block multiplying an n x 2 matrix ( ur ui ) on the
50 * right will be the same as multiplying ur + i*ui by wr + i*wi.
51 *
52 * (2) | A**H * VL - VL * W**H | / ( n |A| ulp )
53 *
54 * Here VL is the matrix of unit left eigenvectors, A**H is the
55 * conjugate transpose of A, and W is as above.
56 *
57 * (3) | |VR(i)| - 1 | / ulp and largest component real
58 *
59 * VR(i) denotes the i-th column of VR.
60 *
61 * (4) | |VL(i)| - 1 | / ulp and largest component real
62 *
63 * VL(i) denotes the i-th column of VL.
64 *
65 * (5) 0 if W(full) = W(partial), 1/ulp otherwise
66 *
67 * W(full) denotes the eigenvalues computed when VR, VL, RCONDV
68 * and RCONDE are also computed, and W(partial) denotes the
69 * eigenvalues computed when only some of VR, VL, RCONDV, and
70 * RCONDE are computed.
71 *
72 * (6) 0 if VR(full) = VR(partial), 1/ulp otherwise
73 *
74 * VR(full) denotes the right eigenvectors computed when VL, RCONDV
75 * and RCONDE are computed, and VR(partial) denotes the result
76 * when only some of VL and RCONDV are computed.
77 *
78 * (7) 0 if VL(full) = VL(partial), 1/ulp otherwise
79 *
80 * VL(full) denotes the left eigenvectors computed when VR, RCONDV
81 * and RCONDE are computed, and VL(partial) denotes the result
82 * when only some of VR and RCONDV are computed.
83 *
84 * (8) 0 if SCALE, ILO, IHI, ABNRM (full) =
85 * SCALE, ILO, IHI, ABNRM (partial)
86 * 1/ulp otherwise
87 *
88 * SCALE, ILO, IHI and ABNRM describe how the matrix is balanced.
89 * (full) is when VR, VL, RCONDE and RCONDV are also computed, and
90 * (partial) is when some are not computed.
91 *
92 * (9) 0 if RCONDV(full) = RCONDV(partial), 1/ulp otherwise
93 *
94 * RCONDV(full) denotes the reciprocal condition numbers of the
95 * right eigenvectors computed when VR, VL and RCONDE are also
96 * computed. RCONDV(partial) denotes the reciprocal condition
97 * numbers when only some of VR, VL and RCONDE are computed.
98 *
99 * (10) |RCONDV - RCDVIN| / cond(RCONDV)
100 *
101 * RCONDV is the reciprocal right eigenvector condition number
102 * computed by DGEEVX and RCDVIN (the precomputed true value)
103 * is supplied as input. cond(RCONDV) is the condition number of
104 * RCONDV, and takes errors in computing RCONDV into account, so
105 * that the resulting quantity should be O(ULP). cond(RCONDV) is
106 * essentially given by norm(A)/RCONDE.
107 *
108 * (11) |RCONDE - RCDEIN| / cond(RCONDE)
109 *
110 * RCONDE is the reciprocal eigenvalue condition number
111 * computed by DGEEVX and RCDEIN (the precomputed true value)
112 * is supplied as input. cond(RCONDE) is the condition number
113 * of RCONDE, and takes errors in computing RCONDE into account,
114 * so that the resulting quantity should be O(ULP). cond(RCONDE)
115 * is essentially given by norm(A)/RCONDV.
116 *
117 * Arguments
118 * =========
119 *
120 * COMP (input) LOGICAL
121 * COMP describes which input tests to perform:
122 * = .FALSE. if the computed condition numbers are not to
123 * be tested against RCDVIN and RCDEIN
124 * = .TRUE. if they are to be compared
125 *
126 * BALANC (input) CHARACTER
127 * Describes the balancing option to be tested.
128 * = 'N' for no permuting or diagonal scaling
129 * = 'P' for permuting but no diagonal scaling
130 * = 'S' for no permuting but diagonal scaling
131 * = 'B' for permuting and diagonal scaling
132 *
133 * JTYPE (input) INTEGER
134 * Type of input matrix. Used to label output if error occurs.
135 *
136 * THRESH (input) DOUBLE PRECISION
137 * A test will count as "failed" if the "error", computed as
138 * described above, exceeds THRESH. Note that the error
139 * is scaled to be O(1), so THRESH should be a reasonably
140 * small multiple of 1, e.g., 10 or 100. In particular,
141 * it should not depend on the precision (single vs. double)
142 * or the size of the matrix. It must be at least zero.
143 *
144 * ISEED (input) INTEGER array, dimension (4)
145 * If COMP = .FALSE., the random number generator seed
146 * used to produce matrix.
147 * If COMP = .TRUE., ISEED(1) = the number of the example.
148 * Used to label output if error occurs.
149 *
150 * NOUNIT (input) INTEGER
151 * The FORTRAN unit number for printing out error messages
152 * (e.g., if a routine returns INFO not equal to 0.)
153 *
154 * N (input) INTEGER
155 * The dimension of A. N must be at least 0.
156 *
157 * A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
158 * Used to hold the matrix whose eigenvalues are to be
159 * computed.
160 *
161 * LDA (input) INTEGER
162 * The leading dimension of A, and H. LDA must be at
163 * least 1 and at least N.
164 *
165 * H (workspace) DOUBLE PRECISION array, dimension (LDA,N)
166 * Another copy of the test matrix A, modified by DGEEVX.
167 *
168 * WR (workspace) DOUBLE PRECISION array, dimension (N)
169 * WI (workspace) DOUBLE PRECISION array, dimension (N)
170 * The real and imaginary parts of the eigenvalues of A.
171 * On exit, WR + WI*i are the eigenvalues of the matrix in A.
172 *
173 * WR1 (workspace) DOUBLE PRECISION array, dimension (N)
174 * WI1 (workspace) DOUBLE PRECISION array, dimension (N)
175 * Like WR, WI, these arrays contain the eigenvalues of A,
176 * but those computed when DGEEVX only computes a partial
177 * eigendecomposition, i.e. not the eigenvalues and left
178 * and right eigenvectors.
179 *
180 * VL (workspace) DOUBLE PRECISION array, dimension (LDVL,N)
181 * VL holds the computed left eigenvectors.
182 *
183 * LDVL (input) INTEGER
184 * Leading dimension of VL. Must be at least max(1,N).
185 *
186 * VR (workspace) DOUBLE PRECISION array, dimension (LDVR,N)
187 * VR holds the computed right eigenvectors.
188 *
189 * LDVR (input) INTEGER
190 * Leading dimension of VR. Must be at least max(1,N).
191 *
192 * LRE (workspace) DOUBLE PRECISION array, dimension (LDLRE,N)
193 * LRE holds the computed right or left eigenvectors.
194 *
195 * LDLRE (input) INTEGER
196 * Leading dimension of LRE. Must be at least max(1,N).
197 *
198 * RCONDV (workspace) DOUBLE PRECISION array, dimension (N)
199 * RCONDV holds the computed reciprocal condition numbers
200 * for eigenvectors.
201 *
202 * RCNDV1 (workspace) DOUBLE PRECISION array, dimension (N)
203 * RCNDV1 holds more computed reciprocal condition numbers
204 * for eigenvectors.
205 *
206 * RCDVIN (input) DOUBLE PRECISION array, dimension (N)
207 * When COMP = .TRUE. RCDVIN holds the precomputed reciprocal
208 * condition numbers for eigenvectors to be compared with
209 * RCONDV.
210 *
211 * RCONDE (workspace) DOUBLE PRECISION array, dimension (N)
212 * RCONDE holds the computed reciprocal condition numbers
213 * for eigenvalues.
214 *
215 * RCNDE1 (workspace) DOUBLE PRECISION array, dimension (N)
216 * RCNDE1 holds more computed reciprocal condition numbers
217 * for eigenvalues.
218 *
219 * RCDEIN (input) DOUBLE PRECISION array, dimension (N)
220 * When COMP = .TRUE. RCDEIN holds the precomputed reciprocal
221 * condition numbers for eigenvalues to be compared with
222 * RCONDE.
223 *
224 * SCALE (workspace) DOUBLE PRECISION array, dimension (N)
225 * Holds information describing balancing of matrix.
226 *
227 * SCALE1 (workspace) DOUBLE PRECISION array, dimension (N)
228 * Holds information describing balancing of matrix.
229 *
230 * RESULT (output) DOUBLE PRECISION array, dimension (11)
231 * The values computed by the 11 tests described above.
232 * The values are currently limited to 1/ulp, to avoid
233 * overflow.
234 *
235 * WORK (workspace) DOUBLE PRECISION array, dimension (LWORK)
236 *
237 * LWORK (input) INTEGER
238 * The number of entries in WORK. This must be at least
239 * 3*N, and 6*N+N**2 if tests 9, 10 or 11 are to be performed.
240 *
241 * IWORK (workspace) INTEGER array, dimension (2*N)
242 *
243 * INFO (output) INTEGER
244 * If 0, successful exit.
245 * If <0, input parameter -INFO had an incorrect value.
246 * If >0, DGEEVX returned an error code, the absolute
247 * value of which is returned.
248 *
249 * =====================================================================
250 *
251 *
252 * .. Parameters ..
253 DOUBLE PRECISION ZERO, ONE, TWO
254 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0 )
255 DOUBLE PRECISION EPSIN
256 PARAMETER ( EPSIN = 5.9605D-8 )
257 * ..
258 * .. Local Scalars ..
259 LOGICAL BALOK, NOBAL
260 CHARACTER SENSE
261 INTEGER I, IHI, IHI1, IINFO, ILO, ILO1, ISENS, ISENSM,
262 $ J, JJ, KMIN
263 DOUBLE PRECISION ABNRM, ABNRM1, EPS, SMLNUM, TNRM, TOL, TOLIN,
264 $ ULP, ULPINV, V, VIMIN, VMAX, VMX, VRMIN, VRMX,
265 $ VTST
266 * ..
267 * .. Local Arrays ..
268 CHARACTER SENS( 2 )
269 DOUBLE PRECISION DUM( 1 ), RES( 2 )
270 * ..
271 * .. External Functions ..
272 LOGICAL LSAME
273 DOUBLE PRECISION DLAMCH, DLAPY2, DNRM2
274 EXTERNAL LSAME, DLAMCH, DLAPY2, DNRM2
275 * ..
276 * .. External Subroutines ..
277 EXTERNAL DGEEVX, DGET22, DLACPY, XERBLA
278 * ..
279 * .. Intrinsic Functions ..
280 INTRINSIC ABS, DBLE, MAX, MIN
281 * ..
282 * .. Data statements ..
283 DATA SENS / 'N', 'V' /
284 * ..
285 * .. Executable Statements ..
286 *
287 * Check for errors
288 *
289 NOBAL = LSAME( BALANC, 'N' )
290 BALOK = NOBAL .OR. LSAME( BALANC, 'P' ) .OR.
291 $ LSAME( BALANC, 'S' ) .OR. LSAME( BALANC, 'B' )
292 INFO = 0
293 IF( .NOT.BALOK ) THEN
294 INFO = -2
295 ELSE IF( THRESH.LT.ZERO ) THEN
296 INFO = -4
297 ELSE IF( NOUNIT.LE.0 ) THEN
298 INFO = -6
299 ELSE IF( N.LT.0 ) THEN
300 INFO = -7
301 ELSE IF( LDA.LT.1 .OR. LDA.LT.N ) THEN
302 INFO = -9
303 ELSE IF( LDVL.LT.1 .OR. LDVL.LT.N ) THEN
304 INFO = -16
305 ELSE IF( LDVR.LT.1 .OR. LDVR.LT.N ) THEN
306 INFO = -18
307 ELSE IF( LDLRE.LT.1 .OR. LDLRE.LT.N ) THEN
308 INFO = -20
309 ELSE IF( LWORK.LT.3*N .OR. ( COMP .AND. LWORK.LT.6*N+N*N ) ) THEN
310 INFO = -31
311 END IF
312 *
313 IF( INFO.NE.0 ) THEN
314 CALL XERBLA( 'DGET23', -INFO )
315 RETURN
316 END IF
317 *
318 * Quick return if nothing to do
319 *
320 DO 10 I = 1, 11
321 RESULT( I ) = -ONE
322 10 CONTINUE
323 *
324 IF( N.EQ.0 )
325 $ RETURN
326 *
327 * More Important constants
328 *
329 ULP = DLAMCH( 'Precision' )
330 SMLNUM = DLAMCH( 'S' )
331 ULPINV = ONE / ULP
332 *
333 * Compute eigenvalues and eigenvectors, and test them
334 *
335 IF( LWORK.GE.6*N+N*N ) THEN
336 SENSE = 'B'
337 ISENSM = 2
338 ELSE
339 SENSE = 'E'
340 ISENSM = 1
341 END IF
342 CALL DLACPY( 'F', N, N, A, LDA, H, LDA )
343 CALL DGEEVX( BALANC, 'V', 'V', SENSE, N, H, LDA, WR, WI, VL, LDVL,
344 $ VR, LDVR, ILO, IHI, SCALE, ABNRM, RCONDE, RCONDV,
345 $ WORK, LWORK, IWORK, IINFO )
346 IF( IINFO.NE.0 ) THEN
347 RESULT( 1 ) = ULPINV
348 IF( JTYPE.NE.22 ) THEN
349 WRITE( NOUNIT, FMT = 9998 )'DGEEVX1', IINFO, N, JTYPE,
350 $ BALANC, ISEED
351 ELSE
352 WRITE( NOUNIT, FMT = 9999 )'DGEEVX1', IINFO, N, ISEED( 1 )
353 END IF
354 INFO = ABS( IINFO )
355 RETURN
356 END IF
357 *
358 * Do Test (1)
359 *
360 CALL DGET22( 'N', 'N', 'N', N, A, LDA, VR, LDVR, WR, WI, WORK,
361 $ RES )
362 RESULT( 1 ) = RES( 1 )
363 *
364 * Do Test (2)
365 *
366 CALL DGET22( 'T', 'N', 'T', N, A, LDA, VL, LDVL, WR, WI, WORK,
367 $ RES )
368 RESULT( 2 ) = RES( 1 )
369 *
370 * Do Test (3)
371 *
372 DO 30 J = 1, N
373 TNRM = ONE
374 IF( WI( J ).EQ.ZERO ) THEN
375 TNRM = DNRM2( N, VR( 1, J ), 1 )
376 ELSE IF( WI( J ).GT.ZERO ) THEN
377 TNRM = DLAPY2( DNRM2( N, VR( 1, J ), 1 ),
378 $ DNRM2( N, VR( 1, J+1 ), 1 ) )
379 END IF
380 RESULT( 3 ) = MAX( RESULT( 3 ),
381 $ MIN( ULPINV, ABS( TNRM-ONE ) / ULP ) )
382 IF( WI( J ).GT.ZERO ) THEN
383 VMX = ZERO
384 VRMX = ZERO
385 DO 20 JJ = 1, N
386 VTST = DLAPY2( VR( JJ, J ), VR( JJ, J+1 ) )
387 IF( VTST.GT.VMX )
388 $ VMX = VTST
389 IF( VR( JJ, J+1 ).EQ.ZERO .AND. ABS( VR( JJ, J ) ).GT.
390 $ VRMX )VRMX = ABS( VR( JJ, J ) )
391 20 CONTINUE
392 IF( VRMX / VMX.LT.ONE-TWO*ULP )
393 $ RESULT( 3 ) = ULPINV
394 END IF
395 30 CONTINUE
396 *
397 * Do Test (4)
398 *
399 DO 50 J = 1, N
400 TNRM = ONE
401 IF( WI( J ).EQ.ZERO ) THEN
402 TNRM = DNRM2( N, VL( 1, J ), 1 )
403 ELSE IF( WI( J ).GT.ZERO ) THEN
404 TNRM = DLAPY2( DNRM2( N, VL( 1, J ), 1 ),
405 $ DNRM2( N, VL( 1, J+1 ), 1 ) )
406 END IF
407 RESULT( 4 ) = MAX( RESULT( 4 ),
408 $ MIN( ULPINV, ABS( TNRM-ONE ) / ULP ) )
409 IF( WI( J ).GT.ZERO ) THEN
410 VMX = ZERO
411 VRMX = ZERO
412 DO 40 JJ = 1, N
413 VTST = DLAPY2( VL( JJ, J ), VL( JJ, J+1 ) )
414 IF( VTST.GT.VMX )
415 $ VMX = VTST
416 IF( VL( JJ, J+1 ).EQ.ZERO .AND. ABS( VL( JJ, J ) ).GT.
417 $ VRMX )VRMX = ABS( VL( JJ, J ) )
418 40 CONTINUE
419 IF( VRMX / VMX.LT.ONE-TWO*ULP )
420 $ RESULT( 4 ) = ULPINV
421 END IF
422 50 CONTINUE
423 *
424 * Test for all options of computing condition numbers
425 *
426 DO 200 ISENS = 1, ISENSM
427 *
428 SENSE = SENS( ISENS )
429 *
430 * Compute eigenvalues only, and test them
431 *
432 CALL DLACPY( 'F', N, N, A, LDA, H, LDA )
433 CALL DGEEVX( BALANC, 'N', 'N', SENSE, N, H, LDA, WR1, WI1, DUM,
434 $ 1, DUM, 1, ILO1, IHI1, SCALE1, ABNRM1, RCNDE1,
435 $ RCNDV1, WORK, LWORK, IWORK, IINFO )
436 IF( IINFO.NE.0 ) THEN
437 RESULT( 1 ) = ULPINV
438 IF( JTYPE.NE.22 ) THEN
439 WRITE( NOUNIT, FMT = 9998 )'DGEEVX2', IINFO, N, JTYPE,
440 $ BALANC, ISEED
441 ELSE
442 WRITE( NOUNIT, FMT = 9999 )'DGEEVX2', IINFO, N,
443 $ ISEED( 1 )
444 END IF
445 INFO = ABS( IINFO )
446 GO TO 190
447 END IF
448 *
449 * Do Test (5)
450 *
451 DO 60 J = 1, N
452 IF( WR( J ).NE.WR1( J ) .OR. WI( J ).NE.WI1( J ) )
453 $ RESULT( 5 ) = ULPINV
454 60 CONTINUE
455 *
456 * Do Test (8)
457 *
458 IF( .NOT.NOBAL ) THEN
459 DO 70 J = 1, N
460 IF( SCALE( J ).NE.SCALE1( J ) )
461 $ RESULT( 8 ) = ULPINV
462 70 CONTINUE
463 IF( ILO.NE.ILO1 )
464 $ RESULT( 8 ) = ULPINV
465 IF( IHI.NE.IHI1 )
466 $ RESULT( 8 ) = ULPINV
467 IF( ABNRM.NE.ABNRM1 )
468 $ RESULT( 8 ) = ULPINV
469 END IF
470 *
471 * Do Test (9)
472 *
473 IF( ISENS.EQ.2 .AND. N.GT.1 ) THEN
474 DO 80 J = 1, N
475 IF( RCONDV( J ).NE.RCNDV1( J ) )
476 $ RESULT( 9 ) = ULPINV
477 80 CONTINUE
478 END IF
479 *
480 * Compute eigenvalues and right eigenvectors, and test them
481 *
482 CALL DLACPY( 'F', N, N, A, LDA, H, LDA )
483 CALL DGEEVX( BALANC, 'N', 'V', SENSE, N, H, LDA, WR1, WI1, DUM,
484 $ 1, LRE, LDLRE, ILO1, IHI1, SCALE1, ABNRM1, RCNDE1,
485 $ RCNDV1, WORK, LWORK, IWORK, IINFO )
486 IF( IINFO.NE.0 ) THEN
487 RESULT( 1 ) = ULPINV
488 IF( JTYPE.NE.22 ) THEN
489 WRITE( NOUNIT, FMT = 9998 )'DGEEVX3', IINFO, N, JTYPE,
490 $ BALANC, ISEED
491 ELSE
492 WRITE( NOUNIT, FMT = 9999 )'DGEEVX3', IINFO, N,
493 $ ISEED( 1 )
494 END IF
495 INFO = ABS( IINFO )
496 GO TO 190
497 END IF
498 *
499 * Do Test (5) again
500 *
501 DO 90 J = 1, N
502 IF( WR( J ).NE.WR1( J ) .OR. WI( J ).NE.WI1( J ) )
503 $ RESULT( 5 ) = ULPINV
504 90 CONTINUE
505 *
506 * Do Test (6)
507 *
508 DO 110 J = 1, N
509 DO 100 JJ = 1, N
510 IF( VR( J, JJ ).NE.LRE( J, JJ ) )
511 $ RESULT( 6 ) = ULPINV
512 100 CONTINUE
513 110 CONTINUE
514 *
515 * Do Test (8) again
516 *
517 IF( .NOT.NOBAL ) THEN
518 DO 120 J = 1, N
519 IF( SCALE( J ).NE.SCALE1( J ) )
520 $ RESULT( 8 ) = ULPINV
521 120 CONTINUE
522 IF( ILO.NE.ILO1 )
523 $ RESULT( 8 ) = ULPINV
524 IF( IHI.NE.IHI1 )
525 $ RESULT( 8 ) = ULPINV
526 IF( ABNRM.NE.ABNRM1 )
527 $ RESULT( 8 ) = ULPINV
528 END IF
529 *
530 * Do Test (9) again
531 *
532 IF( ISENS.EQ.2 .AND. N.GT.1 ) THEN
533 DO 130 J = 1, N
534 IF( RCONDV( J ).NE.RCNDV1( J ) )
535 $ RESULT( 9 ) = ULPINV
536 130 CONTINUE
537 END IF
538 *
539 * Compute eigenvalues and left eigenvectors, and test them
540 *
541 CALL DLACPY( 'F', N, N, A, LDA, H, LDA )
542 CALL DGEEVX( BALANC, 'V', 'N', SENSE, N, H, LDA, WR1, WI1, LRE,
543 $ LDLRE, DUM, 1, ILO1, IHI1, SCALE1, ABNRM1, RCNDE1,
544 $ RCNDV1, WORK, LWORK, IWORK, IINFO )
545 IF( IINFO.NE.0 ) THEN
546 RESULT( 1 ) = ULPINV
547 IF( JTYPE.NE.22 ) THEN
548 WRITE( NOUNIT, FMT = 9998 )'DGEEVX4', IINFO, N, JTYPE,
549 $ BALANC, ISEED
550 ELSE
551 WRITE( NOUNIT, FMT = 9999 )'DGEEVX4', IINFO, N,
552 $ ISEED( 1 )
553 END IF
554 INFO = ABS( IINFO )
555 GO TO 190
556 END IF
557 *
558 * Do Test (5) again
559 *
560 DO 140 J = 1, N
561 IF( WR( J ).NE.WR1( J ) .OR. WI( J ).NE.WI1( J ) )
562 $ RESULT( 5 ) = ULPINV
563 140 CONTINUE
564 *
565 * Do Test (7)
566 *
567 DO 160 J = 1, N
568 DO 150 JJ = 1, N
569 IF( VL( J, JJ ).NE.LRE( J, JJ ) )
570 $ RESULT( 7 ) = ULPINV
571 150 CONTINUE
572 160 CONTINUE
573 *
574 * Do Test (8) again
575 *
576 IF( .NOT.NOBAL ) THEN
577 DO 170 J = 1, N
578 IF( SCALE( J ).NE.SCALE1( J ) )
579 $ RESULT( 8 ) = ULPINV
580 170 CONTINUE
581 IF( ILO.NE.ILO1 )
582 $ RESULT( 8 ) = ULPINV
583 IF( IHI.NE.IHI1 )
584 $ RESULT( 8 ) = ULPINV
585 IF( ABNRM.NE.ABNRM1 )
586 $ RESULT( 8 ) = ULPINV
587 END IF
588 *
589 * Do Test (9) again
590 *
591 IF( ISENS.EQ.2 .AND. N.GT.1 ) THEN
592 DO 180 J = 1, N
593 IF( RCONDV( J ).NE.RCNDV1( J ) )
594 $ RESULT( 9 ) = ULPINV
595 180 CONTINUE
596 END IF
597 *
598 190 CONTINUE
599 *
600 200 CONTINUE
601 *
602 * If COMP, compare condition numbers to precomputed ones
603 *
604 IF( COMP ) THEN
605 CALL DLACPY( 'F', N, N, A, LDA, H, LDA )
606 CALL DGEEVX( 'N', 'V', 'V', 'B', N, H, LDA, WR, WI, VL, LDVL,
607 $ VR, LDVR, ILO, IHI, SCALE, ABNRM, RCONDE, RCONDV,
608 $ WORK, LWORK, IWORK, IINFO )
609 IF( IINFO.NE.0 ) THEN
610 RESULT( 1 ) = ULPINV
611 WRITE( NOUNIT, FMT = 9999 )'DGEEVX5', IINFO, N, ISEED( 1 )
612 INFO = ABS( IINFO )
613 GO TO 250
614 END IF
615 *
616 * Sort eigenvalues and condition numbers lexicographically
617 * to compare with inputs
618 *
619 DO 220 I = 1, N - 1
620 KMIN = I
621 VRMIN = WR( I )
622 VIMIN = WI( I )
623 DO 210 J = I + 1, N
624 IF( WR( J ).LT.VRMIN ) THEN
625 KMIN = J
626 VRMIN = WR( J )
627 VIMIN = WI( J )
628 END IF
629 210 CONTINUE
630 WR( KMIN ) = WR( I )
631 WI( KMIN ) = WI( I )
632 WR( I ) = VRMIN
633 WI( I ) = VIMIN
634 VRMIN = RCONDE( KMIN )
635 RCONDE( KMIN ) = RCONDE( I )
636 RCONDE( I ) = VRMIN
637 VRMIN = RCONDV( KMIN )
638 RCONDV( KMIN ) = RCONDV( I )
639 RCONDV( I ) = VRMIN
640 220 CONTINUE
641 *
642 * Compare condition numbers for eigenvectors
643 * taking their condition numbers into account
644 *
645 RESULT( 10 ) = ZERO
646 EPS = MAX( EPSIN, ULP )
647 V = MAX( DBLE( N )*EPS*ABNRM, SMLNUM )
648 IF( ABNRM.EQ.ZERO )
649 $ V = ONE
650 DO 230 I = 1, N
651 IF( V.GT.RCONDV( I )*RCONDE( I ) ) THEN
652 TOL = RCONDV( I )
653 ELSE
654 TOL = V / RCONDE( I )
655 END IF
656 IF( V.GT.RCDVIN( I )*RCDEIN( I ) ) THEN
657 TOLIN = RCDVIN( I )
658 ELSE
659 TOLIN = V / RCDEIN( I )
660 END IF
661 TOL = MAX( TOL, SMLNUM / EPS )
662 TOLIN = MAX( TOLIN, SMLNUM / EPS )
663 IF( EPS*( RCDVIN( I )-TOLIN ).GT.RCONDV( I )+TOL ) THEN
664 VMAX = ONE / EPS
665 ELSE IF( RCDVIN( I )-TOLIN.GT.RCONDV( I )+TOL ) THEN
666 VMAX = ( RCDVIN( I )-TOLIN ) / ( RCONDV( I )+TOL )
667 ELSE IF( RCDVIN( I )+TOLIN.LT.EPS*( RCONDV( I )-TOL ) ) THEN
668 VMAX = ONE / EPS
669 ELSE IF( RCDVIN( I )+TOLIN.LT.RCONDV( I )-TOL ) THEN
670 VMAX = ( RCONDV( I )-TOL ) / ( RCDVIN( I )+TOLIN )
671 ELSE
672 VMAX = ONE
673 END IF
674 RESULT( 10 ) = MAX( RESULT( 10 ), VMAX )
675 230 CONTINUE
676 *
677 * Compare condition numbers for eigenvalues
678 * taking their condition numbers into account
679 *
680 RESULT( 11 ) = ZERO
681 DO 240 I = 1, N
682 IF( V.GT.RCONDV( I ) ) THEN
683 TOL = ONE
684 ELSE
685 TOL = V / RCONDV( I )
686 END IF
687 IF( V.GT.RCDVIN( I ) ) THEN
688 TOLIN = ONE
689 ELSE
690 TOLIN = V / RCDVIN( I )
691 END IF
692 TOL = MAX( TOL, SMLNUM / EPS )
693 TOLIN = MAX( TOLIN, SMLNUM / EPS )
694 IF( EPS*( RCDEIN( I )-TOLIN ).GT.RCONDE( I )+TOL ) THEN
695 VMAX = ONE / EPS
696 ELSE IF( RCDEIN( I )-TOLIN.GT.RCONDE( I )+TOL ) THEN
697 VMAX = ( RCDEIN( I )-TOLIN ) / ( RCONDE( I )+TOL )
698 ELSE IF( RCDEIN( I )+TOLIN.LT.EPS*( RCONDE( I )-TOL ) ) THEN
699 VMAX = ONE / EPS
700 ELSE IF( RCDEIN( I )+TOLIN.LT.RCONDE( I )-TOL ) THEN
701 VMAX = ( RCONDE( I )-TOL ) / ( RCDEIN( I )+TOLIN )
702 ELSE
703 VMAX = ONE
704 END IF
705 RESULT( 11 ) = MAX( RESULT( 11 ), VMAX )
706 240 CONTINUE
707 250 CONTINUE
708 *
709 END IF
710 *
711 9999 FORMAT( ' DGET23: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
712 $ I6, ', INPUT EXAMPLE NUMBER = ', I4 )
713 9998 FORMAT( ' DGET23: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
714 $ I6, ', JTYPE=', I6, ', BALANC = ', A, ', ISEED=(',
715 $ 3( I5, ',' ), I5, ')' )
716 *
717 RETURN
718 *
719 * End of DGET23
720 *
721 END