1 SUBROUTINE CDRGVX( NSIZE, THRESH, NIN, NOUT, A, LDA, B, AI, BI,
2 $ ALPHA, BETA, VL, VR, ILO, IHI, LSCALE, RSCALE,
3 $ S, STRU, DIF, DIFTRU, WORK, LWORK, RWORK,
4 $ IWORK, LIWORK, RESULT, BWORK, INFO )
5 *
6 * -- LAPACK test routine (version 3.1) --
7 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
8 * November 2006
9 *
10 * .. Scalar Arguments ..
11 INTEGER IHI, ILO, INFO, LDA, LIWORK, LWORK, NIN, NOUT,
12 $ NSIZE
13 REAL THRESH
14 * ..
15 * .. Array Arguments ..
16 LOGICAL BWORK( * )
17 INTEGER IWORK( * )
18 REAL DIF( * ), DIFTRU( * ), LSCALE( * ),
19 $ RESULT( 4 ), RSCALE( * ), RWORK( * ), S( * ),
20 $ STRU( * )
21 COMPLEX A( LDA, * ), AI( LDA, * ), ALPHA( * ),
22 $ B( LDA, * ), BETA( * ), BI( LDA, * ),
23 $ VL( LDA, * ), VR( LDA, * ), WORK( * )
24 * ..
25 *
26 * Purpose
27 * =======
28 *
29 * CDRGVX checks the nonsymmetric generalized eigenvalue problem
30 * expert driver CGGEVX.
31 *
32 * CGGEVX computes the generalized eigenvalues, (optionally) the left
33 * and/or right eigenvectors, (optionally) computes a balancing
34 * transformation to improve the conditioning, and (optionally)
35 * reciprocal condition numbers for the eigenvalues and eigenvectors.
36 *
37 * When CDRGVX is called with NSIZE > 0, two types of test matrix pairs
38 * are generated by the subroutine SLATM6 and test the driver CGGEVX.
39 * The test matrices have the known exact condition numbers for
40 * eigenvalues. For the condition numbers of the eigenvectors
41 * corresponding the first and last eigenvalues are also know
42 * ``exactly'' (see CLATM6).
43 * For each matrix pair, the following tests will be performed and
44 * compared with the threshhold THRESH.
45 *
46 * (1) max over all left eigenvalue/-vector pairs (beta/alpha,l) of
47 *
48 * | l**H * (beta A - alpha B) | / ( ulp max( |beta A|, |alpha B| ) )
49 *
50 * where l**H is the conjugate tranpose of l.
51 *
52 * (2) max over all right eigenvalue/-vector pairs (beta/alpha,r) of
53 *
54 * | (beta A - alpha B) r | / ( ulp max( |beta A|, |alpha B| ) )
55 *
56 * (3) The condition number S(i) of eigenvalues computed by CGGEVX
57 * differs less than a factor THRESH from the exact S(i) (see
58 * CLATM6).
59 *
60 * (4) DIF(i) computed by CTGSNA differs less than a factor 10*THRESH
61 * from the exact value (for the 1st and 5th vectors only).
62 *
63 * Test Matrices
64 * =============
65 *
66 * Two kinds of test matrix pairs
67 * (A, B) = inverse(YH) * (Da, Db) * inverse(X)
68 * are used in the tests:
69 *
70 * 1: Da = 1+a 0 0 0 0 Db = 1 0 0 0 0
71 * 0 2+a 0 0 0 0 1 0 0 0
72 * 0 0 3+a 0 0 0 0 1 0 0
73 * 0 0 0 4+a 0 0 0 0 1 0
74 * 0 0 0 0 5+a , 0 0 0 0 1 , and
75 *
76 * 2: Da = 1 -1 0 0 0 Db = 1 0 0 0 0
77 * 1 1 0 0 0 0 1 0 0 0
78 * 0 0 1 0 0 0 0 1 0 0
79 * 0 0 0 1+a 1+b 0 0 0 1 0
80 * 0 0 0 -1-b 1+a , 0 0 0 0 1 .
81 *
82 * In both cases the same inverse(YH) and inverse(X) are used to compute
83 * (A, B), giving the exact eigenvectors to (A,B) as (YH, X):
84 *
85 * YH: = 1 0 -y y -y X = 1 0 -x -x x
86 * 0 1 -y y -y 0 1 x -x -x
87 * 0 0 1 0 0 0 0 1 0 0
88 * 0 0 0 1 0 0 0 0 1 0
89 * 0 0 0 0 1, 0 0 0 0 1 , where
90 *
91 * a, b, x and y will have all values independently of each other from
92 * { sqrt(sqrt(ULP)), 0.1, 1, 10, 1/sqrt(sqrt(ULP)) }.
93 *
94 * Arguments
95 * =========
96 *
97 * NSIZE (input) INTEGER
98 * The number of sizes of matrices to use. NSIZE must be at
99 * least zero. If it is zero, no randomly generated matrices
100 * are tested, but any test matrices read from NIN will be
101 * tested. If it is not zero, then N = 5.
102 *
103 * THRESH (input) REAL
104 * A test will count as "failed" if the "error", computed as
105 * described above, exceeds THRESH. Note that the error
106 * is scaled to be O(1), so THRESH should be a reasonably
107 * small multiple of 1, e.g., 10 or 100. In particular,
108 * it should not depend on the precision (single vs. double)
109 * or the size of the matrix. It must be at least zero.
110 *
111 * NIN (input) INTEGER
112 * The FORTRAN unit number for reading in the data file of
113 * problems to solve.
114 *
115 * NOUT (input) INTEGER
116 * The FORTRAN unit number for printing out error messages
117 * (e.g., if a routine returns IINFO not equal to 0.)
118 *
119 * A (workspace) COMPLEX array, dimension (LDA, NSIZE)
120 * Used to hold the matrix whose eigenvalues are to be
121 * computed. On exit, A contains the last matrix actually used.
122 *
123 * LDA (input) INTEGER
124 * The leading dimension of A, B, AI, BI, Ao, and Bo.
125 * It must be at least 1 and at least NSIZE.
126 *
127 * B (workspace) COMPLEX array, dimension (LDA, NSIZE)
128 * Used to hold the matrix whose eigenvalues are to be
129 * computed. On exit, B contains the last matrix actually used.
130 *
131 * AI (workspace) COMPLEX array, dimension (LDA, NSIZE)
132 * Copy of A, modified by CGGEVX.
133 *
134 * BI (workspace) COMPLEX array, dimension (LDA, NSIZE)
135 * Copy of B, modified by CGGEVX.
136 *
137 * ALPHA (workspace) COMPLEX array, dimension (NSIZE)
138 * BETA (workspace) COMPLEX array, dimension (NSIZE)
139 * On exit, ALPHA/BETA are the eigenvalues.
140 *
141 * VL (workspace) COMPLEX array, dimension (LDA, NSIZE)
142 * VL holds the left eigenvectors computed by CGGEVX.
143 *
144 * VR (workspace) COMPLEX array, dimension (LDA, NSIZE)
145 * VR holds the right eigenvectors computed by CGGEVX.
146 *
147 * ILO (output/workspace) INTEGER
148 *
149 * IHI (output/workspace) INTEGER
150 *
151 * LSCALE (output/workspace) REAL array, dimension (N)
152 *
153 * RSCALE (output/workspace) REAL array, dimension (N)
154 *
155 * S (output/workspace) REAL array, dimension (N)
156 *
157 * STRU (output/workspace) REAL array, dimension (N)
158 *
159 * DIF (output/workspace) REAL array, dimension (N)
160 *
161 * DIFTRU (output/workspace) REAL array, dimension (N)
162 *
163 * WORK (workspace) COMPLEX array, dimension (LWORK)
164 *
165 * LWORK (input) INTEGER
166 * Leading dimension of WORK. LWORK >= 2*N*N + 2*N
167 *
168 * RWORK (workspace) REAL array, dimension (6*N)
169 *
170 * IWORK (workspace) INTEGER array, dimension (LIWORK)
171 *
172 * LIWORK (input) INTEGER
173 * Leading dimension of IWORK. LIWORK >= N+2.
174 *
175 * RESULT (output/workspace) REAL array, dimension (4)
176 *
177 * BWORK (workspace) LOGICAL array, dimension (N)
178 *
179 * INFO (output) INTEGER
180 * = 0: successful exit
181 * < 0: if INFO = -i, the i-th argument had an illegal value.
182 * > 0: A routine returned an error code.
183 *
184 * =====================================================================
185 *
186 * .. Parameters ..
187 REAL ZERO, ONE, TEN, TNTH
188 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0, TEN = 1.0E+1,
189 $ TNTH = 1.0E-1 )
190 * ..
191 * .. Local Scalars ..
192 INTEGER I, IPTYPE, IWA, IWB, IWX, IWY, J, LINFO,
193 $ MAXWRK, MINWRK, N, NERRS, NMAX, NPTKNT, NTESTT
194 REAL ABNORM, ANORM, BNORM, RATIO1, RATIO2, THRSH2,
195 $ ULP, ULPINV
196 * ..
197 * .. Local Arrays ..
198 COMPLEX WEIGHT( 5 )
199 * ..
200 * .. External Functions ..
201 INTEGER ILAENV
202 REAL CLANGE, SLAMCH
203 EXTERNAL ILAENV, CLANGE, SLAMCH
204 * ..
205 * .. External Subroutines ..
206 EXTERNAL ALASVM, CGET52, CGGEVX, CLACPY, CLATM6, XERBLA
207 * ..
208 * .. Intrinsic Functions ..
209 INTRINSIC ABS, CMPLX, MAX, SQRT
210 * ..
211 * .. Executable Statements ..
212 *
213 * Check for errors
214 *
215 INFO = 0
216 *
217 NMAX = 5
218 *
219 IF( NSIZE.LT.0 ) THEN
220 INFO = -1
221 ELSE IF( THRESH.LT.ZERO ) THEN
222 INFO = -2
223 ELSE IF( NIN.LE.0 ) THEN
224 INFO = -3
225 ELSE IF( NOUT.LE.0 ) THEN
226 INFO = -4
227 ELSE IF( LDA.LT.1 .OR. LDA.LT.NMAX ) THEN
228 INFO = -6
229 ELSE IF( LIWORK.LT.NMAX+2 ) THEN
230 INFO = -26
231 END IF
232 *
233 * Compute workspace
234 * (Note: Comments in the code beginning "Workspace:" describe the
235 * minimal amount of workspace needed at that point in the code,
236 * as well as the preferred amount for good performance.
237 * NB refers to the optimal block size for the immediately
238 * following subroutine, as returned by ILAENV.)
239 *
240 MINWRK = 1
241 IF( INFO.EQ.0 .AND. LWORK.GE.1 ) THEN
242 MINWRK = 2*NMAX*( NMAX+1 )
243 MAXWRK = NMAX*( 1+ILAENV( 1, 'CGEQRF', ' ', NMAX, 1, NMAX,
244 $ 0 ) )
245 MAXWRK = MAX( MAXWRK, 2*NMAX*( NMAX+1 ) )
246 WORK( 1 ) = MAXWRK
247 END IF
248 *
249 IF( LWORK.LT.MINWRK )
250 $ INFO = -23
251 *
252 IF( INFO.NE.0 ) THEN
253 CALL XERBLA( 'CDRGVX', -INFO )
254 RETURN
255 END IF
256 *
257 N = 5
258 ULP = SLAMCH( 'P' )
259 ULPINV = ONE / ULP
260 THRSH2 = TEN*THRESH
261 NERRS = 0
262 NPTKNT = 0
263 NTESTT = 0
264 *
265 IF( NSIZE.EQ.0 )
266 $ GO TO 90
267 *
268 * Parameters used for generating test matrices.
269 *
270 WEIGHT( 1 ) = CMPLX( SQRT( SQRT( ULP ) ), ZERO )
271 WEIGHT( 2 ) = CMPLX( TNTH, ZERO )
272 WEIGHT( 3 ) = ONE
273 WEIGHT( 4 ) = ONE / WEIGHT( 2 )
274 WEIGHT( 5 ) = ONE / WEIGHT( 1 )
275 *
276 DO 80 IPTYPE = 1, 2
277 DO 70 IWA = 1, 5
278 DO 60 IWB = 1, 5
279 DO 50 IWX = 1, 5
280 DO 40 IWY = 1, 5
281 *
282 * generated a pair of test matrix
283 *
284 CALL CLATM6( IPTYPE, 5, A, LDA, B, VR, LDA, VL,
285 $ LDA, WEIGHT( IWA ), WEIGHT( IWB ),
286 $ WEIGHT( IWX ), WEIGHT( IWY ), STRU,
287 $ DIFTRU )
288 *
289 * Compute eigenvalues/eigenvectors of (A, B).
290 * Compute eigenvalue/eigenvector condition numbers
291 * using computed eigenvectors.
292 *
293 CALL CLACPY( 'F', N, N, A, LDA, AI, LDA )
294 CALL CLACPY( 'F', N, N, B, LDA, BI, LDA )
295 *
296 CALL CGGEVX( 'N', 'V', 'V', 'B', N, AI, LDA, BI,
297 $ LDA, ALPHA, BETA, VL, LDA, VR, LDA,
298 $ ILO, IHI, LSCALE, RSCALE, ANORM,
299 $ BNORM, S, DIF, WORK, LWORK, RWORK,
300 $ IWORK, BWORK, LINFO )
301 IF( LINFO.NE.0 ) THEN
302 WRITE( NOUT, FMT = 9999 )'CGGEVX', LINFO, N,
303 $ IPTYPE, IWA, IWB, IWX, IWY
304 GO TO 30
305 END IF
306 *
307 * Compute the norm(A, B)
308 *
309 CALL CLACPY( 'Full', N, N, AI, LDA, WORK, N )
310 CALL CLACPY( 'Full', N, N, BI, LDA, WORK( N*N+1 ),
311 $ N )
312 ABNORM = CLANGE( 'Fro', N, 2*N, WORK, N, RWORK )
313 *
314 * Tests (1) and (2)
315 *
316 RESULT( 1 ) = ZERO
317 CALL CGET52( .TRUE., N, A, LDA, B, LDA, VL, LDA,
318 $ ALPHA, BETA, WORK, RWORK,
319 $ RESULT( 1 ) )
320 IF( RESULT( 2 ).GT.THRESH ) THEN
321 WRITE( NOUT, FMT = 9998 )'Left', 'CGGEVX',
322 $ RESULT( 2 ), N, IPTYPE, IWA, IWB, IWX, IWY
323 END IF
324 *
325 RESULT( 2 ) = ZERO
326 CALL CGET52( .FALSE., N, A, LDA, B, LDA, VR, LDA,
327 $ ALPHA, BETA, WORK, RWORK,
328 $ RESULT( 2 ) )
329 IF( RESULT( 3 ).GT.THRESH ) THEN
330 WRITE( NOUT, FMT = 9998 )'Right', 'CGGEVX',
331 $ RESULT( 3 ), N, IPTYPE, IWA, IWB, IWX, IWY
332 END IF
333 *
334 * Test (3)
335 *
336 RESULT( 3 ) = ZERO
337 DO 10 I = 1, N
338 IF( S( I ).EQ.ZERO ) THEN
339 IF( STRU( I ).GT.ABNORM*ULP )
340 $ RESULT( 3 ) = ULPINV
341 ELSE IF( STRU( I ).EQ.ZERO ) THEN
342 IF( S( I ).GT.ABNORM*ULP )
343 $ RESULT( 3 ) = ULPINV
344 ELSE
345 RWORK( I ) = MAX( ABS( STRU( I ) / S( I ) ),
346 $ ABS( S( I ) / STRU( I ) ) )
347 RESULT( 3 ) = MAX( RESULT( 3 ), RWORK( I ) )
348 END IF
349 10 CONTINUE
350 *
351 * Test (4)
352 *
353 RESULT( 4 ) = ZERO
354 IF( DIF( 1 ).EQ.ZERO ) THEN
355 IF( DIFTRU( 1 ).GT.ABNORM*ULP )
356 $ RESULT( 4 ) = ULPINV
357 ELSE IF( DIFTRU( 1 ).EQ.ZERO ) THEN
358 IF( DIF( 1 ).GT.ABNORM*ULP )
359 $ RESULT( 4 ) = ULPINV
360 ELSE IF( DIF( 5 ).EQ.ZERO ) THEN
361 IF( DIFTRU( 5 ).GT.ABNORM*ULP )
362 $ RESULT( 4 ) = ULPINV
363 ELSE IF( DIFTRU( 5 ).EQ.ZERO ) THEN
364 IF( DIF( 5 ).GT.ABNORM*ULP )
365 $ RESULT( 4 ) = ULPINV
366 ELSE
367 RATIO1 = MAX( ABS( DIFTRU( 1 ) / DIF( 1 ) ),
368 $ ABS( DIF( 1 ) / DIFTRU( 1 ) ) )
369 RATIO2 = MAX( ABS( DIFTRU( 5 ) / DIF( 5 ) ),
370 $ ABS( DIF( 5 ) / DIFTRU( 5 ) ) )
371 RESULT( 4 ) = MAX( RATIO1, RATIO2 )
372 END IF
373 *
374 NTESTT = NTESTT + 4
375 *
376 * Print out tests which fail.
377 *
378 DO 20 J = 1, 4
379 IF( ( RESULT( J ).GE.THRSH2 .AND. J.GE.4 ) .OR.
380 $ ( RESULT( J ).GE.THRESH .AND. J.LE.3 ) )
381 $ THEN
382 *
383 * If this is the first test to fail,
384 * print a header to the data file.
385 *
386 IF( NERRS.EQ.0 ) THEN
387 WRITE( NOUT, FMT = 9997 )'CXV'
388 *
389 * Print out messages for built-in examples
390 *
391 * Matrix types
392 *
393 WRITE( NOUT, FMT = 9995 )
394 WRITE( NOUT, FMT = 9994 )
395 WRITE( NOUT, FMT = 9993 )
396 *
397 * Tests performed
398 *
399 WRITE( NOUT, FMT = 9992 )'''',
400 $ 'transpose', ''''
401 *
402 END IF
403 NERRS = NERRS + 1
404 IF( RESULT( J ).LT.10000.0 ) THEN
405 WRITE( NOUT, FMT = 9991 )IPTYPE, IWA,
406 $ IWB, IWX, IWY, J, RESULT( J )
407 ELSE
408 WRITE( NOUT, FMT = 9990 )IPTYPE, IWA,
409 $ IWB, IWX, IWY, J, RESULT( J )
410 END IF
411 END IF
412 20 CONTINUE
413 *
414 30 CONTINUE
415 *
416 40 CONTINUE
417 50 CONTINUE
418 60 CONTINUE
419 70 CONTINUE
420 80 CONTINUE
421 *
422 GO TO 150
423 *
424 90 CONTINUE
425 *
426 * Read in data from file to check accuracy of condition estimation
427 * Read input data until N=0
428 *
429 READ( NIN, FMT = *, END = 150 )N
430 IF( N.EQ.0 )
431 $ GO TO 150
432 DO 100 I = 1, N
433 READ( NIN, FMT = * )( A( I, J ), J = 1, N )
434 100 CONTINUE
435 DO 110 I = 1, N
436 READ( NIN, FMT = * )( B( I, J ), J = 1, N )
437 110 CONTINUE
438 READ( NIN, FMT = * )( STRU( I ), I = 1, N )
439 READ( NIN, FMT = * )( DIFTRU( I ), I = 1, N )
440 *
441 NPTKNT = NPTKNT + 1
442 *
443 * Compute eigenvalues/eigenvectors of (A, B).
444 * Compute eigenvalue/eigenvector condition numbers
445 * using computed eigenvectors.
446 *
447 CALL CLACPY( 'F', N, N, A, LDA, AI, LDA )
448 CALL CLACPY( 'F', N, N, B, LDA, BI, LDA )
449 *
450 CALL CGGEVX( 'N', 'V', 'V', 'B', N, AI, LDA, BI, LDA, ALPHA, BETA,
451 $ VL, LDA, VR, LDA, ILO, IHI, LSCALE, RSCALE, ANORM,
452 $ BNORM, S, DIF, WORK, LWORK, RWORK, IWORK, BWORK,
453 $ LINFO )
454 *
455 IF( LINFO.NE.0 ) THEN
456 WRITE( NOUT, FMT = 9987 )'CGGEVX', LINFO, N, NPTKNT
457 GO TO 140
458 END IF
459 *
460 * Compute the norm(A, B)
461 *
462 CALL CLACPY( 'Full', N, N, AI, LDA, WORK, N )
463 CALL CLACPY( 'Full', N, N, BI, LDA, WORK( N*N+1 ), N )
464 ABNORM = CLANGE( 'Fro', N, 2*N, WORK, N, RWORK )
465 *
466 * Tests (1) and (2)
467 *
468 RESULT( 1 ) = ZERO
469 CALL CGET52( .TRUE., N, A, LDA, B, LDA, VL, LDA, ALPHA, BETA,
470 $ WORK, RWORK, RESULT( 1 ) )
471 IF( RESULT( 2 ).GT.THRESH ) THEN
472 WRITE( NOUT, FMT = 9986 )'Left', 'CGGEVX', RESULT( 2 ), N,
473 $ NPTKNT
474 END IF
475 *
476 RESULT( 2 ) = ZERO
477 CALL CGET52( .FALSE., N, A, LDA, B, LDA, VR, LDA, ALPHA, BETA,
478 $ WORK, RWORK, RESULT( 2 ) )
479 IF( RESULT( 3 ).GT.THRESH ) THEN
480 WRITE( NOUT, FMT = 9986 )'Right', 'CGGEVX', RESULT( 3 ), N,
481 $ NPTKNT
482 END IF
483 *
484 * Test (3)
485 *
486 RESULT( 3 ) = ZERO
487 DO 120 I = 1, N
488 IF( S( I ).EQ.ZERO ) THEN
489 IF( STRU( I ).GT.ABNORM*ULP )
490 $ RESULT( 3 ) = ULPINV
491 ELSE IF( STRU( I ).EQ.ZERO ) THEN
492 IF( S( I ).GT.ABNORM*ULP )
493 $ RESULT( 3 ) = ULPINV
494 ELSE
495 RWORK( I ) = MAX( ABS( STRU( I ) / S( I ) ),
496 $ ABS( S( I ) / STRU( I ) ) )
497 RESULT( 3 ) = MAX( RESULT( 3 ), RWORK( I ) )
498 END IF
499 120 CONTINUE
500 *
501 * Test (4)
502 *
503 RESULT( 4 ) = ZERO
504 IF( DIF( 1 ).EQ.ZERO ) THEN
505 IF( DIFTRU( 1 ).GT.ABNORM*ULP )
506 $ RESULT( 4 ) = ULPINV
507 ELSE IF( DIFTRU( 1 ).EQ.ZERO ) THEN
508 IF( DIF( 1 ).GT.ABNORM*ULP )
509 $ RESULT( 4 ) = ULPINV
510 ELSE IF( DIF( 5 ).EQ.ZERO ) THEN
511 IF( DIFTRU( 5 ).GT.ABNORM*ULP )
512 $ RESULT( 4 ) = ULPINV
513 ELSE IF( DIFTRU( 5 ).EQ.ZERO ) THEN
514 IF( DIF( 5 ).GT.ABNORM*ULP )
515 $ RESULT( 4 ) = ULPINV
516 ELSE
517 RATIO1 = MAX( ABS( DIFTRU( 1 ) / DIF( 1 ) ),
518 $ ABS( DIF( 1 ) / DIFTRU( 1 ) ) )
519 RATIO2 = MAX( ABS( DIFTRU( 5 ) / DIF( 5 ) ),
520 $ ABS( DIF( 5 ) / DIFTRU( 5 ) ) )
521 RESULT( 4 ) = MAX( RATIO1, RATIO2 )
522 END IF
523 *
524 NTESTT = NTESTT + 4
525 *
526 * Print out tests which fail.
527 *
528 DO 130 J = 1, 4
529 IF( RESULT( J ).GE.THRSH2 ) THEN
530 *
531 * If this is the first test to fail,
532 * print a header to the data file.
533 *
534 IF( NERRS.EQ.0 ) THEN
535 WRITE( NOUT, FMT = 9997 )'CXV'
536 *
537 * Print out messages for built-in examples
538 *
539 * Matrix types
540 *
541 WRITE( NOUT, FMT = 9996 )
542 *
543 * Tests performed
544 *
545 WRITE( NOUT, FMT = 9992 )'''', 'transpose', ''''
546 *
547 END IF
548 NERRS = NERRS + 1
549 IF( RESULT( J ).LT.10000.0 ) THEN
550 WRITE( NOUT, FMT = 9989 )NPTKNT, N, J, RESULT( J )
551 ELSE
552 WRITE( NOUT, FMT = 9988 )NPTKNT, N, J, RESULT( J )
553 END IF
554 END IF
555 130 CONTINUE
556 *
557 140 CONTINUE
558 *
559 GO TO 90
560 150 CONTINUE
561 *
562 * Summary
563 *
564 CALL ALASVM( 'CXV', NOUT, NERRS, NTESTT, 0 )
565 *
566 WORK( 1 ) = MAXWRK
567 *
568 RETURN
569 *
570 9999 FORMAT( ' CDRGVX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
571 $ I6, ', JTYPE=', I6, ')' )
572 *
573 9998 FORMAT( ' CDRGVX: ', A, ' Eigenvectors from ', A, ' incorrectly ',
574 $ 'normalized.', / ' Bits of error=', 0P, G10.3, ',', 9X,
575 $ 'N=', I6, ', JTYPE=', I6, ', IWA=', I5, ', IWB=', I5,
576 $ ', IWX=', I5, ', IWY=', I5 )
577 *
578 9997 FORMAT( / 1X, A3, ' -- Complex Expert Eigenvalue/vector',
579 $ ' problem driver' )
580 *
581 9996 FORMAT( 'Input Example' )
582 *
583 9995 FORMAT( ' Matrix types: ', / )
584 *
585 9994 FORMAT( ' TYPE 1: Da is diagonal, Db is identity, ',
586 $ / ' A = Y^(-H) Da X^(-1), B = Y^(-H) Db X^(-1) ',
587 $ / ' YH and X are left and right eigenvectors. ', / )
588 *
589 9993 FORMAT( ' TYPE 2: Da is quasi-diagonal, Db is identity, ',
590 $ / ' A = Y^(-H) Da X^(-1), B = Y^(-H) Db X^(-1) ',
591 $ / ' YH and X are left and right eigenvectors. ', / )
592 *
593 9992 FORMAT( / ' Tests performed: ', / 4X,
594 $ ' a is alpha, b is beta, l is a left eigenvector, ', / 4X,
595 $ ' r is a right eigenvector and ', A, ' means ', A, '.',
596 $ / ' 1 = max | ( b A - a B )', A, ' l | / const.',
597 $ / ' 2 = max | ( b A - a B ) r | / const.',
598 $ / ' 3 = max ( Sest/Stru, Stru/Sest ) ',
599 $ ' over all eigenvalues', /
600 $ ' 4 = max( DIFest/DIFtru, DIFtru/DIFest ) ',
601 $ ' over the 1st and 5th eigenvectors', / )
602 *
603 9991 FORMAT( ' Type=', I2, ',', ' IWA=', I2, ', IWB=', I2, ', IWX=',
604 $ I2, ', IWY=', I2, ', result ', I2, ' is', 0P, F8.2 )
605 *
606 9990 FORMAT( ' Type=', I2, ',', ' IWA=', I2, ', IWB=', I2, ', IWX=',
607 $ I2, ', IWY=', I2, ', result ', I2, ' is', 1P, E10.3 )
608 *
609 9989 FORMAT( ' Input example #', I2, ', matrix order=', I4, ',',
610 $ ' result ', I2, ' is', 0P, F8.2 )
611 *
612 9988 FORMAT( ' Input example #', I2, ', matrix order=', I4, ',',
613 $ ' result ', I2, ' is', 1P, E10.3 )
614 *
615 9987 FORMAT( ' CDRGVX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
616 $ I6, ', Input example #', I2, ')' )
617 *
618 9986 FORMAT( ' CDRGVX: ', A, ' Eigenvectors from ', A, ' incorrectly ',
619 $ 'normalized.', / ' Bits of error=', 0P, G10.3, ',', 9X,
620 $ 'N=', I6, ', Input Example #', I2, ')' )
621 *
622 * End of CDRGVX
623 *
624 END
2 $ ALPHA, BETA, VL, VR, ILO, IHI, LSCALE, RSCALE,
3 $ S, STRU, DIF, DIFTRU, WORK, LWORK, RWORK,
4 $ IWORK, LIWORK, RESULT, BWORK, INFO )
5 *
6 * -- LAPACK test routine (version 3.1) --
7 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
8 * November 2006
9 *
10 * .. Scalar Arguments ..
11 INTEGER IHI, ILO, INFO, LDA, LIWORK, LWORK, NIN, NOUT,
12 $ NSIZE
13 REAL THRESH
14 * ..
15 * .. Array Arguments ..
16 LOGICAL BWORK( * )
17 INTEGER IWORK( * )
18 REAL DIF( * ), DIFTRU( * ), LSCALE( * ),
19 $ RESULT( 4 ), RSCALE( * ), RWORK( * ), S( * ),
20 $ STRU( * )
21 COMPLEX A( LDA, * ), AI( LDA, * ), ALPHA( * ),
22 $ B( LDA, * ), BETA( * ), BI( LDA, * ),
23 $ VL( LDA, * ), VR( LDA, * ), WORK( * )
24 * ..
25 *
26 * Purpose
27 * =======
28 *
29 * CDRGVX checks the nonsymmetric generalized eigenvalue problem
30 * expert driver CGGEVX.
31 *
32 * CGGEVX computes the generalized eigenvalues, (optionally) the left
33 * and/or right eigenvectors, (optionally) computes a balancing
34 * transformation to improve the conditioning, and (optionally)
35 * reciprocal condition numbers for the eigenvalues and eigenvectors.
36 *
37 * When CDRGVX is called with NSIZE > 0, two types of test matrix pairs
38 * are generated by the subroutine SLATM6 and test the driver CGGEVX.
39 * The test matrices have the known exact condition numbers for
40 * eigenvalues. For the condition numbers of the eigenvectors
41 * corresponding the first and last eigenvalues are also know
42 * ``exactly'' (see CLATM6).
43 * For each matrix pair, the following tests will be performed and
44 * compared with the threshhold THRESH.
45 *
46 * (1) max over all left eigenvalue/-vector pairs (beta/alpha,l) of
47 *
48 * | l**H * (beta A - alpha B) | / ( ulp max( |beta A|, |alpha B| ) )
49 *
50 * where l**H is the conjugate tranpose of l.
51 *
52 * (2) max over all right eigenvalue/-vector pairs (beta/alpha,r) of
53 *
54 * | (beta A - alpha B) r | / ( ulp max( |beta A|, |alpha B| ) )
55 *
56 * (3) The condition number S(i) of eigenvalues computed by CGGEVX
57 * differs less than a factor THRESH from the exact S(i) (see
58 * CLATM6).
59 *
60 * (4) DIF(i) computed by CTGSNA differs less than a factor 10*THRESH
61 * from the exact value (for the 1st and 5th vectors only).
62 *
63 * Test Matrices
64 * =============
65 *
66 * Two kinds of test matrix pairs
67 * (A, B) = inverse(YH) * (Da, Db) * inverse(X)
68 * are used in the tests:
69 *
70 * 1: Da = 1+a 0 0 0 0 Db = 1 0 0 0 0
71 * 0 2+a 0 0 0 0 1 0 0 0
72 * 0 0 3+a 0 0 0 0 1 0 0
73 * 0 0 0 4+a 0 0 0 0 1 0
74 * 0 0 0 0 5+a , 0 0 0 0 1 , and
75 *
76 * 2: Da = 1 -1 0 0 0 Db = 1 0 0 0 0
77 * 1 1 0 0 0 0 1 0 0 0
78 * 0 0 1 0 0 0 0 1 0 0
79 * 0 0 0 1+a 1+b 0 0 0 1 0
80 * 0 0 0 -1-b 1+a , 0 0 0 0 1 .
81 *
82 * In both cases the same inverse(YH) and inverse(X) are used to compute
83 * (A, B), giving the exact eigenvectors to (A,B) as (YH, X):
84 *
85 * YH: = 1 0 -y y -y X = 1 0 -x -x x
86 * 0 1 -y y -y 0 1 x -x -x
87 * 0 0 1 0 0 0 0 1 0 0
88 * 0 0 0 1 0 0 0 0 1 0
89 * 0 0 0 0 1, 0 0 0 0 1 , where
90 *
91 * a, b, x and y will have all values independently of each other from
92 * { sqrt(sqrt(ULP)), 0.1, 1, 10, 1/sqrt(sqrt(ULP)) }.
93 *
94 * Arguments
95 * =========
96 *
97 * NSIZE (input) INTEGER
98 * The number of sizes of matrices to use. NSIZE must be at
99 * least zero. If it is zero, no randomly generated matrices
100 * are tested, but any test matrices read from NIN will be
101 * tested. If it is not zero, then N = 5.
102 *
103 * THRESH (input) REAL
104 * A test will count as "failed" if the "error", computed as
105 * described above, exceeds THRESH. Note that the error
106 * is scaled to be O(1), so THRESH should be a reasonably
107 * small multiple of 1, e.g., 10 or 100. In particular,
108 * it should not depend on the precision (single vs. double)
109 * or the size of the matrix. It must be at least zero.
110 *
111 * NIN (input) INTEGER
112 * The FORTRAN unit number for reading in the data file of
113 * problems to solve.
114 *
115 * NOUT (input) INTEGER
116 * The FORTRAN unit number for printing out error messages
117 * (e.g., if a routine returns IINFO not equal to 0.)
118 *
119 * A (workspace) COMPLEX array, dimension (LDA, NSIZE)
120 * Used to hold the matrix whose eigenvalues are to be
121 * computed. On exit, A contains the last matrix actually used.
122 *
123 * LDA (input) INTEGER
124 * The leading dimension of A, B, AI, BI, Ao, and Bo.
125 * It must be at least 1 and at least NSIZE.
126 *
127 * B (workspace) COMPLEX array, dimension (LDA, NSIZE)
128 * Used to hold the matrix whose eigenvalues are to be
129 * computed. On exit, B contains the last matrix actually used.
130 *
131 * AI (workspace) COMPLEX array, dimension (LDA, NSIZE)
132 * Copy of A, modified by CGGEVX.
133 *
134 * BI (workspace) COMPLEX array, dimension (LDA, NSIZE)
135 * Copy of B, modified by CGGEVX.
136 *
137 * ALPHA (workspace) COMPLEX array, dimension (NSIZE)
138 * BETA (workspace) COMPLEX array, dimension (NSIZE)
139 * On exit, ALPHA/BETA are the eigenvalues.
140 *
141 * VL (workspace) COMPLEX array, dimension (LDA, NSIZE)
142 * VL holds the left eigenvectors computed by CGGEVX.
143 *
144 * VR (workspace) COMPLEX array, dimension (LDA, NSIZE)
145 * VR holds the right eigenvectors computed by CGGEVX.
146 *
147 * ILO (output/workspace) INTEGER
148 *
149 * IHI (output/workspace) INTEGER
150 *
151 * LSCALE (output/workspace) REAL array, dimension (N)
152 *
153 * RSCALE (output/workspace) REAL array, dimension (N)
154 *
155 * S (output/workspace) REAL array, dimension (N)
156 *
157 * STRU (output/workspace) REAL array, dimension (N)
158 *
159 * DIF (output/workspace) REAL array, dimension (N)
160 *
161 * DIFTRU (output/workspace) REAL array, dimension (N)
162 *
163 * WORK (workspace) COMPLEX array, dimension (LWORK)
164 *
165 * LWORK (input) INTEGER
166 * Leading dimension of WORK. LWORK >= 2*N*N + 2*N
167 *
168 * RWORK (workspace) REAL array, dimension (6*N)
169 *
170 * IWORK (workspace) INTEGER array, dimension (LIWORK)
171 *
172 * LIWORK (input) INTEGER
173 * Leading dimension of IWORK. LIWORK >= N+2.
174 *
175 * RESULT (output/workspace) REAL array, dimension (4)
176 *
177 * BWORK (workspace) LOGICAL array, dimension (N)
178 *
179 * INFO (output) INTEGER
180 * = 0: successful exit
181 * < 0: if INFO = -i, the i-th argument had an illegal value.
182 * > 0: A routine returned an error code.
183 *
184 * =====================================================================
185 *
186 * .. Parameters ..
187 REAL ZERO, ONE, TEN, TNTH
188 PARAMETER ( ZERO = 0.0E+0, ONE = 1.0E+0, TEN = 1.0E+1,
189 $ TNTH = 1.0E-1 )
190 * ..
191 * .. Local Scalars ..
192 INTEGER I, IPTYPE, IWA, IWB, IWX, IWY, J, LINFO,
193 $ MAXWRK, MINWRK, N, NERRS, NMAX, NPTKNT, NTESTT
194 REAL ABNORM, ANORM, BNORM, RATIO1, RATIO2, THRSH2,
195 $ ULP, ULPINV
196 * ..
197 * .. Local Arrays ..
198 COMPLEX WEIGHT( 5 )
199 * ..
200 * .. External Functions ..
201 INTEGER ILAENV
202 REAL CLANGE, SLAMCH
203 EXTERNAL ILAENV, CLANGE, SLAMCH
204 * ..
205 * .. External Subroutines ..
206 EXTERNAL ALASVM, CGET52, CGGEVX, CLACPY, CLATM6, XERBLA
207 * ..
208 * .. Intrinsic Functions ..
209 INTRINSIC ABS, CMPLX, MAX, SQRT
210 * ..
211 * .. Executable Statements ..
212 *
213 * Check for errors
214 *
215 INFO = 0
216 *
217 NMAX = 5
218 *
219 IF( NSIZE.LT.0 ) THEN
220 INFO = -1
221 ELSE IF( THRESH.LT.ZERO ) THEN
222 INFO = -2
223 ELSE IF( NIN.LE.0 ) THEN
224 INFO = -3
225 ELSE IF( NOUT.LE.0 ) THEN
226 INFO = -4
227 ELSE IF( LDA.LT.1 .OR. LDA.LT.NMAX ) THEN
228 INFO = -6
229 ELSE IF( LIWORK.LT.NMAX+2 ) THEN
230 INFO = -26
231 END IF
232 *
233 * Compute workspace
234 * (Note: Comments in the code beginning "Workspace:" describe the
235 * minimal amount of workspace needed at that point in the code,
236 * as well as the preferred amount for good performance.
237 * NB refers to the optimal block size for the immediately
238 * following subroutine, as returned by ILAENV.)
239 *
240 MINWRK = 1
241 IF( INFO.EQ.0 .AND. LWORK.GE.1 ) THEN
242 MINWRK = 2*NMAX*( NMAX+1 )
243 MAXWRK = NMAX*( 1+ILAENV( 1, 'CGEQRF', ' ', NMAX, 1, NMAX,
244 $ 0 ) )
245 MAXWRK = MAX( MAXWRK, 2*NMAX*( NMAX+1 ) )
246 WORK( 1 ) = MAXWRK
247 END IF
248 *
249 IF( LWORK.LT.MINWRK )
250 $ INFO = -23
251 *
252 IF( INFO.NE.0 ) THEN
253 CALL XERBLA( 'CDRGVX', -INFO )
254 RETURN
255 END IF
256 *
257 N = 5
258 ULP = SLAMCH( 'P' )
259 ULPINV = ONE / ULP
260 THRSH2 = TEN*THRESH
261 NERRS = 0
262 NPTKNT = 0
263 NTESTT = 0
264 *
265 IF( NSIZE.EQ.0 )
266 $ GO TO 90
267 *
268 * Parameters used for generating test matrices.
269 *
270 WEIGHT( 1 ) = CMPLX( SQRT( SQRT( ULP ) ), ZERO )
271 WEIGHT( 2 ) = CMPLX( TNTH, ZERO )
272 WEIGHT( 3 ) = ONE
273 WEIGHT( 4 ) = ONE / WEIGHT( 2 )
274 WEIGHT( 5 ) = ONE / WEIGHT( 1 )
275 *
276 DO 80 IPTYPE = 1, 2
277 DO 70 IWA = 1, 5
278 DO 60 IWB = 1, 5
279 DO 50 IWX = 1, 5
280 DO 40 IWY = 1, 5
281 *
282 * generated a pair of test matrix
283 *
284 CALL CLATM6( IPTYPE, 5, A, LDA, B, VR, LDA, VL,
285 $ LDA, WEIGHT( IWA ), WEIGHT( IWB ),
286 $ WEIGHT( IWX ), WEIGHT( IWY ), STRU,
287 $ DIFTRU )
288 *
289 * Compute eigenvalues/eigenvectors of (A, B).
290 * Compute eigenvalue/eigenvector condition numbers
291 * using computed eigenvectors.
292 *
293 CALL CLACPY( 'F', N, N, A, LDA, AI, LDA )
294 CALL CLACPY( 'F', N, N, B, LDA, BI, LDA )
295 *
296 CALL CGGEVX( 'N', 'V', 'V', 'B', N, AI, LDA, BI,
297 $ LDA, ALPHA, BETA, VL, LDA, VR, LDA,
298 $ ILO, IHI, LSCALE, RSCALE, ANORM,
299 $ BNORM, S, DIF, WORK, LWORK, RWORK,
300 $ IWORK, BWORK, LINFO )
301 IF( LINFO.NE.0 ) THEN
302 WRITE( NOUT, FMT = 9999 )'CGGEVX', LINFO, N,
303 $ IPTYPE, IWA, IWB, IWX, IWY
304 GO TO 30
305 END IF
306 *
307 * Compute the norm(A, B)
308 *
309 CALL CLACPY( 'Full', N, N, AI, LDA, WORK, N )
310 CALL CLACPY( 'Full', N, N, BI, LDA, WORK( N*N+1 ),
311 $ N )
312 ABNORM = CLANGE( 'Fro', N, 2*N, WORK, N, RWORK )
313 *
314 * Tests (1) and (2)
315 *
316 RESULT( 1 ) = ZERO
317 CALL CGET52( .TRUE., N, A, LDA, B, LDA, VL, LDA,
318 $ ALPHA, BETA, WORK, RWORK,
319 $ RESULT( 1 ) )
320 IF( RESULT( 2 ).GT.THRESH ) THEN
321 WRITE( NOUT, FMT = 9998 )'Left', 'CGGEVX',
322 $ RESULT( 2 ), N, IPTYPE, IWA, IWB, IWX, IWY
323 END IF
324 *
325 RESULT( 2 ) = ZERO
326 CALL CGET52( .FALSE., N, A, LDA, B, LDA, VR, LDA,
327 $ ALPHA, BETA, WORK, RWORK,
328 $ RESULT( 2 ) )
329 IF( RESULT( 3 ).GT.THRESH ) THEN
330 WRITE( NOUT, FMT = 9998 )'Right', 'CGGEVX',
331 $ RESULT( 3 ), N, IPTYPE, IWA, IWB, IWX, IWY
332 END IF
333 *
334 * Test (3)
335 *
336 RESULT( 3 ) = ZERO
337 DO 10 I = 1, N
338 IF( S( I ).EQ.ZERO ) THEN
339 IF( STRU( I ).GT.ABNORM*ULP )
340 $ RESULT( 3 ) = ULPINV
341 ELSE IF( STRU( I ).EQ.ZERO ) THEN
342 IF( S( I ).GT.ABNORM*ULP )
343 $ RESULT( 3 ) = ULPINV
344 ELSE
345 RWORK( I ) = MAX( ABS( STRU( I ) / S( I ) ),
346 $ ABS( S( I ) / STRU( I ) ) )
347 RESULT( 3 ) = MAX( RESULT( 3 ), RWORK( I ) )
348 END IF
349 10 CONTINUE
350 *
351 * Test (4)
352 *
353 RESULT( 4 ) = ZERO
354 IF( DIF( 1 ).EQ.ZERO ) THEN
355 IF( DIFTRU( 1 ).GT.ABNORM*ULP )
356 $ RESULT( 4 ) = ULPINV
357 ELSE IF( DIFTRU( 1 ).EQ.ZERO ) THEN
358 IF( DIF( 1 ).GT.ABNORM*ULP )
359 $ RESULT( 4 ) = ULPINV
360 ELSE IF( DIF( 5 ).EQ.ZERO ) THEN
361 IF( DIFTRU( 5 ).GT.ABNORM*ULP )
362 $ RESULT( 4 ) = ULPINV
363 ELSE IF( DIFTRU( 5 ).EQ.ZERO ) THEN
364 IF( DIF( 5 ).GT.ABNORM*ULP )
365 $ RESULT( 4 ) = ULPINV
366 ELSE
367 RATIO1 = MAX( ABS( DIFTRU( 1 ) / DIF( 1 ) ),
368 $ ABS( DIF( 1 ) / DIFTRU( 1 ) ) )
369 RATIO2 = MAX( ABS( DIFTRU( 5 ) / DIF( 5 ) ),
370 $ ABS( DIF( 5 ) / DIFTRU( 5 ) ) )
371 RESULT( 4 ) = MAX( RATIO1, RATIO2 )
372 END IF
373 *
374 NTESTT = NTESTT + 4
375 *
376 * Print out tests which fail.
377 *
378 DO 20 J = 1, 4
379 IF( ( RESULT( J ).GE.THRSH2 .AND. J.GE.4 ) .OR.
380 $ ( RESULT( J ).GE.THRESH .AND. J.LE.3 ) )
381 $ THEN
382 *
383 * If this is the first test to fail,
384 * print a header to the data file.
385 *
386 IF( NERRS.EQ.0 ) THEN
387 WRITE( NOUT, FMT = 9997 )'CXV'
388 *
389 * Print out messages for built-in examples
390 *
391 * Matrix types
392 *
393 WRITE( NOUT, FMT = 9995 )
394 WRITE( NOUT, FMT = 9994 )
395 WRITE( NOUT, FMT = 9993 )
396 *
397 * Tests performed
398 *
399 WRITE( NOUT, FMT = 9992 )'''',
400 $ 'transpose', ''''
401 *
402 END IF
403 NERRS = NERRS + 1
404 IF( RESULT( J ).LT.10000.0 ) THEN
405 WRITE( NOUT, FMT = 9991 )IPTYPE, IWA,
406 $ IWB, IWX, IWY, J, RESULT( J )
407 ELSE
408 WRITE( NOUT, FMT = 9990 )IPTYPE, IWA,
409 $ IWB, IWX, IWY, J, RESULT( J )
410 END IF
411 END IF
412 20 CONTINUE
413 *
414 30 CONTINUE
415 *
416 40 CONTINUE
417 50 CONTINUE
418 60 CONTINUE
419 70 CONTINUE
420 80 CONTINUE
421 *
422 GO TO 150
423 *
424 90 CONTINUE
425 *
426 * Read in data from file to check accuracy of condition estimation
427 * Read input data until N=0
428 *
429 READ( NIN, FMT = *, END = 150 )N
430 IF( N.EQ.0 )
431 $ GO TO 150
432 DO 100 I = 1, N
433 READ( NIN, FMT = * )( A( I, J ), J = 1, N )
434 100 CONTINUE
435 DO 110 I = 1, N
436 READ( NIN, FMT = * )( B( I, J ), J = 1, N )
437 110 CONTINUE
438 READ( NIN, FMT = * )( STRU( I ), I = 1, N )
439 READ( NIN, FMT = * )( DIFTRU( I ), I = 1, N )
440 *
441 NPTKNT = NPTKNT + 1
442 *
443 * Compute eigenvalues/eigenvectors of (A, B).
444 * Compute eigenvalue/eigenvector condition numbers
445 * using computed eigenvectors.
446 *
447 CALL CLACPY( 'F', N, N, A, LDA, AI, LDA )
448 CALL CLACPY( 'F', N, N, B, LDA, BI, LDA )
449 *
450 CALL CGGEVX( 'N', 'V', 'V', 'B', N, AI, LDA, BI, LDA, ALPHA, BETA,
451 $ VL, LDA, VR, LDA, ILO, IHI, LSCALE, RSCALE, ANORM,
452 $ BNORM, S, DIF, WORK, LWORK, RWORK, IWORK, BWORK,
453 $ LINFO )
454 *
455 IF( LINFO.NE.0 ) THEN
456 WRITE( NOUT, FMT = 9987 )'CGGEVX', LINFO, N, NPTKNT
457 GO TO 140
458 END IF
459 *
460 * Compute the norm(A, B)
461 *
462 CALL CLACPY( 'Full', N, N, AI, LDA, WORK, N )
463 CALL CLACPY( 'Full', N, N, BI, LDA, WORK( N*N+1 ), N )
464 ABNORM = CLANGE( 'Fro', N, 2*N, WORK, N, RWORK )
465 *
466 * Tests (1) and (2)
467 *
468 RESULT( 1 ) = ZERO
469 CALL CGET52( .TRUE., N, A, LDA, B, LDA, VL, LDA, ALPHA, BETA,
470 $ WORK, RWORK, RESULT( 1 ) )
471 IF( RESULT( 2 ).GT.THRESH ) THEN
472 WRITE( NOUT, FMT = 9986 )'Left', 'CGGEVX', RESULT( 2 ), N,
473 $ NPTKNT
474 END IF
475 *
476 RESULT( 2 ) = ZERO
477 CALL CGET52( .FALSE., N, A, LDA, B, LDA, VR, LDA, ALPHA, BETA,
478 $ WORK, RWORK, RESULT( 2 ) )
479 IF( RESULT( 3 ).GT.THRESH ) THEN
480 WRITE( NOUT, FMT = 9986 )'Right', 'CGGEVX', RESULT( 3 ), N,
481 $ NPTKNT
482 END IF
483 *
484 * Test (3)
485 *
486 RESULT( 3 ) = ZERO
487 DO 120 I = 1, N
488 IF( S( I ).EQ.ZERO ) THEN
489 IF( STRU( I ).GT.ABNORM*ULP )
490 $ RESULT( 3 ) = ULPINV
491 ELSE IF( STRU( I ).EQ.ZERO ) THEN
492 IF( S( I ).GT.ABNORM*ULP )
493 $ RESULT( 3 ) = ULPINV
494 ELSE
495 RWORK( I ) = MAX( ABS( STRU( I ) / S( I ) ),
496 $ ABS( S( I ) / STRU( I ) ) )
497 RESULT( 3 ) = MAX( RESULT( 3 ), RWORK( I ) )
498 END IF
499 120 CONTINUE
500 *
501 * Test (4)
502 *
503 RESULT( 4 ) = ZERO
504 IF( DIF( 1 ).EQ.ZERO ) THEN
505 IF( DIFTRU( 1 ).GT.ABNORM*ULP )
506 $ RESULT( 4 ) = ULPINV
507 ELSE IF( DIFTRU( 1 ).EQ.ZERO ) THEN
508 IF( DIF( 1 ).GT.ABNORM*ULP )
509 $ RESULT( 4 ) = ULPINV
510 ELSE IF( DIF( 5 ).EQ.ZERO ) THEN
511 IF( DIFTRU( 5 ).GT.ABNORM*ULP )
512 $ RESULT( 4 ) = ULPINV
513 ELSE IF( DIFTRU( 5 ).EQ.ZERO ) THEN
514 IF( DIF( 5 ).GT.ABNORM*ULP )
515 $ RESULT( 4 ) = ULPINV
516 ELSE
517 RATIO1 = MAX( ABS( DIFTRU( 1 ) / DIF( 1 ) ),
518 $ ABS( DIF( 1 ) / DIFTRU( 1 ) ) )
519 RATIO2 = MAX( ABS( DIFTRU( 5 ) / DIF( 5 ) ),
520 $ ABS( DIF( 5 ) / DIFTRU( 5 ) ) )
521 RESULT( 4 ) = MAX( RATIO1, RATIO2 )
522 END IF
523 *
524 NTESTT = NTESTT + 4
525 *
526 * Print out tests which fail.
527 *
528 DO 130 J = 1, 4
529 IF( RESULT( J ).GE.THRSH2 ) THEN
530 *
531 * If this is the first test to fail,
532 * print a header to the data file.
533 *
534 IF( NERRS.EQ.0 ) THEN
535 WRITE( NOUT, FMT = 9997 )'CXV'
536 *
537 * Print out messages for built-in examples
538 *
539 * Matrix types
540 *
541 WRITE( NOUT, FMT = 9996 )
542 *
543 * Tests performed
544 *
545 WRITE( NOUT, FMT = 9992 )'''', 'transpose', ''''
546 *
547 END IF
548 NERRS = NERRS + 1
549 IF( RESULT( J ).LT.10000.0 ) THEN
550 WRITE( NOUT, FMT = 9989 )NPTKNT, N, J, RESULT( J )
551 ELSE
552 WRITE( NOUT, FMT = 9988 )NPTKNT, N, J, RESULT( J )
553 END IF
554 END IF
555 130 CONTINUE
556 *
557 140 CONTINUE
558 *
559 GO TO 90
560 150 CONTINUE
561 *
562 * Summary
563 *
564 CALL ALASVM( 'CXV', NOUT, NERRS, NTESTT, 0 )
565 *
566 WORK( 1 ) = MAXWRK
567 *
568 RETURN
569 *
570 9999 FORMAT( ' CDRGVX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
571 $ I6, ', JTYPE=', I6, ')' )
572 *
573 9998 FORMAT( ' CDRGVX: ', A, ' Eigenvectors from ', A, ' incorrectly ',
574 $ 'normalized.', / ' Bits of error=', 0P, G10.3, ',', 9X,
575 $ 'N=', I6, ', JTYPE=', I6, ', IWA=', I5, ', IWB=', I5,
576 $ ', IWX=', I5, ', IWY=', I5 )
577 *
578 9997 FORMAT( / 1X, A3, ' -- Complex Expert Eigenvalue/vector',
579 $ ' problem driver' )
580 *
581 9996 FORMAT( 'Input Example' )
582 *
583 9995 FORMAT( ' Matrix types: ', / )
584 *
585 9994 FORMAT( ' TYPE 1: Da is diagonal, Db is identity, ',
586 $ / ' A = Y^(-H) Da X^(-1), B = Y^(-H) Db X^(-1) ',
587 $ / ' YH and X are left and right eigenvectors. ', / )
588 *
589 9993 FORMAT( ' TYPE 2: Da is quasi-diagonal, Db is identity, ',
590 $ / ' A = Y^(-H) Da X^(-1), B = Y^(-H) Db X^(-1) ',
591 $ / ' YH and X are left and right eigenvectors. ', / )
592 *
593 9992 FORMAT( / ' Tests performed: ', / 4X,
594 $ ' a is alpha, b is beta, l is a left eigenvector, ', / 4X,
595 $ ' r is a right eigenvector and ', A, ' means ', A, '.',
596 $ / ' 1 = max | ( b A - a B )', A, ' l | / const.',
597 $ / ' 2 = max | ( b A - a B ) r | / const.',
598 $ / ' 3 = max ( Sest/Stru, Stru/Sest ) ',
599 $ ' over all eigenvalues', /
600 $ ' 4 = max( DIFest/DIFtru, DIFtru/DIFest ) ',
601 $ ' over the 1st and 5th eigenvectors', / )
602 *
603 9991 FORMAT( ' Type=', I2, ',', ' IWA=', I2, ', IWB=', I2, ', IWX=',
604 $ I2, ', IWY=', I2, ', result ', I2, ' is', 0P, F8.2 )
605 *
606 9990 FORMAT( ' Type=', I2, ',', ' IWA=', I2, ', IWB=', I2, ', IWX=',
607 $ I2, ', IWY=', I2, ', result ', I2, ' is', 1P, E10.3 )
608 *
609 9989 FORMAT( ' Input example #', I2, ', matrix order=', I4, ',',
610 $ ' result ', I2, ' is', 0P, F8.2 )
611 *
612 9988 FORMAT( ' Input example #', I2, ', matrix order=', I4, ',',
613 $ ' result ', I2, ' is', 1P, E10.3 )
614 *
615 9987 FORMAT( ' CDRGVX: ', A, ' returned INFO=', I6, '.', / 9X, 'N=',
616 $ I6, ', Input example #', I2, ')' )
617 *
618 9986 FORMAT( ' CDRGVX: ', A, ' Eigenvectors from ', A, ' incorrectly ',
619 $ 'normalized.', / ' Bits of error=', 0P, G10.3, ',', 9X,
620 $ 'N=', I6, ', Input Example #', I2, ')' )
621 *
622 * End of CDRGVX
623 *
624 END