LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
dgegv.f
Go to the documentation of this file.
1 *> \brief <b> DGEEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices</b>
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download DGEGV + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgegv.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgegv.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgegv.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DGEGV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR, ALPHAI,
22 * BETA, VL, LDVL, VR, LDVR, WORK, LWORK, INFO )
23 *
24 * .. Scalar Arguments ..
25 * CHARACTER JOBVL, JOBVR
26 * INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, N
27 * ..
28 * .. Array Arguments ..
29 * DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
30 * $ B( LDB, * ), BETA( * ), VL( LDVL, * ),
31 * $ VR( LDVR, * ), WORK( * )
32 * ..
33 *
34 *
35 *> \par Purpose:
36 * =============
37 *>
38 *> \verbatim
39 *>
40 *> This routine is deprecated and has been replaced by routine DGGEV.
41 *>
42 *> DGEGV computes the eigenvalues and, optionally, the left and/or right
43 *> eigenvectors of a real matrix pair (A,B).
44 *> Given two square matrices A and B,
45 *> the generalized nonsymmetric eigenvalue problem (GNEP) is to find the
46 *> eigenvalues lambda and corresponding (non-zero) eigenvectors x such
47 *> that
48 *>
49 *> A*x = lambda*B*x.
50 *>
51 *> An alternate form is to find the eigenvalues mu and corresponding
52 *> eigenvectors y such that
53 *>
54 *> mu*A*y = B*y.
55 *>
56 *> These two forms are equivalent with mu = 1/lambda and x = y if
57 *> neither lambda nor mu is zero. In order to deal with the case that
58 *> lambda or mu is zero or small, two values alpha and beta are returned
59 *> for each eigenvalue, such that lambda = alpha/beta and
60 *> mu = beta/alpha.
61 *>
62 *> The vectors x and y in the above equations are right eigenvectors of
63 *> the matrix pair (A,B). Vectors u and v satisfying
64 *>
65 *> u**H*A = lambda*u**H*B or mu*v**H*A = v**H*B
66 *>
67 *> are left eigenvectors of (A,B).
68 *>
69 *> Note: this routine performs "full balancing" on A and B
70 *> \endverbatim
71 *
72 * Arguments:
73 * ==========
74 *
75 *> \param[in] JOBVL
76 *> \verbatim
77 *> JOBVL is CHARACTER*1
78 *> = 'N': do not compute the left generalized eigenvectors;
79 *> = 'V': compute the left generalized eigenvectors (returned
80 *> in VL).
81 *> \endverbatim
82 *>
83 *> \param[in] JOBVR
84 *> \verbatim
85 *> JOBVR is CHARACTER*1
86 *> = 'N': do not compute the right generalized eigenvectors;
87 *> = 'V': compute the right generalized eigenvectors (returned
88 *> in VR).
89 *> \endverbatim
90 *>
91 *> \param[in] N
92 *> \verbatim
93 *> N is INTEGER
94 *> The order of the matrices A, B, VL, and VR. N >= 0.
95 *> \endverbatim
96 *>
97 *> \param[in,out] A
98 *> \verbatim
99 *> A is DOUBLE PRECISION array, dimension (LDA, N)
100 *> On entry, the matrix A.
101 *> If JOBVL = 'V' or JOBVR = 'V', then on exit A
102 *> contains the real Schur form of A from the generalized Schur
103 *> factorization of the pair (A,B) after balancing.
104 *> If no eigenvectors were computed, then only the diagonal
105 *> blocks from the Schur form will be correct. See DGGHRD and
106 *> DHGEQZ for details.
107 *> \endverbatim
108 *>
109 *> \param[in] LDA
110 *> \verbatim
111 *> LDA is INTEGER
112 *> The leading dimension of A. LDA >= max(1,N).
113 *> \endverbatim
114 *>
115 *> \param[in,out] B
116 *> \verbatim
117 *> B is DOUBLE PRECISION array, dimension (LDB, N)
118 *> On entry, the matrix B.
119 *> If JOBVL = 'V' or JOBVR = 'V', then on exit B contains the
120 *> upper triangular matrix obtained from B in the generalized
121 *> Schur factorization of the pair (A,B) after balancing.
122 *> If no eigenvectors were computed, then only those elements of
123 *> B corresponding to the diagonal blocks from the Schur form of
124 *> A will be correct. See DGGHRD and DHGEQZ for details.
125 *> \endverbatim
126 *>
127 *> \param[in] LDB
128 *> \verbatim
129 *> LDB is INTEGER
130 *> The leading dimension of B. LDB >= max(1,N).
131 *> \endverbatim
132 *>
133 *> \param[out] ALPHAR
134 *> \verbatim
135 *> ALPHAR is DOUBLE PRECISION array, dimension (N)
136 *> The real parts of each scalar alpha defining an eigenvalue of
137 *> GNEP.
138 *> \endverbatim
139 *>
140 *> \param[out] ALPHAI
141 *> \verbatim
142 *> ALPHAI is DOUBLE PRECISION array, dimension (N)
143 *> The imaginary parts of each scalar alpha defining an
144 *> eigenvalue of GNEP. If ALPHAI(j) is zero, then the j-th
145 *> eigenvalue is real; if positive, then the j-th and
146 *> (j+1)-st eigenvalues are a complex conjugate pair, with
147 *> ALPHAI(j+1) = -ALPHAI(j).
148 *> \endverbatim
149 *>
150 *> \param[out] BETA
151 *> \verbatim
152 *> BETA is DOUBLE PRECISION array, dimension (N)
153 *> The scalars beta that define the eigenvalues of GNEP.
154 *>
155 *> Together, the quantities alpha = (ALPHAR(j),ALPHAI(j)) and
156 *> beta = BETA(j) represent the j-th eigenvalue of the matrix
157 *> pair (A,B), in one of the forms lambda = alpha/beta or
158 *> mu = beta/alpha. Since either lambda or mu may overflow,
159 *> they should not, in general, be computed.
160 *> \endverbatim
161 *>
162 *> \param[out] VL
163 *> \verbatim
164 *> VL is DOUBLE PRECISION array, dimension (LDVL,N)
165 *> If JOBVL = 'V', the left eigenvectors u(j) are stored
166 *> in the columns of VL, in the same order as their eigenvalues.
167 *> If the j-th eigenvalue is real, then u(j) = VL(:,j).
168 *> If the j-th and (j+1)-st eigenvalues form a complex conjugate
169 *> pair, then
170 *> u(j) = VL(:,j) + i*VL(:,j+1)
171 *> and
172 *> u(j+1) = VL(:,j) - i*VL(:,j+1).
173 *>
174 *> Each eigenvector is scaled so that its largest component has
175 *> abs(real part) + abs(imag. part) = 1, except for eigenvectors
176 *> corresponding to an eigenvalue with alpha = beta = 0, which
177 *> are set to zero.
178 *> Not referenced if JOBVL = 'N'.
179 *> \endverbatim
180 *>
181 *> \param[in] LDVL
182 *> \verbatim
183 *> LDVL is INTEGER
184 *> The leading dimension of the matrix VL. LDVL >= 1, and
185 *> if JOBVL = 'V', LDVL >= N.
186 *> \endverbatim
187 *>
188 *> \param[out] VR
189 *> \verbatim
190 *> VR is DOUBLE PRECISION array, dimension (LDVR,N)
191 *> If JOBVR = 'V', the right eigenvectors x(j) are stored
192 *> in the columns of VR, in the same order as their eigenvalues.
193 *> If the j-th eigenvalue is real, then x(j) = VR(:,j).
194 *> If the j-th and (j+1)-st eigenvalues form a complex conjugate
195 *> pair, then
196 *> x(j) = VR(:,j) + i*VR(:,j+1)
197 *> and
198 *> x(j+1) = VR(:,j) - i*VR(:,j+1).
199 *>
200 *> Each eigenvector is scaled so that its largest component has
201 *> abs(real part) + abs(imag. part) = 1, except for eigenvalues
202 *> corresponding to an eigenvalue with alpha = beta = 0, which
203 *> are set to zero.
204 *> Not referenced if JOBVR = 'N'.
205 *> \endverbatim
206 *>
207 *> \param[in] LDVR
208 *> \verbatim
209 *> LDVR is INTEGER
210 *> The leading dimension of the matrix VR. LDVR >= 1, and
211 *> if JOBVR = 'V', LDVR >= N.
212 *> \endverbatim
213 *>
214 *> \param[out] WORK
215 *> \verbatim
216 *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
217 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
218 *> \endverbatim
219 *>
220 *> \param[in] LWORK
221 *> \verbatim
222 *> LWORK is INTEGER
223 *> The dimension of the array WORK. LWORK >= max(1,8*N).
224 *> For good performance, LWORK must generally be larger.
225 *> To compute the optimal value of LWORK, call ILAENV to get
226 *> blocksizes (for DGEQRF, DORMQR, and DORGQR.) Then compute:
227 *> NB -- MAX of the blocksizes for DGEQRF, DORMQR, and DORGQR;
228 *> The optimal LWORK is:
229 *> 2*N + MAX( 6*N, N*(NB+1) ).
230 *>
231 *> If LWORK = -1, then a workspace query is assumed; the routine
232 *> only calculates the optimal size of the WORK array, returns
233 *> this value as the first entry of the WORK array, and no error
234 *> message related to LWORK is issued by XERBLA.
235 *> \endverbatim
236 *>
237 *> \param[out] INFO
238 *> \verbatim
239 *> INFO is INTEGER
240 *> = 0: successful exit
241 *> < 0: if INFO = -i, the i-th argument had an illegal value.
242 *> = 1,...,N:
243 *> The QZ iteration failed. No eigenvectors have been
244 *> calculated, but ALPHAR(j), ALPHAI(j), and BETA(j)
245 *> should be correct for j=INFO+1,...,N.
246 *> > N: errors that usually indicate LAPACK problems:
247 *> =N+1: error return from DGGBAL
248 *> =N+2: error return from DGEQRF
249 *> =N+3: error return from DORMQR
250 *> =N+4: error return from DORGQR
251 *> =N+5: error return from DGGHRD
252 *> =N+6: error return from DHGEQZ (other than failed
253 *> iteration)
254 *> =N+7: error return from DTGEVC
255 *> =N+8: error return from DGGBAK (computing VL)
256 *> =N+9: error return from DGGBAK (computing VR)
257 *> =N+10: error return from DLASCL (various calls)
258 *> \endverbatim
259 *
260 * Authors:
261 * ========
262 *
263 *> \author Univ. of Tennessee
264 *> \author Univ. of California Berkeley
265 *> \author Univ. of Colorado Denver
266 *> \author NAG Ltd.
267 *
268 *> \date November 2011
269 *
270 *> \ingroup doubleGEeigen
271 *
272 *> \par Further Details:
273 * =====================
274 *>
275 *> \verbatim
276 *>
277 *> Balancing
278 *> ---------
279 *>
280 *> This driver calls DGGBAL to both permute and scale rows and columns
281 *> of A and B. The permutations PL and PR are chosen so that PL*A*PR
282 *> and PL*B*R will be upper triangular except for the diagonal blocks
283 *> A(i:j,i:j) and B(i:j,i:j), with i and j as close together as
284 *> possible. The diagonal scaling matrices DL and DR are chosen so
285 *> that the pair DL*PL*A*PR*DR, DL*PL*B*PR*DR have elements close to
286 *> one (except for the elements that start out zero.)
287 *>
288 *> After the eigenvalues and eigenvectors of the balanced matrices
289 *> have been computed, DGGBAK transforms the eigenvectors back to what
290 *> they would have been (in perfect arithmetic) if they had not been
291 *> balanced.
292 *>
293 *> Contents of A and B on Exit
294 *> -------- -- - --- - -- ----
295 *>
296 *> If any eigenvectors are computed (either JOBVL='V' or JOBVR='V' or
297 *> both), then on exit the arrays A and B will contain the real Schur
298 *> form[*] of the "balanced" versions of A and B. If no eigenvectors
299 *> are computed, then only the diagonal blocks will be correct.
300 *>
301 *> [*] See DHGEQZ, DGEGS, or read the book "Matrix Computations",
302 *> by Golub & van Loan, pub. by Johns Hopkins U. Press.
303 *> \endverbatim
304 *>
305 * =====================================================================
306  SUBROUTINE dgegv( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR, ALPHAI,
307  $ beta, vl, ldvl, vr, ldvr, work, lwork, info )
308 *
309 * -- LAPACK driver routine (version 3.4.0) --
310 * -- LAPACK is a software package provided by Univ. of Tennessee, --
311 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
312 * November 2011
313 *
314 * .. Scalar Arguments ..
315  CHARACTER jobvl, jobvr
316  INTEGER info, lda, ldb, ldvl, ldvr, lwork, n
317 * ..
318 * .. Array Arguments ..
319  DOUBLE PRECISION a( lda, * ), alphai( * ), alphar( * ),
320  $ b( ldb, * ), beta( * ), vl( ldvl, * ),
321  $ vr( ldvr, * ), work( * )
322 * ..
323 *
324 * =====================================================================
325 *
326 * .. Parameters ..
327  DOUBLE PRECISION zero, one
328  parameter( zero = 0.0d0, one = 1.0d0 )
329 * ..
330 * .. Local Scalars ..
331  LOGICAL ilimit, ilv, ilvl, ilvr, lquery
332  CHARACTER chtemp
333  INTEGER icols, ihi, iinfo, ijobvl, ijobvr, ileft, ilo,
334  $ in, iright, irows, itau, iwork, jc, jr, lopt,
335  $ lwkmin, lwkopt, nb, nb1, nb2, nb3
336  DOUBLE PRECISION absai, absar, absb, anrm, anrm1, anrm2, bnrm,
337  $ bnrm1, bnrm2, eps, onepls, safmax, safmin,
338  $ salfai, salfar, sbeta, scale, temp
339 * ..
340 * .. Local Arrays ..
341  LOGICAL ldumma( 1 )
342 * ..
343 * .. External Subroutines ..
344  EXTERNAL dgeqrf, dggbak, dggbal, dgghrd, dhgeqz, dlacpy,
346 * ..
347 * .. External Functions ..
348  LOGICAL lsame
349  INTEGER ilaenv
350  DOUBLE PRECISION dlamch, dlange
351  EXTERNAL lsame, ilaenv, dlamch, dlange
352 * ..
353 * .. Intrinsic Functions ..
354  INTRINSIC abs, int, max
355 * ..
356 * .. Executable Statements ..
357 *
358 * Decode the input arguments
359 *
360  IF( lsame( jobvl, 'N' ) ) THEN
361  ijobvl = 1
362  ilvl = .false.
363  ELSE IF( lsame( jobvl, 'V' ) ) THEN
364  ijobvl = 2
365  ilvl = .true.
366  ELSE
367  ijobvl = -1
368  ilvl = .false.
369  END IF
370 *
371  IF( lsame( jobvr, 'N' ) ) THEN
372  ijobvr = 1
373  ilvr = .false.
374  ELSE IF( lsame( jobvr, 'V' ) ) THEN
375  ijobvr = 2
376  ilvr = .true.
377  ELSE
378  ijobvr = -1
379  ilvr = .false.
380  END IF
381  ilv = ilvl .OR. ilvr
382 *
383 * Test the input arguments
384 *
385  lwkmin = max( 8*n, 1 )
386  lwkopt = lwkmin
387  work( 1 ) = lwkopt
388  lquery = ( lwork.EQ.-1 )
389  info = 0
390  IF( ijobvl.LE.0 ) THEN
391  info = -1
392  ELSE IF( ijobvr.LE.0 ) THEN
393  info = -2
394  ELSE IF( n.LT.0 ) THEN
395  info = -3
396  ELSE IF( lda.LT.max( 1, n ) ) THEN
397  info = -5
398  ELSE IF( ldb.LT.max( 1, n ) ) THEN
399  info = -7
400  ELSE IF( ldvl.LT.1 .OR. ( ilvl .AND. ldvl.LT.n ) ) THEN
401  info = -12
402  ELSE IF( ldvr.LT.1 .OR. ( ilvr .AND. ldvr.LT.n ) ) THEN
403  info = -14
404  ELSE IF( lwork.LT.lwkmin .AND. .NOT.lquery ) THEN
405  info = -16
406  END IF
407 *
408  IF( info.EQ.0 ) THEN
409  nb1 = ilaenv( 1, 'DGEQRF', ' ', n, n, -1, -1 )
410  nb2 = ilaenv( 1, 'DORMQR', ' ', n, n, n, -1 )
411  nb3 = ilaenv( 1, 'DORGQR', ' ', n, n, n, -1 )
412  nb = max( nb1, nb2, nb3 )
413  lopt = 2*n + max( 6*n, n*( nb+1 ) )
414  work( 1 ) = lopt
415  END IF
416 *
417  IF( info.NE.0 ) THEN
418  CALL xerbla( 'DGEGV ', -info )
419  RETURN
420  ELSE IF( lquery ) THEN
421  RETURN
422  END IF
423 *
424 * Quick return if possible
425 *
426  IF( n.EQ.0 )
427  $ RETURN
428 *
429 * Get machine constants
430 *
431  eps = dlamch( 'E' )*dlamch( 'B' )
432  safmin = dlamch( 'S' )
433  safmin = safmin + safmin
434  safmax = one / safmin
435  onepls = one + ( 4*eps )
436 *
437 * Scale A
438 *
439  anrm = dlange( 'M', n, n, a, lda, work )
440  anrm1 = anrm
441  anrm2 = one
442  IF( anrm.LT.one ) THEN
443  IF( safmax*anrm.LT.one ) THEN
444  anrm1 = safmin
445  anrm2 = safmax*anrm
446  END IF
447  END IF
448 *
449  IF( anrm.GT.zero ) THEN
450  CALL dlascl( 'G', -1, -1, anrm, one, n, n, a, lda, iinfo )
451  IF( iinfo.NE.0 ) THEN
452  info = n + 10
453  RETURN
454  END IF
455  END IF
456 *
457 * Scale B
458 *
459  bnrm = dlange( 'M', n, n, b, ldb, work )
460  bnrm1 = bnrm
461  bnrm2 = one
462  IF( bnrm.LT.one ) THEN
463  IF( safmax*bnrm.LT.one ) THEN
464  bnrm1 = safmin
465  bnrm2 = safmax*bnrm
466  END IF
467  END IF
468 *
469  IF( bnrm.GT.zero ) THEN
470  CALL dlascl( 'G', -1, -1, bnrm, one, n, n, b, ldb, iinfo )
471  IF( iinfo.NE.0 ) THEN
472  info = n + 10
473  RETURN
474  END IF
475  END IF
476 *
477 * Permute the matrix to make it more nearly triangular
478 * Workspace layout: (8*N words -- "work" requires 6*N words)
479 * left_permutation, right_permutation, work...
480 *
481  ileft = 1
482  iright = n + 1
483  iwork = iright + n
484  CALL dggbal( 'P', n, a, lda, b, ldb, ilo, ihi, work( ileft ),
485  $ work( iright ), work( iwork ), iinfo )
486  IF( iinfo.NE.0 ) THEN
487  info = n + 1
488  go to 120
489  END IF
490 *
491 * Reduce B to triangular form, and initialize VL and/or VR
492 * Workspace layout: ("work..." must have at least N words)
493 * left_permutation, right_permutation, tau, work...
494 *
495  irows = ihi + 1 - ilo
496  IF( ilv ) THEN
497  icols = n + 1 - ilo
498  ELSE
499  icols = irows
500  END IF
501  itau = iwork
502  iwork = itau + irows
503  CALL dgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
504  $ work( iwork ), lwork+1-iwork, iinfo )
505  IF( iinfo.GE.0 )
506  $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
507  IF( iinfo.NE.0 ) THEN
508  info = n + 2
509  go to 120
510  END IF
511 *
512  CALL dormqr( 'L', 'T', irows, icols, irows, b( ilo, ilo ), ldb,
513  $ work( itau ), a( ilo, ilo ), lda, work( iwork ),
514  $ lwork+1-iwork, iinfo )
515  IF( iinfo.GE.0 )
516  $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
517  IF( iinfo.NE.0 ) THEN
518  info = n + 3
519  go to 120
520  END IF
521 *
522  IF( ilvl ) THEN
523  CALL dlaset( 'Full', n, n, zero, one, vl, ldvl )
524  CALL dlacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
525  $ vl( ilo+1, ilo ), ldvl )
526  CALL dorgqr( irows, irows, irows, vl( ilo, ilo ), ldvl,
527  $ work( itau ), work( iwork ), lwork+1-iwork,
528  $ iinfo )
529  IF( iinfo.GE.0 )
530  $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
531  IF( iinfo.NE.0 ) THEN
532  info = n + 4
533  go to 120
534  END IF
535  END IF
536 *
537  IF( ilvr )
538  $ CALL dlaset( 'Full', n, n, zero, one, vr, ldvr )
539 *
540 * Reduce to generalized Hessenberg form
541 *
542  IF( ilv ) THEN
543 *
544 * Eigenvectors requested -- work on whole matrix.
545 *
546  CALL dgghrd( jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb, vl,
547  $ ldvl, vr, ldvr, iinfo )
548  ELSE
549  CALL dgghrd( 'N', 'N', irows, 1, irows, a( ilo, ilo ), lda,
550  $ b( ilo, ilo ), ldb, vl, ldvl, vr, ldvr, iinfo )
551  END IF
552  IF( iinfo.NE.0 ) THEN
553  info = n + 5
554  go to 120
555  END IF
556 *
557 * Perform QZ algorithm
558 * Workspace layout: ("work..." must have at least 1 word)
559 * left_permutation, right_permutation, work...
560 *
561  iwork = itau
562  IF( ilv ) THEN
563  chtemp = 'S'
564  ELSE
565  chtemp = 'E'
566  END IF
567  CALL dhgeqz( chtemp, jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb,
568  $ alphar, alphai, beta, vl, ldvl, vr, ldvr,
569  $ work( iwork ), lwork+1-iwork, iinfo )
570  IF( iinfo.GE.0 )
571  $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
572  IF( iinfo.NE.0 ) THEN
573  IF( iinfo.GT.0 .AND. iinfo.LE.n ) THEN
574  info = iinfo
575  ELSE IF( iinfo.GT.n .AND. iinfo.LE.2*n ) THEN
576  info = iinfo - n
577  ELSE
578  info = n + 6
579  END IF
580  go to 120
581  END IF
582 *
583  IF( ilv ) THEN
584 *
585 * Compute Eigenvectors (DTGEVC requires 6*N words of workspace)
586 *
587  IF( ilvl ) THEN
588  IF( ilvr ) THEN
589  chtemp = 'B'
590  ELSE
591  chtemp = 'L'
592  END IF
593  ELSE
594  chtemp = 'R'
595  END IF
596 *
597  CALL dtgevc( chtemp, 'B', ldumma, n, a, lda, b, ldb, vl, ldvl,
598  $ vr, ldvr, n, in, work( iwork ), iinfo )
599  IF( iinfo.NE.0 ) THEN
600  info = n + 7
601  go to 120
602  END IF
603 *
604 * Undo balancing on VL and VR, rescale
605 *
606  IF( ilvl ) THEN
607  CALL dggbak( 'P', 'L', n, ilo, ihi, work( ileft ),
608  $ work( iright ), n, vl, ldvl, iinfo )
609  IF( iinfo.NE.0 ) THEN
610  info = n + 8
611  go to 120
612  END IF
613  DO 50 jc = 1, n
614  IF( alphai( jc ).LT.zero )
615  $ go to 50
616  temp = zero
617  IF( alphai( jc ).EQ.zero ) THEN
618  DO 10 jr = 1, n
619  temp = max( temp, abs( vl( jr, jc ) ) )
620  10 CONTINUE
621  ELSE
622  DO 20 jr = 1, n
623  temp = max( temp, abs( vl( jr, jc ) )+
624  $ abs( vl( jr, jc+1 ) ) )
625  20 CONTINUE
626  END IF
627  IF( temp.LT.safmin )
628  $ go to 50
629  temp = one / temp
630  IF( alphai( jc ).EQ.zero ) THEN
631  DO 30 jr = 1, n
632  vl( jr, jc ) = vl( jr, jc )*temp
633  30 CONTINUE
634  ELSE
635  DO 40 jr = 1, n
636  vl( jr, jc ) = vl( jr, jc )*temp
637  vl( jr, jc+1 ) = vl( jr, jc+1 )*temp
638  40 CONTINUE
639  END IF
640  50 CONTINUE
641  END IF
642  IF( ilvr ) THEN
643  CALL dggbak( 'P', 'R', n, ilo, ihi, work( ileft ),
644  $ work( iright ), n, vr, ldvr, iinfo )
645  IF( iinfo.NE.0 ) THEN
646  info = n + 9
647  go to 120
648  END IF
649  DO 100 jc = 1, n
650  IF( alphai( jc ).LT.zero )
651  $ go to 100
652  temp = zero
653  IF( alphai( jc ).EQ.zero ) THEN
654  DO 60 jr = 1, n
655  temp = max( temp, abs( vr( jr, jc ) ) )
656  60 CONTINUE
657  ELSE
658  DO 70 jr = 1, n
659  temp = max( temp, abs( vr( jr, jc ) )+
660  $ abs( vr( jr, jc+1 ) ) )
661  70 CONTINUE
662  END IF
663  IF( temp.LT.safmin )
664  $ go to 100
665  temp = one / temp
666  IF( alphai( jc ).EQ.zero ) THEN
667  DO 80 jr = 1, n
668  vr( jr, jc ) = vr( jr, jc )*temp
669  80 CONTINUE
670  ELSE
671  DO 90 jr = 1, n
672  vr( jr, jc ) = vr( jr, jc )*temp
673  vr( jr, jc+1 ) = vr( jr, jc+1 )*temp
674  90 CONTINUE
675  END IF
676  100 CONTINUE
677  END IF
678 *
679 * End of eigenvector calculation
680 *
681  END IF
682 *
683 * Undo scaling in alpha, beta
684 *
685 * Note: this does not give the alpha and beta for the unscaled
686 * problem.
687 *
688 * Un-scaling is limited to avoid underflow in alpha and beta
689 * if they are significant.
690 *
691  DO 110 jc = 1, n
692  absar = abs( alphar( jc ) )
693  absai = abs( alphai( jc ) )
694  absb = abs( beta( jc ) )
695  salfar = anrm*alphar( jc )
696  salfai = anrm*alphai( jc )
697  sbeta = bnrm*beta( jc )
698  ilimit = .false.
699  scale = one
700 *
701 * Check for significant underflow in ALPHAI
702 *
703  IF( abs( salfai ).LT.safmin .AND. absai.GE.
704  $ max( safmin, eps*absar, eps*absb ) ) THEN
705  ilimit = .true.
706  scale = ( onepls*safmin / anrm1 ) /
707  $ max( onepls*safmin, anrm2*absai )
708 *
709  ELSE IF( salfai.EQ.zero ) THEN
710 *
711 * If insignificant underflow in ALPHAI, then make the
712 * conjugate eigenvalue real.
713 *
714  IF( alphai( jc ).LT.zero .AND. jc.GT.1 ) THEN
715  alphai( jc-1 ) = zero
716  ELSE IF( alphai( jc ).GT.zero .AND. jc.LT.n ) THEN
717  alphai( jc+1 ) = zero
718  END IF
719  END IF
720 *
721 * Check for significant underflow in ALPHAR
722 *
723  IF( abs( salfar ).LT.safmin .AND. absar.GE.
724  $ max( safmin, eps*absai, eps*absb ) ) THEN
725  ilimit = .true.
726  scale = max( scale, ( onepls*safmin / anrm1 ) /
727  $ max( onepls*safmin, anrm2*absar ) )
728  END IF
729 *
730 * Check for significant underflow in BETA
731 *
732  IF( abs( sbeta ).LT.safmin .AND. absb.GE.
733  $ max( safmin, eps*absar, eps*absai ) ) THEN
734  ilimit = .true.
735  scale = max( scale, ( onepls*safmin / bnrm1 ) /
736  $ max( onepls*safmin, bnrm2*absb ) )
737  END IF
738 *
739 * Check for possible overflow when limiting scaling
740 *
741  IF( ilimit ) THEN
742  temp = ( scale*safmin )*max( abs( salfar ), abs( salfai ),
743  $ abs( sbeta ) )
744  IF( temp.GT.one )
745  $ scale = scale / temp
746  IF( scale.LT.one )
747  $ ilimit = .false.
748  END IF
749 *
750 * Recompute un-scaled ALPHAR, ALPHAI, BETA if necessary.
751 *
752  IF( ilimit ) THEN
753  salfar = ( scale*alphar( jc ) )*anrm
754  salfai = ( scale*alphai( jc ) )*anrm
755  sbeta = ( scale*beta( jc ) )*bnrm
756  END IF
757  alphar( jc ) = salfar
758  alphai( jc ) = salfai
759  beta( jc ) = sbeta
760  110 CONTINUE
761 *
762  120 CONTINUE
763  work( 1 ) = lwkopt
764 *
765  RETURN
766 *
767 * End of DGEGV
768 *
769  END