1*> \brief <b> DGGEV3 computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices (blocked algorithm)</b>
2*
3*  =========== DOCUMENTATION ===========
4*
5* Online html documentation available at
6*            http://www.netlib.org/lapack/explore-html/
7*
8*> \htmlonly
9*> Download DGGEV3 + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dggev3.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dggev3.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dggev3.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE DGGEV3( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR,
22*      $                   ALPHAI, BETA, VL, LDVL, VR, LDVR, WORK, LWORK,
23*      $                   INFO )
24*
25*       .. Scalar Arguments ..
26*       CHARACTER          JOBVL, JOBVR
27*       INTEGER            INFO, LDA, LDB, LDVL, LDVR, LWORK, N
28*       ..
29*       .. Array Arguments ..
30*       DOUBLE PRECISION   A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
31*      $                   B( LDB, * ), BETA( * ), VL( LDVL, * ),
32*      $                   VR( LDVR, * ), WORK( * )
33*       ..
34*
35*
36*> \par Purpose:
37*  =============
38*>
39*> \verbatim
40*>
41*> DGGEV3 computes for a pair of N-by-N real nonsymmetric matrices (A,B)
42*> the generalized eigenvalues, and optionally, the left and/or right
43*> generalized eigenvectors.
44*>
45*> A generalized eigenvalue for a pair of matrices (A,B) is a scalar
46*> lambda or a ratio alpha/beta = lambda, such that A - lambda*B is
47*> singular. It is usually represented as the pair (alpha,beta), as
48*> there is a reasonable interpretation for beta=0, and even for both
49*> being zero.
50*>
51*> The right eigenvector v(j) corresponding to the eigenvalue lambda(j)
52*> of (A,B) satisfies
53*>
54*>                  A * v(j) = lambda(j) * B * v(j).
55*>
56*> The left eigenvector u(j) corresponding to the eigenvalue lambda(j)
57*> of (A,B) satisfies
58*>
59*>                  u(j)**H * A  = lambda(j) * u(j)**H * B .
60*>
61*> where u(j)**H is the conjugate-transpose of u(j).
62*>
63*> \endverbatim
64*
65*  Arguments:
66*  ==========
67*
68*> \param[in] JOBVL
69*> \verbatim
70*>          JOBVL is CHARACTER*1
71*>          = 'N':  do not compute the left generalized eigenvectors;
72*>          = 'V':  compute the left generalized eigenvectors.
73*> \endverbatim
74*>
75*> \param[in] JOBVR
76*> \verbatim
77*>          JOBVR is CHARACTER*1
78*>          = 'N':  do not compute the right generalized eigenvectors;
79*>          = 'V':  compute the right generalized eigenvectors.
80*> \endverbatim
81*>
82*> \param[in] N
83*> \verbatim
84*>          N is INTEGER
85*>          The order of the matrices A, B, VL, and VR.  N >= 0.
86*> \endverbatim
87*>
88*> \param[in,out] A
89*> \verbatim
90*>          A is DOUBLE PRECISION array, dimension (LDA, N)
91*>          On entry, the matrix A in the pair (A,B).
92*>          On exit, A has been overwritten.
93*> \endverbatim
94*>
95*> \param[in] LDA
96*> \verbatim
97*>          LDA is INTEGER
98*>          The leading dimension of A.  LDA >= max(1,N).
99*> \endverbatim
100*>
101*> \param[in,out] B
102*> \verbatim
103*>          B is DOUBLE PRECISION array, dimension (LDB, N)
104*>          On entry, the matrix B in the pair (A,B).
105*>          On exit, B has been overwritten.
106*> \endverbatim
107*>
108*> \param[in] LDB
109*> \verbatim
110*>          LDB is INTEGER
111*>          The leading dimension of B.  LDB >= max(1,N).
112*> \endverbatim
113*>
114*> \param[out] ALPHAR
115*> \verbatim
116*>          ALPHAR is DOUBLE PRECISION array, dimension (N)
117*> \endverbatim
118*>
119*> \param[out] ALPHAI
120*> \verbatim
121*>          ALPHAI is DOUBLE PRECISION array, dimension (N)
122*> \endverbatim
123*>
124*> \param[out] BETA
125*> \verbatim
126*>          BETA is DOUBLE PRECISION array, dimension (N)
127*>          On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will
128*>          be the generalized eigenvalues.  If ALPHAI(j) is zero, then
129*>          the j-th eigenvalue is real; if positive, then the j-th and
130*>          (j+1)-st eigenvalues are a complex conjugate pair, with
131*>          ALPHAI(j+1) negative.
132*>
133*>          Note: the quotients ALPHAR(j)/BETA(j) and ALPHAI(j)/BETA(j)
134*>          may easily over- or underflow, and BETA(j) may even be zero.
135*>          Thus, the user should avoid naively computing the ratio
136*>          alpha/beta.  However, ALPHAR and ALPHAI will be always less
137*>          than and usually comparable with norm(A) in magnitude, and
138*>          BETA always less than and usually comparable with norm(B).
139*> \endverbatim
140*>
141*> \param[out] VL
142*> \verbatim
143*>          VL is DOUBLE PRECISION array, dimension (LDVL,N)
144*>          If JOBVL = 'V', the left eigenvectors u(j) are stored one
145*>          after another in the columns of VL, in the same order as
146*>          their eigenvalues. If the j-th eigenvalue is real, then
147*>          u(j) = VL(:,j), the j-th column of VL. If the j-th and
148*>          (j+1)-th eigenvalues form a complex conjugate pair, then
149*>          u(j) = VL(:,j)+i*VL(:,j+1) and u(j+1) = VL(:,j)-i*VL(:,j+1).
150*>          Each eigenvector is scaled so the largest component has
151*>          abs(real part)+abs(imag. part)=1.
152*>          Not referenced if JOBVL = 'N'.
153*> \endverbatim
154*>
155*> \param[in] LDVL
156*> \verbatim
157*>          LDVL is INTEGER
158*>          The leading dimension of the matrix VL. LDVL >= 1, and
159*>          if JOBVL = 'V', LDVL >= N.
160*> \endverbatim
161*>
162*> \param[out] VR
163*> \verbatim
164*>          VR is DOUBLE PRECISION array, dimension (LDVR,N)
165*>          If JOBVR = 'V', the right eigenvectors v(j) are stored one
166*>          after another in the columns of VR, in the same order as
167*>          their eigenvalues. If the j-th eigenvalue is real, then
168*>          v(j) = VR(:,j), the j-th column of VR. If the j-th and
169*>          (j+1)-th eigenvalues form a complex conjugate pair, then
170*>          v(j) = VR(:,j)+i*VR(:,j+1) and v(j+1) = VR(:,j)-i*VR(:,j+1).
171*>          Each eigenvector is scaled so the largest component has
172*>          abs(real part)+abs(imag. part)=1.
173*>          Not referenced if JOBVR = 'N'.
174*> \endverbatim
175*>
176*> \param[in] LDVR
177*> \verbatim
178*>          LDVR is INTEGER
179*>          The leading dimension of the matrix VR. LDVR >= 1, and
180*>          if JOBVR = 'V', LDVR >= N.
181*> \endverbatim
182*>
183*> \param[out] WORK
184*> \verbatim
185*>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
186*>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
187*> \endverbatim
188*>
189*> \param[in] LWORK
190*> \verbatim
191*>          LWORK is INTEGER
192*>
193*>          If LWORK = -1, then a workspace query is assumed; the routine
194*>          only calculates the optimal size of the WORK array, returns
195*>          this value as the first entry of the WORK array, and no error
196*>          message related to LWORK is issued by XERBLA.
197*> \endverbatim
198*>
199*> \param[out] INFO
200*> \verbatim
201*>          INFO is INTEGER
202*>          = 0:  successful exit
203*>          < 0:  if INFO = -i, the i-th argument had an illegal value.
204*>          = 1,...,N:
205*>                The QZ iteration failed.  No eigenvectors have been
206*>                calculated, but ALPHAR(j), ALPHAI(j), and BETA(j)
207*>                should be correct for j=INFO+1,...,N.
208*>          > N:  =N+1: other than QZ iteration failed in DLAQZ0.
209*>                =N+2: error return from DTGEVC.
210*> \endverbatim
211*
212*  Authors:
213*  ========
214*
215*> \author Univ. of Tennessee
216*> \author Univ. of California Berkeley
217*> \author Univ. of Colorado Denver
218*> \author NAG Ltd.
219*
220*> \ingroup doubleGEeigen
221*
222*  =====================================================================
223      SUBROUTINE DGGEV3( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR,
224     $                   ALPHAI, BETA, VL, LDVL, VR, LDVR, WORK, LWORK,
225     $                   INFO )
226*
227*  -- LAPACK driver routine --
228*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
229*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
230*
231*     .. Scalar Arguments ..
232      CHARACTER          JOBVL, JOBVR
233      INTEGER            INFO, LDA, LDB, LDVL, LDVR, LWORK, N
234*     ..
235*     .. Array Arguments ..
236      DOUBLE PRECISION   A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
237     $                   B( LDB, * ), BETA( * ), VL( LDVL, * ),
238     $                   VR( LDVR, * ), WORK( * )
239*     ..
240*
241*  =====================================================================
242*
243*     .. Parameters ..
244      DOUBLE PRECISION   ZERO, ONE
245      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0 )
246*     ..
247*     .. Local Scalars ..
248      LOGICAL            ILASCL, ILBSCL, ILV, ILVL, ILVR, LQUERY
249      CHARACTER          CHTEMP
250      INTEGER            ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT, ILO,
251     $                   IN, IRIGHT, IROWS, ITAU, IWRK, JC, JR, LWKOPT
252      DOUBLE PRECISION   ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
253     $                   SMLNUM, TEMP
254*     ..
255*     .. Local Arrays ..
256      LOGICAL            LDUMMA( 1 )
257*     ..
258*     .. External Subroutines ..
259      EXTERNAL           DGEQRF, DGGBAK, DGGBAL, DGGHD3, DLAQZ0, DLABAD,
260     $                   DLACPY, DLASCL, DLASET, DORGQR, DORMQR, DTGEVC,
261     $                   XERBLA
262*     ..
263*     .. External Functions ..
264      LOGICAL            LSAME
265      DOUBLE PRECISION   DLAMCH, DLANGE
266      EXTERNAL           LSAME, DLAMCH, DLANGE
267*     ..
268*     .. Intrinsic Functions ..
269      INTRINSIC          ABS, MAX, SQRT
270*     ..
271*     .. Executable Statements ..
272*
273*     Decode the input arguments
274*
275      IF( LSAME( JOBVL, 'N' ) ) THEN
276         IJOBVL = 1
277         ILVL = .FALSE.
278      ELSE IF( LSAME( JOBVL, 'V' ) ) THEN
279         IJOBVL = 2
280         ILVL = .TRUE.
281      ELSE
282         IJOBVL = -1
283         ILVL = .FALSE.
284      END IF
285*
286      IF( LSAME( JOBVR, 'N' ) ) THEN
287         IJOBVR = 1
288         ILVR = .FALSE.
289      ELSE IF( LSAME( JOBVR, 'V' ) ) THEN
290         IJOBVR = 2
291         ILVR = .TRUE.
292      ELSE
293         IJOBVR = -1
294         ILVR = .FALSE.
295      END IF
296      ILV = ILVL .OR. ILVR
297*
298*     Test the input arguments
299*
300      INFO = 0
301      LQUERY = ( LWORK.EQ.-1 )
302      IF( IJOBVL.LE.0 ) THEN
303         INFO = -1
304      ELSE IF( IJOBVR.LE.0 ) THEN
305         INFO = -2
306      ELSE IF( N.LT.0 ) THEN
307         INFO = -3
308      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
309         INFO = -5
310      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
311         INFO = -7
312      ELSE IF( LDVL.LT.1 .OR. ( ILVL .AND. LDVL.LT.N ) ) THEN
313         INFO = -12
314      ELSE IF( LDVR.LT.1 .OR. ( ILVR .AND. LDVR.LT.N ) ) THEN
315         INFO = -14
316      ELSE IF( LWORK.LT.MAX( 1, 8*N ) .AND. .NOT.LQUERY ) THEN
317         INFO = -16
318      END IF
319*
320*     Compute workspace
321*
322      IF( INFO.EQ.0 ) THEN
323         CALL DGEQRF( N, N, B, LDB, WORK, WORK, -1, IERR )
324         LWKOPT = MAX(1, 8*N, 3*N+INT( WORK( 1 ) ) )
325         CALL DORMQR( 'L', 'T', N, N, N, B, LDB, WORK, A, LDA, WORK, -1,
326     $                IERR )
327         LWKOPT = MAX( LWKOPT, 3*N+INT( WORK ( 1 ) ) )
328         IF( ILVL ) THEN
329            CALL DORGQR( N, N, N, VL, LDVL, WORK, WORK, -1, IERR )
330            LWKOPT = MAX( LWKOPT, 3*N+INT( WORK ( 1 ) ) )
331         END IF
332         IF( ILV ) THEN
333            CALL DGGHD3( JOBVL, JOBVR, N, 1, N, A, LDA, B, LDB, VL,
334     $                   LDVL, VR, LDVR, WORK, -1, IERR )
335            LWKOPT = MAX( LWKOPT, 3*N+INT( WORK ( 1 ) ) )
336            CALL DLAQZ0( 'S', JOBVL, JOBVR, N, 1, N, A, LDA, B, LDB,
337     $                   ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR,
338     $                   WORK, -1, 0, IERR )
339            LWKOPT = MAX( LWKOPT, 2*N+INT( WORK ( 1 ) ) )
340         ELSE
341            CALL DGGHD3( 'N', 'N', N, 1, N, A, LDA, B, LDB, VL, LDVL,
342     $                   VR, LDVR, WORK, -1, IERR )
343            LWKOPT = MAX( LWKOPT, 3*N+INT( WORK ( 1 ) ) )
344            CALL DLAQZ0( 'E', JOBVL, JOBVR, N, 1, N, A, LDA, B, LDB,
345     $                   ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR,
346     $                   WORK, -1, 0, IERR )
347            LWKOPT = MAX( LWKOPT, 2*N+INT( WORK ( 1 ) ) )
348         END IF
349
350         WORK( 1 ) = LWKOPT
351      END IF
352*
353      IF( INFO.NE.0 ) THEN
354         CALL XERBLA( 'DGGEV3 ', -INFO )
355         RETURN
356      ELSE IF( LQUERY ) THEN
357         RETURN
358      END IF
359*
360*     Quick return if possible
361*
362      IF( N.EQ.0 )
363     $   RETURN
364*
365*     Get machine constants
366*
367      EPS = DLAMCH( 'P' )
368      SMLNUM = DLAMCH( 'S' )
369      BIGNUM = ONE / SMLNUM
370      CALL DLABAD( SMLNUM, BIGNUM )
371      SMLNUM = SQRT( SMLNUM ) / EPS
372      BIGNUM = ONE / SMLNUM
373*
374*     Scale A if max element outside range [SMLNUM,BIGNUM]
375*
376      ANRM = DLANGE( 'M', N, N, A, LDA, WORK )
377      ILASCL = .FALSE.
378      IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
379         ANRMTO = SMLNUM
380         ILASCL = .TRUE.
381      ELSE IF( ANRM.GT.BIGNUM ) THEN
382         ANRMTO = BIGNUM
383         ILASCL = .TRUE.
384      END IF
385      IF( ILASCL )
386     $   CALL DLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
387*
388*     Scale B if max element outside range [SMLNUM,BIGNUM]
389*
390      BNRM = DLANGE( 'M', N, N, B, LDB, WORK )
391      ILBSCL = .FALSE.
392      IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
393         BNRMTO = SMLNUM
394         ILBSCL = .TRUE.
395      ELSE IF( BNRM.GT.BIGNUM ) THEN
396         BNRMTO = BIGNUM
397         ILBSCL = .TRUE.
398      END IF
399      IF( ILBSCL )
400     $   CALL DLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
401*
402*     Permute the matrices A, B to isolate eigenvalues if possible
403*
404      ILEFT = 1
405      IRIGHT = N + 1
406      IWRK = IRIGHT + N
407      CALL DGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, WORK( ILEFT ),
408     $             WORK( IRIGHT ), WORK( IWRK ), IERR )
409*
410*     Reduce B to triangular form (QR decomposition of B)
411*
412      IROWS = IHI + 1 - ILO
413      IF( ILV ) THEN
414         ICOLS = N + 1 - ILO
415      ELSE
416         ICOLS = IROWS
417      END IF
418      ITAU = IWRK
419      IWRK = ITAU + IROWS
420      CALL DGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
421     $             WORK( IWRK ), LWORK+1-IWRK, IERR )
422*
423*     Apply the orthogonal transformation to matrix A
424*
425      CALL DORMQR( 'L', 'T', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
426     $             WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ),
427     $             LWORK+1-IWRK, IERR )
428*
429*     Initialize VL
430*
431      IF( ILVL ) THEN
432         CALL DLASET( 'Full', N, N, ZERO, ONE, VL, LDVL )
433         IF( IROWS.GT.1 ) THEN
434            CALL DLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
435     $                   VL( ILO+1, ILO ), LDVL )
436         END IF
437         CALL DORGQR( IROWS, IROWS, IROWS, VL( ILO, ILO ), LDVL,
438     $                WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
439      END IF
440*
441*     Initialize VR
442*
443      IF( ILVR )
444     $   CALL DLASET( 'Full', N, N, ZERO, ONE, VR, LDVR )
445*
446*     Reduce to generalized Hessenberg form
447*
448      IF( ILV ) THEN
449*
450*        Eigenvectors requested -- work on whole matrix.
451*
452         CALL DGGHD3( JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, VL,
453     $                LDVL, VR, LDVR, WORK( IWRK ), LWORK+1-IWRK, IERR )
454      ELSE
455         CALL DGGHD3( 'N', 'N', IROWS, 1, IROWS, A( ILO, ILO ), LDA,
456     $                B( ILO, ILO ), LDB, VL, LDVL, VR, LDVR,
457     $                WORK( IWRK ), LWORK+1-IWRK, IERR )
458      END IF
459*
460*     Perform QZ algorithm (Compute eigenvalues, and optionally, the
461*     Schur forms and Schur vectors)
462*
463      IWRK = ITAU
464      IF( ILV ) THEN
465         CHTEMP = 'S'
466      ELSE
467         CHTEMP = 'E'
468      END IF
469      CALL DLAQZ0( CHTEMP, JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB,
470     $             ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR,
471     $             WORK( IWRK ), LWORK+1-IWRK, 0, IERR )
472      IF( IERR.NE.0 ) THEN
473         IF( IERR.GT.0 .AND. IERR.LE.N ) THEN
474            INFO = IERR
475         ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN
476            INFO = IERR - N
477         ELSE
478            INFO = N + 1
479         END IF
480         GO TO 110
481      END IF
482*
483*     Compute Eigenvectors
484*
485      IF( ILV ) THEN
486         IF( ILVL ) THEN
487            IF( ILVR ) THEN
488               CHTEMP = 'B'
489            ELSE
490               CHTEMP = 'L'
491            END IF
492         ELSE
493            CHTEMP = 'R'
494         END IF
495         CALL DTGEVC( CHTEMP, 'B', LDUMMA, N, A, LDA, B, LDB, VL, LDVL,
496     $                VR, LDVR, N, IN, WORK( IWRK ), IERR )
497         IF( IERR.NE.0 ) THEN
498            INFO = N + 2
499            GO TO 110
500         END IF
501*
502*        Undo balancing on VL and VR and normalization
503*
504         IF( ILVL ) THEN
505            CALL DGGBAK( 'P', 'L', N, ILO, IHI, WORK( ILEFT ),
506     $                   WORK( IRIGHT ), N, VL, LDVL, IERR )
507            DO 50 JC = 1, N
508               IF( ALPHAI( JC ).LT.ZERO )
509     $            GO TO 50
510               TEMP = ZERO
511               IF( ALPHAI( JC ).EQ.ZERO ) THEN
512                  DO 10 JR = 1, N
513                     TEMP = MAX( TEMP, ABS( VL( JR, JC ) ) )
514   10             CONTINUE
515               ELSE
516                  DO 20 JR = 1, N
517                     TEMP = MAX( TEMP, ABS( VL( JR, JC ) )+
518     $                      ABS( VL( JR, JC+1 ) ) )
519   20             CONTINUE
520               END IF
521               IF( TEMP.LT.SMLNUM )
522     $            GO TO 50
523               TEMP = ONE / TEMP
524               IF( ALPHAI( JC ).EQ.ZERO ) THEN
525                  DO 30 JR = 1, N
526                     VL( JR, JC ) = VL( JR, JC )*TEMP
527   30             CONTINUE
528               ELSE
529                  DO 40 JR = 1, N
530                     VL( JR, JC ) = VL( JR, JC )*TEMP
531                     VL( JR, JC+1 ) = VL( JR, JC+1 )*TEMP
532   40             CONTINUE
533               END IF
534   50       CONTINUE
535         END IF
536         IF( ILVR ) THEN
537            CALL DGGBAK( 'P', 'R', N, ILO, IHI, WORK( ILEFT ),
538     $                   WORK( IRIGHT ), N, VR, LDVR, IERR )
539            DO 100 JC = 1, N
540               IF( ALPHAI( JC ).LT.ZERO )
541     $            GO TO 100
542               TEMP = ZERO
543               IF( ALPHAI( JC ).EQ.ZERO ) THEN
544                  DO 60 JR = 1, N
545                     TEMP = MAX( TEMP, ABS( VR( JR, JC ) ) )
546   60             CONTINUE
547               ELSE
548                  DO 70 JR = 1, N
549                     TEMP = MAX( TEMP, ABS( VR( JR, JC ) )+
550     $                      ABS( VR( JR, JC+1 ) ) )
551   70             CONTINUE
552               END IF
553               IF( TEMP.LT.SMLNUM )
554     $            GO TO 100
555               TEMP = ONE / TEMP
556               IF( ALPHAI( JC ).EQ.ZERO ) THEN
557                  DO 80 JR = 1, N
558                     VR( JR, JC ) = VR( JR, JC )*TEMP
559   80             CONTINUE
560               ELSE
561                  DO 90 JR = 1, N
562                     VR( JR, JC ) = VR( JR, JC )*TEMP
563                     VR( JR, JC+1 ) = VR( JR, JC+1 )*TEMP
564   90             CONTINUE
565               END IF
566  100       CONTINUE
567         END IF
568*
569*        End of eigenvector calculation
570*
571      END IF
572*
573*     Undo scaling if necessary
574*
575  110 CONTINUE
576*
577      IF( ILASCL ) THEN
578         CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAR, N, IERR )
579         CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAI, N, IERR )
580      END IF
581*
582      IF( ILBSCL ) THEN
583         CALL DLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
584      END IF
585*
586      WORK( 1 ) = LWKOPT
587      RETURN
588*
589*     End of DGGEV3
590*
591      END
592