LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
ddrgvx.f
Go to the documentation of this file.
1 *> \brief \b DDRGVX
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 * Definition:
9 * ===========
10 *
11 * SUBROUTINE DDRGVX( NSIZE, THRESH, NIN, NOUT, A, LDA, B, AI, BI,
12 * ALPHAR, ALPHAI, BETA, VL, VR, ILO, IHI, LSCALE,
13 * RSCALE, S, DTRU, DIF, DIFTRU, WORK, LWORK,
14 * IWORK, LIWORK, RESULT, BWORK, INFO )
15 *
16 * .. Scalar Arguments ..
17 * INTEGER IHI, ILO, INFO, LDA, LIWORK, LWORK, NIN, NOUT,
18 * $ NSIZE
19 * DOUBLE PRECISION THRESH
20 * ..
21 * .. Array Arguments ..
22 * LOGICAL BWORK( * )
23 * INTEGER IWORK( * )
24 * DOUBLE PRECISION A( LDA, * ), AI( LDA, * ), ALPHAI( * ),
25 * $ ALPHAR( * ), B( LDA, * ), BETA( * ),
26 * $ BI( LDA, * ), DIF( * ), DIFTRU( * ), DTRU( * ),
27 * $ LSCALE( * ), RESULT( 4 ), RSCALE( * ), S( * ),
28 * $ VL( LDA, * ), VR( LDA, * ), WORK( * )
29 * ..
30 *
31 *
32 *> \par Purpose:
33 * =============
34 *>
35 *> \verbatim
36 *>
37 *> DDRGVX checks the nonsymmetric generalized eigenvalue problem
38 *> expert driver DGGEVX.
39 *>
40 *> DGGEVX computes the generalized eigenvalues, (optionally) the left
41 *> and/or right eigenvectors, (optionally) computes a balancing
42 *> transformation to improve the conditioning, and (optionally)
43 *> reciprocal condition numbers for the eigenvalues and eigenvectors.
44 *>
45 *> When DDRGVX is called with NSIZE > 0, two types of test matrix pairs
46 *> are generated by the subroutine DLATM6 and test the driver DGGEVX.
47 *> The test matrices have the known exact condition numbers for
48 *> eigenvalues. For the condition numbers of the eigenvectors
49 *> corresponding the first and last eigenvalues are also know
50 *> ``exactly'' (see DLATM6).
51 *>
52 *> For each matrix pair, the following tests will be performed and
53 *> compared with the threshold THRESH.
54 *>
55 *> (1) max over all left eigenvalue/-vector pairs (beta/alpha,l) of
56 *>
57 *> | l**H * (beta A - alpha B) | / ( ulp max( |beta A|, |alpha B| ) )
58 *>
59 *> where l**H is the conjugate tranpose of l.
60 *>
61 *> (2) max over all right eigenvalue/-vector pairs (beta/alpha,r) of
62 *>
63 *> | (beta A - alpha B) r | / ( ulp max( |beta A|, |alpha B| ) )
64 *>
65 *> (3) The condition number S(i) of eigenvalues computed by DGGEVX
66 *> differs less than a factor THRESH from the exact S(i) (see
67 *> DLATM6).
68 *>
69 *> (4) DIF(i) computed by DTGSNA differs less than a factor 10*THRESH
70 *> from the exact value (for the 1st and 5th vectors only).
71 *>
72 *> Test Matrices
73 *> =============
74 *>
75 *> Two kinds of test matrix pairs
76 *>
77 *> (A, B) = inverse(YH) * (Da, Db) * inverse(X)
78 *>
79 *> are used in the tests:
80 *>
81 *> 1: Da = 1+a 0 0 0 0 Db = 1 0 0 0 0
82 *> 0 2+a 0 0 0 0 1 0 0 0
83 *> 0 0 3+a 0 0 0 0 1 0 0
84 *> 0 0 0 4+a 0 0 0 0 1 0
85 *> 0 0 0 0 5+a , 0 0 0 0 1 , and
86 *>
87 *> 2: Da = 1 -1 0 0 0 Db = 1 0 0 0 0
88 *> 1 1 0 0 0 0 1 0 0 0
89 *> 0 0 1 0 0 0 0 1 0 0
90 *> 0 0 0 1+a 1+b 0 0 0 1 0
91 *> 0 0 0 -1-b 1+a , 0 0 0 0 1 .
92 *>
93 *> In both cases the same inverse(YH) and inverse(X) are used to compute
94 *> (A, B), giving the exact eigenvectors to (A,B) as (YH, X):
95 *>
96 *> YH: = 1 0 -y y -y X = 1 0 -x -x x
97 *> 0 1 -y y -y 0 1 x -x -x
98 *> 0 0 1 0 0 0 0 1 0 0
99 *> 0 0 0 1 0 0 0 0 1 0
100 *> 0 0 0 0 1, 0 0 0 0 1 , where
101 *>
102 *> a, b, x and y will have all values independently of each other from
103 *> { sqrt(sqrt(ULP)), 0.1, 1, 10, 1/sqrt(sqrt(ULP)) }.
104 *> \endverbatim
105 *
106 * Arguments:
107 * ==========
108 *
109 *> \param[in] NSIZE
110 *> \verbatim
111 *> NSIZE is INTEGER
112 *> The number of sizes of matrices to use. NSIZE must be at
113 *> least zero. If it is zero, no randomly generated matrices
114 *> are tested, but any test matrices read from NIN will be
115 *> tested.
116 *> \endverbatim
117 *>
118 *> \param[in] THRESH
119 *> \verbatim
120 *> THRESH is DOUBLE PRECISION
121 *> A test will count as "failed" if the "error", computed as
122 *> described above, exceeds THRESH. Note that the error
123 *> is scaled to be O(1), so THRESH should be a reasonably
124 *> small multiple of 1, e.g., 10 or 100. In particular,
125 *> it should not depend on the precision (single vs. double)
126 *> or the size of the matrix. It must be at least zero.
127 *> \endverbatim
128 *>
129 *> \param[in] NIN
130 *> \verbatim
131 *> NIN is INTEGER
132 *> The FORTRAN unit number for reading in the data file of
133 *> problems to solve.
134 *> \endverbatim
135 *>
136 *> \param[in] NOUT
137 *> \verbatim
138 *> NOUT is INTEGER
139 *> The FORTRAN unit number for printing out error messages
140 *> (e.g., if a routine returns IINFO not equal to 0.)
141 *> \endverbatim
142 *>
143 *> \param[out] A
144 *> \verbatim
145 *> A is DOUBLE PRECISION array, dimension (LDA, NSIZE)
146 *> Used to hold the matrix whose eigenvalues are to be
147 *> computed. On exit, A contains the last matrix actually used.
148 *> \endverbatim
149 *>
150 *> \param[in] LDA
151 *> \verbatim
152 *> LDA is INTEGER
153 *> The leading dimension of A, B, AI, BI, Ao, and Bo.
154 *> It must be at least 1 and at least NSIZE.
155 *> \endverbatim
156 *>
157 *> \param[out] B
158 *> \verbatim
159 *> B is DOUBLE PRECISION array, dimension (LDA, NSIZE)
160 *> Used to hold the matrix whose eigenvalues are to be
161 *> computed. On exit, B contains the last matrix actually used.
162 *> \endverbatim
163 *>
164 *> \param[out] AI
165 *> \verbatim
166 *> AI is DOUBLE PRECISION array, dimension (LDA, NSIZE)
167 *> Copy of A, modified by DGGEVX.
168 *> \endverbatim
169 *>
170 *> \param[out] BI
171 *> \verbatim
172 *> BI is DOUBLE PRECISION array, dimension (LDA, NSIZE)
173 *> Copy of B, modified by DGGEVX.
174 *> \endverbatim
175 *>
176 *> \param[out] ALPHAR
177 *> \verbatim
178 *> ALPHAR is DOUBLE PRECISION array, dimension (NSIZE)
179 *> \endverbatim
180 *>
181 *> \param[out] ALPHAI
182 *> \verbatim
183 *> ALPHAI is DOUBLE PRECISION array, dimension (NSIZE)
184 *> \endverbatim
185 *>
186 *> \param[out] BETA
187 *> \verbatim
188 *> BETA is DOUBLE PRECISION array, dimension (NSIZE)
189 *>
190 *> On exit, (ALPHAR + ALPHAI*i)/BETA are the eigenvalues.
191 *> \endverbatim
192 *>
193 *> \param[out] VL
194 *> \verbatim
195 *> VL is DOUBLE PRECISION array, dimension (LDA, NSIZE)
196 *> VL holds the left eigenvectors computed by DGGEVX.
197 *> \endverbatim
198 *>
199 *> \param[out] VR
200 *> \verbatim
201 *> VR is DOUBLE PRECISION array, dimension (LDA, NSIZE)
202 *> VR holds the right eigenvectors computed by DGGEVX.
203 *> \endverbatim
204 *>
205 *> \param[out] ILO
206 *> \verbatim
207 *> ILO is INTEGER
208 *> \endverbatim
209 *>
210 *> \param[out] IHI
211 *> \verbatim
212 *> IHI is INTEGER
213 *> \endverbatim
214 *>
215 *> \param[out] LSCALE
216 *> \verbatim
217 *> LSCALE is DOUBLE PRECISION array, dimension (N)
218 *> \endverbatim
219 *>
220 *> \param[out] RSCALE
221 *> \verbatim
222 *> RSCALE is DOUBLE PRECISION array, dimension (N)
223 *> \endverbatim
224 *>
225 *> \param[out] S
226 *> \verbatim
227 *> S is DOUBLE PRECISION array, dimension (N)
228 *> \endverbatim
229 *>
230 *> \param[out] DTRU
231 *> \verbatim
232 *> DTRU is DOUBLE PRECISION array, dimension (N)
233 *> \endverbatim
234 *>
235 *> \param[out] DIF
236 *> \verbatim
237 *> DIF is DOUBLE PRECISION array, dimension (N)
238 *> \endverbatim
239 *>
240 *> \param[out] DIFTRU
241 *> \verbatim
242 *> DIFTRU is DOUBLE PRECISION array, dimension (N)
243 *> \endverbatim
244 *>
245 *> \param[out] WORK
246 *> \verbatim
247 *> WORK is DOUBLE PRECISION array, dimension (LWORK)
248 *> \endverbatim
249 *>
250 *> \param[in] LWORK
251 *> \verbatim
252 *> LWORK is INTEGER
253 *> Leading dimension of WORK. LWORK >= 2*N*N+12*N+16.
254 *> \endverbatim
255 *>
256 *> \param[out] IWORK
257 *> \verbatim
258 *> IWORK is INTEGER array, dimension (LIWORK)
259 *> \endverbatim
260 *>
261 *> \param[in] LIWORK
262 *> \verbatim
263 *> LIWORK is INTEGER
264 *> Leading dimension of IWORK. Must be at least N+6.
265 *> \endverbatim
266 *>
267 *> \param[out] RESULT
268 *> \verbatim
269 *> RESULT is DOUBLE PRECISION array, dimension (4)
270 *> \endverbatim
271 *>
272 *> \param[out] BWORK
273 *> \verbatim
274 *> BWORK is LOGICAL array, dimension (N)
275 *> \endverbatim
276 *>
277 *> \param[out] INFO
278 *> \verbatim
279 *> INFO is INTEGER
280 *> = 0: successful exit
281 *> < 0: if INFO = -i, the i-th argument had an illegal value.
282 *> > 0: A routine returned an error code.
283 *> \endverbatim
284 *
285 * Authors:
286 * ========
287 *
288 *> \author Univ. of Tennessee
289 *> \author Univ. of California Berkeley
290 *> \author Univ. of Colorado Denver
291 *> \author NAG Ltd.
292 *
293 *> \date April 2012
294 *
295 *> \ingroup double_eig
296 *
297 * =====================================================================
298  SUBROUTINE ddrgvx( NSIZE, THRESH, NIN, NOUT, A, LDA, B, AI, BI,
299  $ alphar, alphai, beta, vl, vr, ilo, ihi, lscale,
300  $ rscale, s, dtru, dif, diftru, work, lwork,
301  $ iwork, liwork, result, bwork, info )
302 *
303 * -- LAPACK test routine (version 3.6.1) --
304 * -- LAPACK is a software package provided by Univ. of Tennessee, --
305 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
306 * April 2012
307 *
308 * .. Scalar Arguments ..
309  INTEGER IHI, ILO, INFO, LDA, LIWORK, LWORK, NIN, NOUT,
310  $ nsize
311  DOUBLE PRECISION THRESH
312 * ..
313 * .. Array Arguments ..
314  LOGICAL BWORK( * )
315  INTEGER IWORK( * )
316  DOUBLE PRECISION A( lda, * ), AI( lda, * ), ALPHAI( * ),
317  $ alphar( * ), b( lda, * ), beta( * ),
318  $ bi( lda, * ), dif( * ), diftru( * ), dtru( * ),
319  $ lscale( * ), result( 4 ), rscale( * ), s( * ),
320  $ vl( lda, * ), vr( lda, * ), work( * )
321 * ..
322 *
323 * =====================================================================
324 *
325 * .. Parameters ..
326  DOUBLE PRECISION ZERO, ONE, TEN, TNTH, HALF
327  parameter ( zero = 0.0d+0, one = 1.0d+0, ten = 1.0d+1,
328  $ tnth = 1.0d-1, half = 0.5d+0 )
329 * ..
330 * .. Local Scalars ..
331  INTEGER I, IPTYPE, IWA, IWB, IWX, IWY, J, LINFO,
332  $ maxwrk, minwrk, n, nerrs, nmax, nptknt, ntestt
333  DOUBLE PRECISION ABNORM, ANORM, BNORM, RATIO1, RATIO2, THRSH2,
334  $ ulp, ulpinv
335 * ..
336 * .. Local Arrays ..
337  DOUBLE PRECISION WEIGHT( 5 )
338 * ..
339 * .. External Functions ..
340  INTEGER ILAENV
341  DOUBLE PRECISION DLAMCH, DLANGE
342  EXTERNAL ilaenv, dlamch, dlange
343 * ..
344 * .. External Subroutines ..
345  EXTERNAL alasvm, dget52, dggevx, dlacpy, dlatm6, xerbla
346 * ..
347 * .. Intrinsic Functions ..
348  INTRINSIC abs, max, sqrt
349 * ..
350 * .. Executable Statements ..
351 *
352 * Check for errors
353 *
354  info = 0
355 *
356  nmax = 5
357 *
358  IF( nsize.LT.0 ) THEN
359  info = -1
360  ELSE IF( thresh.LT.zero ) THEN
361  info = -2
362  ELSE IF( nin.LE.0 ) THEN
363  info = -3
364  ELSE IF( nout.LE.0 ) THEN
365  info = -4
366  ELSE IF( lda.LT.1 .OR. lda.LT.nmax ) THEN
367  info = -6
368  ELSE IF( liwork.LT.nmax+6 ) THEN
369  info = -26
370  END IF
371 *
372 * Compute workspace
373 * (Note: Comments in the code beginning "Workspace:" describe the
374 * minimal amount of workspace needed at that point in the code,
375 * as well as the preferred amount for good performance.
376 * NB refers to the optimal block size for the immediately
377 * following subroutine, as returned by ILAENV.)
378 *
379  minwrk = 1
380  IF( info.EQ.0 .AND. lwork.GE.1 ) THEN
381  minwrk = 2*nmax*nmax + 12*nmax + 16
382  maxwrk = 6*nmax + nmax*ilaenv( 1, 'DGEQRF', ' ', nmax, 1, nmax,
383  $ 0 )
384  maxwrk = max( maxwrk, 2*nmax*nmax+12*nmax+16 )
385  work( 1 ) = maxwrk
386  END IF
387 *
388  IF( lwork.LT.minwrk )
389  $ info = -24
390 *
391  IF( info.NE.0 ) THEN
392  CALL xerbla( 'DDRGVX', -info )
393  RETURN
394  END IF
395 *
396  n = 5
397  ulp = dlamch( 'P' )
398  ulpinv = one / ulp
399  thrsh2 = ten*thresh
400  nerrs = 0
401  nptknt = 0
402  ntestt = 0
403 *
404  IF( nsize.EQ.0 )
405  $ GO TO 90
406 *
407 * Parameters used for generating test matrices.
408 *
409  weight( 1 ) = tnth
410  weight( 2 ) = half
411  weight( 3 ) = one
412  weight( 4 ) = one / weight( 2 )
413  weight( 5 ) = one / weight( 1 )
414 *
415  DO 80 iptype = 1, 2
416  DO 70 iwa = 1, 5
417  DO 60 iwb = 1, 5
418  DO 50 iwx = 1, 5
419  DO 40 iwy = 1, 5
420 *
421 * generated a test matrix pair
422 *
423  CALL dlatm6( iptype, 5, a, lda, b, vr, lda, vl,
424  $ lda, weight( iwa ), weight( iwb ),
425  $ weight( iwx ), weight( iwy ), dtru,
426  $ diftru )
427 *
428 * Compute eigenvalues/eigenvectors of (A, B).
429 * Compute eigenvalue/eigenvector condition numbers
430 * using computed eigenvectors.
431 *
432  CALL dlacpy( 'F', n, n, a, lda, ai, lda )
433  CALL dlacpy( 'F', n, n, b, lda, bi, lda )
434 *
435  CALL dggevx( 'N', 'V', 'V', 'B', n, ai, lda, bi,
436  $ lda, alphar, alphai, beta, vl, lda,
437  $ vr, lda, ilo, ihi, lscale, rscale,
438  $ anorm, bnorm, s, dif, work, lwork,
439  $ iwork, bwork, linfo )
440  IF( linfo.NE.0 ) THEN
441  result( 1 ) = ulpinv
442  WRITE( nout, fmt = 9999 )'DGGEVX', linfo, n,
443  $ iptype
444  GO TO 30
445  END IF
446 *
447 * Compute the norm(A, B)
448 *
449  CALL dlacpy( 'Full', n, n, ai, lda, work, n )
450  CALL dlacpy( 'Full', n, n, bi, lda, work( n*n+1 ),
451  $ n )
452  abnorm = dlange( 'Fro', n, 2*n, work, n, work )
453 *
454 * Tests (1) and (2)
455 *
456  result( 1 ) = zero
457  CALL dget52( .true., n, a, lda, b, lda, vl, lda,
458  $ alphar, alphai, beta, work,
459  $ result( 1 ) )
460  IF( result( 2 ).GT.thresh ) THEN
461  WRITE( nout, fmt = 9998 )'Left', 'DGGEVX',
462  $ result( 2 ), n, iptype, iwa, iwb, iwx, iwy
463  END IF
464 *
465  result( 2 ) = zero
466  CALL dget52( .false., n, a, lda, b, lda, vr, lda,
467  $ alphar, alphai, beta, work,
468  $ result( 2 ) )
469  IF( result( 3 ).GT.thresh ) THEN
470  WRITE( nout, fmt = 9998 )'Right', 'DGGEVX',
471  $ result( 3 ), n, iptype, iwa, iwb, iwx, iwy
472  END IF
473 *
474 * Test (3)
475 *
476  result( 3 ) = zero
477  DO 10 i = 1, n
478  IF( s( i ).EQ.zero ) THEN
479  IF( dtru( i ).GT.abnorm*ulp )
480  $ result( 3 ) = ulpinv
481  ELSE IF( dtru( i ).EQ.zero ) THEN
482  IF( s( i ).GT.abnorm*ulp )
483  $ result( 3 ) = ulpinv
484  ELSE
485  work( i ) = max( abs( dtru( i ) / s( i ) ),
486  $ abs( s( i ) / dtru( i ) ) )
487  result( 3 ) = max( result( 3 ), work( i ) )
488  END IF
489  10 CONTINUE
490 *
491 * Test (4)
492 *
493  result( 4 ) = zero
494  IF( dif( 1 ).EQ.zero ) THEN
495  IF( diftru( 1 ).GT.abnorm*ulp )
496  $ result( 4 ) = ulpinv
497  ELSE IF( diftru( 1 ).EQ.zero ) THEN
498  IF( dif( 1 ).GT.abnorm*ulp )
499  $ result( 4 ) = ulpinv
500  ELSE IF( dif( 5 ).EQ.zero ) THEN
501  IF( diftru( 5 ).GT.abnorm*ulp )
502  $ result( 4 ) = ulpinv
503  ELSE IF( diftru( 5 ).EQ.zero ) THEN
504  IF( dif( 5 ).GT.abnorm*ulp )
505  $ result( 4 ) = ulpinv
506  ELSE
507  ratio1 = max( abs( diftru( 1 ) / dif( 1 ) ),
508  $ abs( dif( 1 ) / diftru( 1 ) ) )
509  ratio2 = max( abs( diftru( 5 ) / dif( 5 ) ),
510  $ abs( dif( 5 ) / diftru( 5 ) ) )
511  result( 4 ) = max( ratio1, ratio2 )
512  END IF
513 *
514  ntestt = ntestt + 4
515 *
516 * Print out tests which fail.
517 *
518  DO 20 j = 1, 4
519  IF( ( result( j ).GE.thrsh2 .AND. j.GE.4 ) .OR.
520  $ ( result( j ).GE.thresh .AND. j.LE.3 ) )
521  $ THEN
522 *
523 * If this is the first test to fail,
524 * print a header to the data file.
525 *
526  IF( nerrs.EQ.0 ) THEN
527  WRITE( nout, fmt = 9997 )'DXV'
528 *
529 * Print out messages for built-in examples
530 *
531 * Matrix types
532 *
533  WRITE( nout, fmt = 9995 )
534  WRITE( nout, fmt = 9994 )
535  WRITE( nout, fmt = 9993 )
536 *
537 * Tests performed
538 *
539  WRITE( nout, fmt = 9992 )'''',
540  $ 'transpose', ''''
541 *
542  END IF
543  nerrs = nerrs + 1
544  IF( result( j ).LT.10000.0d0 ) THEN
545  WRITE( nout, fmt = 9991 )iptype, iwa,
546  $ iwb, iwx, iwy, j, result( j )
547  ELSE
548  WRITE( nout, fmt = 9990 )iptype, iwa,
549  $ iwb, iwx, iwy, j, result( j )
550  END IF
551  END IF
552  20 CONTINUE
553 *
554  30 CONTINUE
555 *
556  40 CONTINUE
557  50 CONTINUE
558  60 CONTINUE
559  70 CONTINUE
560  80 CONTINUE
561 *
562  GO TO 150
563 *
564  90 CONTINUE
565 *
566 * Read in data from file to check accuracy of condition estimation
567 * Read input data until N=0
568 *
569  READ( nin, fmt = *, end = 150 )n
570  IF( n.EQ.0 )
571  $ GO TO 150
572  DO 100 i = 1, n
573  READ( nin, fmt = * )( a( i, j ), j = 1, n )
574  100 CONTINUE
575  DO 110 i = 1, n
576  READ( nin, fmt = * )( b( i, j ), j = 1, n )
577  110 CONTINUE
578  READ( nin, fmt = * )( dtru( i ), i = 1, n )
579  READ( nin, fmt = * )( diftru( i ), i = 1, n )
580 *
581  nptknt = nptknt + 1
582 *
583 * Compute eigenvalues/eigenvectors of (A, B).
584 * Compute eigenvalue/eigenvector condition numbers
585 * using computed eigenvectors.
586 *
587  CALL dlacpy( 'F', n, n, a, lda, ai, lda )
588  CALL dlacpy( 'F', n, n, b, lda, bi, lda )
589 *
590  CALL dggevx( 'N', 'V', 'V', 'B', n, ai, lda, bi, lda, alphar,
591  $ alphai, beta, vl, lda, vr, lda, ilo, ihi, lscale,
592  $ rscale, anorm, bnorm, s, dif, work, lwork, iwork,
593  $ bwork, linfo )
594 *
595  IF( linfo.NE.0 ) THEN
596  result( 1 ) = ulpinv
597  WRITE( nout, fmt = 9987 )'DGGEVX', linfo, n, nptknt
598  GO TO 140
599  END IF
600 *
601 * Compute the norm(A, B)
602 *
603  CALL dlacpy( 'Full', n, n, ai, lda, work, n )
604  CALL dlacpy( 'Full', n, n, bi, lda, work( n*n+1 ), n )
605  abnorm = dlange( 'Fro', n, 2*n, work, n, work )
606 *
607 * Tests (1) and (2)
608 *
609  result( 1 ) = zero
610  CALL dget52( .true., n, a, lda, b, lda, vl, lda, alphar, alphai,
611  $ beta, work, result( 1 ) )
612  IF( result( 2 ).GT.thresh ) THEN
613  WRITE( nout, fmt = 9986 )'Left', 'DGGEVX', result( 2 ), n,
614  $ nptknt
615  END IF
616 *
617  result( 2 ) = zero
618  CALL dget52( .false., n, a, lda, b, lda, vr, lda, alphar, alphai,
619  $ beta, work, result( 2 ) )
620  IF( result( 3 ).GT.thresh ) THEN
621  WRITE( nout, fmt = 9986 )'Right', 'DGGEVX', result( 3 ), n,
622  $ nptknt
623  END IF
624 *
625 * Test (3)
626 *
627  result( 3 ) = zero
628  DO 120 i = 1, n
629  IF( s( i ).EQ.zero ) THEN
630  IF( dtru( i ).GT.abnorm*ulp )
631  $ result( 3 ) = ulpinv
632  ELSE IF( dtru( i ).EQ.zero ) THEN
633  IF( s( i ).GT.abnorm*ulp )
634  $ result( 3 ) = ulpinv
635  ELSE
636  work( i ) = max( abs( dtru( i ) / s( i ) ),
637  $ abs( s( i ) / dtru( i ) ) )
638  result( 3 ) = max( result( 3 ), work( i ) )
639  END IF
640  120 CONTINUE
641 *
642 * Test (4)
643 *
644  result( 4 ) = zero
645  IF( dif( 1 ).EQ.zero ) THEN
646  IF( diftru( 1 ).GT.abnorm*ulp )
647  $ result( 4 ) = ulpinv
648  ELSE IF( diftru( 1 ).EQ.zero ) THEN
649  IF( dif( 1 ).GT.abnorm*ulp )
650  $ result( 4 ) = ulpinv
651  ELSE IF( dif( 5 ).EQ.zero ) THEN
652  IF( diftru( 5 ).GT.abnorm*ulp )
653  $ result( 4 ) = ulpinv
654  ELSE IF( diftru( 5 ).EQ.zero ) THEN
655  IF( dif( 5 ).GT.abnorm*ulp )
656  $ result( 4 ) = ulpinv
657  ELSE
658  ratio1 = max( abs( diftru( 1 ) / dif( 1 ) ),
659  $ abs( dif( 1 ) / diftru( 1 ) ) )
660  ratio2 = max( abs( diftru( 5 ) / dif( 5 ) ),
661  $ abs( dif( 5 ) / diftru( 5 ) ) )
662  result( 4 ) = max( ratio1, ratio2 )
663  END IF
664 *
665  ntestt = ntestt + 4
666 *
667 * Print out tests which fail.
668 *
669  DO 130 j = 1, 4
670  IF( result( j ).GE.thrsh2 ) THEN
671 *
672 * If this is the first test to fail,
673 * print a header to the data file.
674 *
675  IF( nerrs.EQ.0 ) THEN
676  WRITE( nout, fmt = 9997 )'DXV'
677 *
678 * Print out messages for built-in examples
679 *
680 * Matrix types
681 *
682  WRITE( nout, fmt = 9996 )
683 *
684 * Tests performed
685 *
686  WRITE( nout, fmt = 9992 )'''', 'transpose', ''''
687 *
688  END IF
689  nerrs = nerrs + 1
690  IF( result( j ).LT.10000.0d0 ) THEN
691  WRITE( nout, fmt = 9989 )nptknt, n, j, result( j )
692  ELSE
693  WRITE( nout, fmt = 9988 )nptknt, n, j, result( j )
694  END IF
695  END IF
696  130 CONTINUE
697 *
698  140 CONTINUE
699 *
700  GO TO 90
701  150 CONTINUE
702 *
703 * Summary
704 *
705  CALL alasvm( 'DXV', nout, nerrs, ntestt, 0 )
706 *
707  work( 1 ) = maxwrk
708 *
709  RETURN
710 *
711  9999 FORMAT( ' DDRGVX: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
712  $ i6, ', JTYPE=', i6, ')' )
713 *
714  9998 FORMAT( ' DDRGVX: ', a, ' Eigenvectors from ', a, ' incorrectly ',
715  $ 'normalized.', / ' Bits of error=', 0p, g10.3, ',', 9x,
716  $ 'N=', i6, ', JTYPE=', i6, ', IWA=', i5, ', IWB=', i5,
717  $ ', IWX=', i5, ', IWY=', i5 )
718 *
719  9997 FORMAT( / 1x, a3, ' -- Real Expert Eigenvalue/vector',
720  $ ' problem driver' )
721 *
722  9996 FORMAT( ' Input Example' )
723 *
724  9995 FORMAT( ' Matrix types: ', / )
725 *
726  9994 FORMAT( ' TYPE 1: Da is diagonal, Db is identity, ',
727  $ / ' A = Y^(-H) Da X^(-1), B = Y^(-H) Db X^(-1) ',
728  $ / ' YH and X are left and right eigenvectors. ', / )
729 *
730  9993 FORMAT( ' TYPE 2: Da is quasi-diagonal, Db is identity, ',
731  $ / ' A = Y^(-H) Da X^(-1), B = Y^(-H) Db X^(-1) ',
732  $ / ' YH and X are left and right eigenvectors. ', / )
733 *
734  9992 FORMAT( / ' Tests performed: ', / 4x,
735  $ ' a is alpha, b is beta, l is a left eigenvector, ', / 4x,
736  $ ' r is a right eigenvector and ', a, ' means ', a, '.',
737  $ / ' 1 = max | ( b A - a B )', a, ' l | / const.',
738  $ / ' 2 = max | ( b A - a B ) r | / const.',
739  $ / ' 3 = max ( Sest/Stru, Stru/Sest ) ',
740  $ ' over all eigenvalues', /
741  $ ' 4 = max( DIFest/DIFtru, DIFtru/DIFest ) ',
742  $ ' over the 1st and 5th eigenvectors', / )
743 *
744  9991 FORMAT( ' Type=', i2, ',', ' IWA=', i2, ', IWB=', i2, ', IWX=',
745  $ i2, ', IWY=', i2, ', result ', i2, ' is', 0p, f8.2 )
746  9990 FORMAT( ' Type=', i2, ',', ' IWA=', i2, ', IWB=', i2, ', IWX=',
747  $ i2, ', IWY=', i2, ', result ', i2, ' is', 1p, d10.3 )
748  9989 FORMAT( ' Input example #', i2, ', matrix order=', i4, ',',
749  $ ' result ', i2, ' is', 0p, f8.2 )
750  9988 FORMAT( ' Input example #', i2, ', matrix order=', i4, ',',
751  $ ' result ', i2, ' is', 1p, d10.3 )
752  9987 FORMAT( ' DDRGVX: ', a, ' returned INFO=', i6, '.', / 9x, 'N=',
753  $ i6, ', Input example #', i2, ')' )
754 *
755  9986 FORMAT( ' DDRGVX: ', a, ' Eigenvectors from ', a, ' incorrectly ',
756  $ 'normalized.', / ' Bits of error=', 0p, g10.3, ',', 9x,
757  $ 'N=', i6, ', Input Example #', i2, ')' )
758 *
759 *
760 * End of DDRGVX
761 *
762  END
subroutine alasvm(TYPE, NOUT, NFAIL, NRUN, NERRS)
ALASVM
Definition: alasvm.f:75
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:105
subroutine dlatm6(TYPE, N, A, LDA, B, X, LDX, Y, LDY, ALPHA, BETA, WX, WY, S, DIF)
DLATM6
Definition: dlatm6.f:178
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine ddrgvx(NSIZE, THRESH, NIN, NOUT, A, LDA, B, AI, BI, ALPHAR, ALPHAI, BETA, VL, VR, ILO, IHI, LSCALE, RSCALE, S, DTRU, DIF, DIFTRU, WORK, LWORK, IWORK, LIWORK, RESULT, BWORK, INFO)
DDRGVX
Definition: ddrgvx.f:302
subroutine dget52(LEFT, N, A, LDA, B, LDB, E, LDE, ALPHAR, ALPHAI, BETA, WORK, RESULT)
DGET52
Definition: dget52.f:201
subroutine dggevx(BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, B, LDB, ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR, ILO, IHI, LSCALE, RSCALE, ABNRM, BBNRM, RCONDE, RCONDV, WORK, LWORK, IWORK, BWORK, INFO)
DGGEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices ...
Definition: dggevx.f:393