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