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