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