LAPACK  3.10.0
LAPACK: Linear Algebra PACKage
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 *> \ingroup doubleGEeigen
269 *
270 *> \par Further Details:
271 * =====================
272 *>
273 *> \verbatim
274 *>
275 *> Balancing
276 *> ---------
277 *>
278 *> This driver calls DGGBAL to both permute and scale rows and columns
279 *> of A and B. The permutations PL and PR are chosen so that PL*A*PR
280 *> and PL*B*R will be upper triangular except for the diagonal blocks
281 *> A(i:j,i:j) and B(i:j,i:j), with i and j as close together as
282 *> possible. The diagonal scaling matrices DL and DR are chosen so
283 *> that the pair DL*PL*A*PR*DR, DL*PL*B*PR*DR have elements close to
284 *> one (except for the elements that start out zero.)
285 *>
286 *> After the eigenvalues and eigenvectors of the balanced matrices
287 *> have been computed, DGGBAK transforms the eigenvectors back to what
288 *> they would have been (in perfect arithmetic) if they had not been
289 *> balanced.
290 *>
291 *> Contents of A and B on Exit
292 *> -------- -- - --- - -- ----
293 *>
294 *> If any eigenvectors are computed (either JOBVL='V' or JOBVR='V' or
295 *> both), then on exit the arrays A and B will contain the real Schur
296 *> form[*] of the "balanced" versions of A and B. If no eigenvectors
297 *> are computed, then only the diagonal blocks will be correct.
298 *>
299 *> [*] See DHGEQZ, DGEGS, or read the book "Matrix Computations",
300 *> by Golub & van Loan, pub. by Johns Hopkins U. Press.
301 *> \endverbatim
302 *>
303 * =====================================================================
304  SUBROUTINE dgegv( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR, ALPHAI,
305  $ BETA, VL, LDVL, VR, LDVR, WORK, LWORK, INFO )
306 *
307 * -- LAPACK driver routine --
308 * -- LAPACK is a software package provided by Univ. of Tennessee, --
309 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
310 *
311 * .. Scalar Arguments ..
312  CHARACTER JOBVL, JOBVR
313  INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, N
314 * ..
315 * .. Array Arguments ..
316  DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
317  $ b( ldb, * ), beta( * ), vl( ldvl, * ),
318  $ vr( ldvr, * ), work( * )
319 * ..
320 *
321 * =====================================================================
322 *
323 * .. Parameters ..
324  DOUBLE PRECISION ZERO, ONE
325  parameter( zero = 0.0d0, one = 1.0d0 )
326 * ..
327 * .. Local Scalars ..
328  LOGICAL ILIMIT, ILV, ILVL, ILVR, LQUERY
329  CHARACTER CHTEMP
330  INTEGER ICOLS, IHI, IINFO, IJOBVL, IJOBVR, ILEFT, ILO,
331  $ in, iright, irows, itau, iwork, jc, jr, lopt,
332  $ lwkmin, lwkopt, nb, nb1, nb2, nb3
333  DOUBLE PRECISION ABSAI, ABSAR, ABSB, ANRM, ANRM1, ANRM2, BNRM,
334  $ bnrm1, bnrm2, eps, onepls, safmax, safmin,
335  $ salfai, salfar, sbeta, scale, temp
336 * ..
337 * .. Local Arrays ..
338  LOGICAL LDUMMA( 1 )
339 * ..
340 * .. External Subroutines ..
341  EXTERNAL dgeqrf, dggbak, dggbal, dgghrd, dhgeqz, dlacpy,
343 * ..
344 * .. External Functions ..
345  LOGICAL LSAME
346  INTEGER ILAENV
347  DOUBLE PRECISION DLAMCH, DLANGE
348  EXTERNAL lsame, ilaenv, dlamch, dlange
349 * ..
350 * .. Intrinsic Functions ..
351  INTRINSIC abs, int, max
352 * ..
353 * .. Executable Statements ..
354 *
355 * Decode the input arguments
356 *
357  IF( lsame( jobvl, 'N' ) ) THEN
358  ijobvl = 1
359  ilvl = .false.
360  ELSE IF( lsame( jobvl, 'V' ) ) THEN
361  ijobvl = 2
362  ilvl = .true.
363  ELSE
364  ijobvl = -1
365  ilvl = .false.
366  END IF
367 *
368  IF( lsame( jobvr, 'N' ) ) THEN
369  ijobvr = 1
370  ilvr = .false.
371  ELSE IF( lsame( jobvr, 'V' ) ) THEN
372  ijobvr = 2
373  ilvr = .true.
374  ELSE
375  ijobvr = -1
376  ilvr = .false.
377  END IF
378  ilv = ilvl .OR. ilvr
379 *
380 * Test the input arguments
381 *
382  lwkmin = max( 8*n, 1 )
383  lwkopt = lwkmin
384  work( 1 ) = lwkopt
385  lquery = ( lwork.EQ.-1 )
386  info = 0
387  IF( ijobvl.LE.0 ) THEN
388  info = -1
389  ELSE IF( ijobvr.LE.0 ) THEN
390  info = -2
391  ELSE IF( n.LT.0 ) THEN
392  info = -3
393  ELSE IF( lda.LT.max( 1, n ) ) THEN
394  info = -5
395  ELSE IF( ldb.LT.max( 1, n ) ) THEN
396  info = -7
397  ELSE IF( ldvl.LT.1 .OR. ( ilvl .AND. ldvl.LT.n ) ) THEN
398  info = -12
399  ELSE IF( ldvr.LT.1 .OR. ( ilvr .AND. ldvr.LT.n ) ) THEN
400  info = -14
401  ELSE IF( lwork.LT.lwkmin .AND. .NOT.lquery ) THEN
402  info = -16
403  END IF
404 *
405  IF( info.EQ.0 ) THEN
406  nb1 = ilaenv( 1, 'DGEQRF', ' ', n, n, -1, -1 )
407  nb2 = ilaenv( 1, 'DORMQR', ' ', n, n, n, -1 )
408  nb3 = ilaenv( 1, 'DORGQR', ' ', n, n, n, -1 )
409  nb = max( nb1, nb2, nb3 )
410  lopt = 2*n + max( 6*n, n*( nb+1 ) )
411  work( 1 ) = lopt
412  END IF
413 *
414  IF( info.NE.0 ) THEN
415  CALL xerbla( 'DGEGV ', -info )
416  RETURN
417  ELSE IF( lquery ) THEN
418  RETURN
419  END IF
420 *
421 * Quick return if possible
422 *
423  IF( n.EQ.0 )
424  $ RETURN
425 *
426 * Get machine constants
427 *
428  eps = dlamch( 'E' )*dlamch( 'B' )
429  safmin = dlamch( 'S' )
430  safmin = safmin + safmin
431  safmax = one / safmin
432  onepls = one + ( 4*eps )
433 *
434 * Scale A
435 *
436  anrm = dlange( 'M', n, n, a, lda, work )
437  anrm1 = anrm
438  anrm2 = one
439  IF( anrm.LT.one ) THEN
440  IF( safmax*anrm.LT.one ) THEN
441  anrm1 = safmin
442  anrm2 = safmax*anrm
443  END IF
444  END IF
445 *
446  IF( anrm.GT.zero ) THEN
447  CALL dlascl( 'G', -1, -1, anrm, one, n, n, a, lda, iinfo )
448  IF( iinfo.NE.0 ) THEN
449  info = n + 10
450  RETURN
451  END IF
452  END IF
453 *
454 * Scale B
455 *
456  bnrm = dlange( 'M', n, n, b, ldb, work )
457  bnrm1 = bnrm
458  bnrm2 = one
459  IF( bnrm.LT.one ) THEN
460  IF( safmax*bnrm.LT.one ) THEN
461  bnrm1 = safmin
462  bnrm2 = safmax*bnrm
463  END IF
464  END IF
465 *
466  IF( bnrm.GT.zero ) THEN
467  CALL dlascl( 'G', -1, -1, bnrm, one, n, n, b, ldb, iinfo )
468  IF( iinfo.NE.0 ) THEN
469  info = n + 10
470  RETURN
471  END IF
472  END IF
473 *
474 * Permute the matrix to make it more nearly triangular
475 * Workspace layout: (8*N words -- "work" requires 6*N words)
476 * left_permutation, right_permutation, work...
477 *
478  ileft = 1
479  iright = n + 1
480  iwork = iright + n
481  CALL dggbal( 'P', n, a, lda, b, ldb, ilo, ihi, work( ileft ),
482  $ work( iright ), work( iwork ), iinfo )
483  IF( iinfo.NE.0 ) THEN
484  info = n + 1
485  GO TO 120
486  END IF
487 *
488 * Reduce B to triangular form, and initialize VL and/or VR
489 * Workspace layout: ("work..." must have at least N words)
490 * left_permutation, right_permutation, tau, work...
491 *
492  irows = ihi + 1 - ilo
493  IF( ilv ) THEN
494  icols = n + 1 - ilo
495  ELSE
496  icols = irows
497  END IF
498  itau = iwork
499  iwork = itau + irows
500  CALL dgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
501  $ work( iwork ), lwork+1-iwork, iinfo )
502  IF( iinfo.GE.0 )
503  $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
504  IF( iinfo.NE.0 ) THEN
505  info = n + 2
506  GO TO 120
507  END IF
508 *
509  CALL dormqr( 'L', 'T', irows, icols, irows, b( ilo, ilo ), ldb,
510  $ work( itau ), a( ilo, ilo ), lda, work( iwork ),
511  $ lwork+1-iwork, iinfo )
512  IF( iinfo.GE.0 )
513  $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
514  IF( iinfo.NE.0 ) THEN
515  info = n + 3
516  GO TO 120
517  END IF
518 *
519  IF( ilvl ) THEN
520  CALL dlaset( 'Full', n, n, zero, one, vl, ldvl )
521  CALL dlacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
522  $ vl( ilo+1, ilo ), ldvl )
523  CALL dorgqr( irows, irows, irows, vl( ilo, ilo ), ldvl,
524  $ work( itau ), work( iwork ), lwork+1-iwork,
525  $ iinfo )
526  IF( iinfo.GE.0 )
527  $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
528  IF( iinfo.NE.0 ) THEN
529  info = n + 4
530  GO TO 120
531  END IF
532  END IF
533 *
534  IF( ilvr )
535  $ CALL dlaset( 'Full', n, n, zero, one, vr, ldvr )
536 *
537 * Reduce to generalized Hessenberg form
538 *
539  IF( ilv ) THEN
540 *
541 * Eigenvectors requested -- work on whole matrix.
542 *
543  CALL dgghrd( jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb, vl,
544  $ ldvl, vr, ldvr, iinfo )
545  ELSE
546  CALL dgghrd( 'N', 'N', irows, 1, irows, a( ilo, ilo ), lda,
547  $ b( ilo, ilo ), ldb, vl, ldvl, vr, ldvr, iinfo )
548  END IF
549  IF( iinfo.NE.0 ) THEN
550  info = n + 5
551  GO TO 120
552  END IF
553 *
554 * Perform QZ algorithm
555 * Workspace layout: ("work..." must have at least 1 word)
556 * left_permutation, right_permutation, work...
557 *
558  iwork = itau
559  IF( ilv ) THEN
560  chtemp = 'S'
561  ELSE
562  chtemp = 'E'
563  END IF
564  CALL dhgeqz( chtemp, jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb,
565  $ alphar, alphai, beta, vl, ldvl, vr, ldvr,
566  $ work( iwork ), lwork+1-iwork, iinfo )
567  IF( iinfo.GE.0 )
568  $ lwkopt = max( lwkopt, int( work( iwork ) )+iwork-1 )
569  IF( iinfo.NE.0 ) THEN
570  IF( iinfo.GT.0 .AND. iinfo.LE.n ) THEN
571  info = iinfo
572  ELSE IF( iinfo.GT.n .AND. iinfo.LE.2*n ) THEN
573  info = iinfo - n
574  ELSE
575  info = n + 6
576  END IF
577  GO TO 120
578  END IF
579 *
580  IF( ilv ) THEN
581 *
582 * Compute Eigenvectors (DTGEVC requires 6*N words of workspace)
583 *
584  IF( ilvl ) THEN
585  IF( ilvr ) THEN
586  chtemp = 'B'
587  ELSE
588  chtemp = 'L'
589  END IF
590  ELSE
591  chtemp = 'R'
592  END IF
593 *
594  CALL dtgevc( chtemp, 'B', ldumma, n, a, lda, b, ldb, vl, ldvl,
595  $ vr, ldvr, n, in, work( iwork ), iinfo )
596  IF( iinfo.NE.0 ) THEN
597  info = n + 7
598  GO TO 120
599  END IF
600 *
601 * Undo balancing on VL and VR, rescale
602 *
603  IF( ilvl ) THEN
604  CALL dggbak( 'P', 'L', n, ilo, ihi, work( ileft ),
605  $ work( iright ), n, vl, ldvl, iinfo )
606  IF( iinfo.NE.0 ) THEN
607  info = n + 8
608  GO TO 120
609  END IF
610  DO 50 jc = 1, n
611  IF( alphai( jc ).LT.zero )
612  $ GO TO 50
613  temp = zero
614  IF( alphai( jc ).EQ.zero ) THEN
615  DO 10 jr = 1, n
616  temp = max( temp, abs( vl( jr, jc ) ) )
617  10 CONTINUE
618  ELSE
619  DO 20 jr = 1, n
620  temp = max( temp, abs( vl( jr, jc ) )+
621  $ abs( vl( jr, jc+1 ) ) )
622  20 CONTINUE
623  END IF
624  IF( temp.LT.safmin )
625  $ GO TO 50
626  temp = one / temp
627  IF( alphai( jc ).EQ.zero ) THEN
628  DO 30 jr = 1, n
629  vl( jr, jc ) = vl( jr, jc )*temp
630  30 CONTINUE
631  ELSE
632  DO 40 jr = 1, n
633  vl( jr, jc ) = vl( jr, jc )*temp
634  vl( jr, jc+1 ) = vl( jr, jc+1 )*temp
635  40 CONTINUE
636  END IF
637  50 CONTINUE
638  END IF
639  IF( ilvr ) THEN
640  CALL dggbak( 'P', 'R', n, ilo, ihi, work( ileft ),
641  $ work( iright ), n, vr, ldvr, iinfo )
642  IF( iinfo.NE.0 ) THEN
643  info = n + 9
644  GO TO 120
645  END IF
646  DO 100 jc = 1, n
647  IF( alphai( jc ).LT.zero )
648  $ GO TO 100
649  temp = zero
650  IF( alphai( jc ).EQ.zero ) THEN
651  DO 60 jr = 1, n
652  temp = max( temp, abs( vr( jr, jc ) ) )
653  60 CONTINUE
654  ELSE
655  DO 70 jr = 1, n
656  temp = max( temp, abs( vr( jr, jc ) )+
657  $ abs( vr( jr, jc+1 ) ) )
658  70 CONTINUE
659  END IF
660  IF( temp.LT.safmin )
661  $ GO TO 100
662  temp = one / temp
663  IF( alphai( jc ).EQ.zero ) THEN
664  DO 80 jr = 1, n
665  vr( jr, jc ) = vr( jr, jc )*temp
666  80 CONTINUE
667  ELSE
668  DO 90 jr = 1, n
669  vr( jr, jc ) = vr( jr, jc )*temp
670  vr( jr, jc+1 ) = vr( jr, jc+1 )*temp
671  90 CONTINUE
672  END IF
673  100 CONTINUE
674  END IF
675 *
676 * End of eigenvector calculation
677 *
678  END IF
679 *
680 * Undo scaling in alpha, beta
681 *
682 * Note: this does not give the alpha and beta for the unscaled
683 * problem.
684 *
685 * Un-scaling is limited to avoid underflow in alpha and beta
686 * if they are significant.
687 *
688  DO 110 jc = 1, n
689  absar = abs( alphar( jc ) )
690  absai = abs( alphai( jc ) )
691  absb = abs( beta( jc ) )
692  salfar = anrm*alphar( jc )
693  salfai = anrm*alphai( jc )
694  sbeta = bnrm*beta( jc )
695  ilimit = .false.
696  scale = one
697 *
698 * Check for significant underflow in ALPHAI
699 *
700  IF( abs( salfai ).LT.safmin .AND. absai.GE.
701  $ max( safmin, eps*absar, eps*absb ) ) THEN
702  ilimit = .true.
703  scale = ( onepls*safmin / anrm1 ) /
704  $ max( onepls*safmin, anrm2*absai )
705 *
706  ELSE IF( salfai.EQ.zero ) THEN
707 *
708 * If insignificant underflow in ALPHAI, then make the
709 * conjugate eigenvalue real.
710 *
711  IF( alphai( jc ).LT.zero .AND. jc.GT.1 ) THEN
712  alphai( jc-1 ) = zero
713  ELSE IF( alphai( jc ).GT.zero .AND. jc.LT.n ) THEN
714  alphai( jc+1 ) = zero
715  END IF
716  END IF
717 *
718 * Check for significant underflow in ALPHAR
719 *
720  IF( abs( salfar ).LT.safmin .AND. absar.GE.
721  $ max( safmin, eps*absai, eps*absb ) ) THEN
722  ilimit = .true.
723  scale = max( scale, ( onepls*safmin / anrm1 ) /
724  $ max( onepls*safmin, anrm2*absar ) )
725  END IF
726 *
727 * Check for significant underflow in BETA
728 *
729  IF( abs( sbeta ).LT.safmin .AND. absb.GE.
730  $ max( safmin, eps*absar, eps*absai ) ) THEN
731  ilimit = .true.
732  scale = max( scale, ( onepls*safmin / bnrm1 ) /
733  $ max( onepls*safmin, bnrm2*absb ) )
734  END IF
735 *
736 * Check for possible overflow when limiting scaling
737 *
738  IF( ilimit ) THEN
739  temp = ( scale*safmin )*max( abs( salfar ), abs( salfai ),
740  $ abs( sbeta ) )
741  IF( temp.GT.one )
742  $ scale = scale / temp
743  IF( scale.LT.one )
744  $ ilimit = .false.
745  END IF
746 *
747 * Recompute un-scaled ALPHAR, ALPHAI, BETA if necessary.
748 *
749  IF( ilimit ) THEN
750  salfar = ( scale*alphar( jc ) )*anrm
751  salfai = ( scale*alphai( jc ) )*anrm
752  sbeta = ( scale*beta( jc ) )*bnrm
753  END IF
754  alphar( jc ) = salfar
755  alphai( jc ) = salfai
756  beta( jc ) = sbeta
757  110 CONTINUE
758 *
759  120 CONTINUE
760  work( 1 ) = lwkopt
761 *
762  RETURN
763 *
764 * End of DGEGV
765 *
766  END
subroutine dlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
DLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: dlascl.f:143
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:103
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values.
Definition: dlaset.f:110
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
subroutine dggbak(JOB, SIDE, N, ILO, IHI, LSCALE, RSCALE, M, V, LDV, INFO)
DGGBAK
Definition: dggbak.f:147
subroutine dggbal(JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE, WORK, INFO)
DGGBAL
Definition: dggbal.f:177
subroutine dhgeqz(JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT, ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, INFO)
DHGEQZ
Definition: dhgeqz.f:304
subroutine dgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
DGEQRF
Definition: dgeqrf.f:145
subroutine dtgevc(SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL, LDVL, VR, LDVR, MM, M, WORK, INFO)
DTGEVC
Definition: dtgevc.f:295
subroutine dgegv(JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR, WORK, LWORK, INFO)
DGEEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices
Definition: dgegv.f:306
subroutine dorgqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
DORGQR
Definition: dorgqr.f:128
subroutine dgghrd(COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, LDQ, Z, LDZ, INFO)
DGGHRD
Definition: dgghrd.f:207
subroutine dormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
DORMQR
Definition: dormqr.f:167