LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
dggevx.f
Go to the documentation of this file.
1 *> \brief <b> DGGEVX 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 DGGEVX + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dggevx.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dggevx.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dggevx.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DGGEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, B, LDB,
22 * ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR, ILO,
23 * IHI, LSCALE, RSCALE, ABNRM, BBNRM, RCONDE,
24 * RCONDV, WORK, LWORK, IWORK, BWORK, INFO )
25 *
26 * .. Scalar Arguments ..
27 * CHARACTER BALANC, JOBVL, JOBVR, SENSE
28 * INTEGER IHI, ILO, INFO, LDA, LDB, LDVL, LDVR, LWORK, N
29 * DOUBLE PRECISION ABNRM, BBNRM
30 * ..
31 * .. Array Arguments ..
32 * LOGICAL BWORK( * )
33 * INTEGER IWORK( * )
34 * DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
35 * $ B( LDB, * ), BETA( * ), LSCALE( * ),
36 * $ RCONDE( * ), RCONDV( * ), RSCALE( * ),
37 * $ VL( LDVL, * ), VR( LDVR, * ), WORK( * )
38 * ..
39 *
40 *
41 *> \par Purpose:
42 * =============
43 *>
44 *> \verbatim
45 *>
46 *> DGGEVX computes for a pair of N-by-N real nonsymmetric matrices (A,B)
47 *> the generalized eigenvalues, and optionally, the left and/or right
48 *> generalized eigenvectors.
49 *>
50 *> Optionally also, it computes a balancing transformation to improve
51 *> the conditioning of the eigenvalues and eigenvectors (ILO, IHI,
52 *> LSCALE, RSCALE, ABNRM, and BBNRM), reciprocal condition numbers for
53 *> the eigenvalues (RCONDE), and reciprocal condition numbers for the
54 *> right eigenvectors (RCONDV).
55 *>
56 *> A generalized eigenvalue for a pair of matrices (A,B) is a scalar
57 *> lambda or a ratio alpha/beta = lambda, such that A - lambda*B is
58 *> singular. It is usually represented as the pair (alpha,beta), as
59 *> there is a reasonable interpretation for beta=0, and even for both
60 *> being zero.
61 *>
62 *> The right eigenvector v(j) corresponding to the eigenvalue lambda(j)
63 *> of (A,B) satisfies
64 *>
65 *> A * v(j) = lambda(j) * B * v(j) .
66 *>
67 *> The left eigenvector u(j) corresponding to the eigenvalue lambda(j)
68 *> of (A,B) satisfies
69 *>
70 *> u(j)**H * A = lambda(j) * u(j)**H * B.
71 *>
72 *> where u(j)**H is the conjugate-transpose of u(j).
73 *>
74 *> \endverbatim
75 *
76 * Arguments:
77 * ==========
78 *
79 *> \param[in] BALANC
80 *> \verbatim
81 *> BALANC is CHARACTER*1
82 *> Specifies the balance option to be performed.
83 *> = 'N': do not diagonally scale or permute;
84 *> = 'P': permute only;
85 *> = 'S': scale only;
86 *> = 'B': both permute and scale.
87 *> Computed reciprocal condition numbers will be for the
88 *> matrices after permuting and/or balancing. Permuting does
89 *> not change condition numbers (in exact arithmetic), but
90 *> balancing does.
91 *> \endverbatim
92 *>
93 *> \param[in] JOBVL
94 *> \verbatim
95 *> JOBVL is CHARACTER*1
96 *> = 'N': do not compute the left generalized eigenvectors;
97 *> = 'V': compute the left generalized eigenvectors.
98 *> \endverbatim
99 *>
100 *> \param[in] JOBVR
101 *> \verbatim
102 *> JOBVR is CHARACTER*1
103 *> = 'N': do not compute the right generalized eigenvectors;
104 *> = 'V': compute the right generalized eigenvectors.
105 *> \endverbatim
106 *>
107 *> \param[in] SENSE
108 *> \verbatim
109 *> SENSE is CHARACTER*1
110 *> Determines which reciprocal condition numbers are computed.
111 *> = 'N': none are computed;
112 *> = 'E': computed for eigenvalues only;
113 *> = 'V': computed for eigenvectors only;
114 *> = 'B': computed for eigenvalues and eigenvectors.
115 *> \endverbatim
116 *>
117 *> \param[in] N
118 *> \verbatim
119 *> N is INTEGER
120 *> The order of the matrices A, B, VL, and VR. N >= 0.
121 *> \endverbatim
122 *>
123 *> \param[in,out] A
124 *> \verbatim
125 *> A is DOUBLE PRECISION array, dimension (LDA, N)
126 *> On entry, the matrix A in the pair (A,B).
127 *> On exit, A has been overwritten. If JOBVL='V' or JOBVR='V'
128 *> or both, then A contains the first part of the real Schur
129 *> form of the "balanced" versions of the input A and B.
130 *> \endverbatim
131 *>
132 *> \param[in] LDA
133 *> \verbatim
134 *> LDA is INTEGER
135 *> The leading dimension of A. LDA >= max(1,N).
136 *> \endverbatim
137 *>
138 *> \param[in,out] B
139 *> \verbatim
140 *> B is DOUBLE PRECISION array, dimension (LDB, N)
141 *> On entry, the matrix B in the pair (A,B).
142 *> On exit, B has been overwritten. If JOBVL='V' or JOBVR='V'
143 *> or both, then B contains the second part of the real Schur
144 *> form of the "balanced" versions of the input A and B.
145 *> \endverbatim
146 *>
147 *> \param[in] LDB
148 *> \verbatim
149 *> LDB is INTEGER
150 *> The leading dimension of B. LDB >= max(1,N).
151 *> \endverbatim
152 *>
153 *> \param[out] ALPHAR
154 *> \verbatim
155 *> ALPHAR is DOUBLE PRECISION array, dimension (N)
156 *> \endverbatim
157 *>
158 *> \param[out] ALPHAI
159 *> \verbatim
160 *> ALPHAI is DOUBLE PRECISION array, dimension (N)
161 *> \endverbatim
162 *>
163 *> \param[out] BETA
164 *> \verbatim
165 *> BETA is DOUBLE PRECISION array, dimension (N)
166 *> On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will
167 *> be the generalized eigenvalues. If ALPHAI(j) is zero, then
168 *> the j-th eigenvalue is real; if positive, then the j-th and
169 *> (j+1)-st eigenvalues are a complex conjugate pair, with
170 *> ALPHAI(j+1) negative.
171 *>
172 *> Note: the quotients ALPHAR(j)/BETA(j) and ALPHAI(j)/BETA(j)
173 *> may easily over- or underflow, and BETA(j) may even be zero.
174 *> Thus, the user should avoid naively computing the ratio
175 *> ALPHA/BETA. However, ALPHAR and ALPHAI will be always less
176 *> than and usually comparable with norm(A) in magnitude, and
177 *> BETA always less than and usually comparable with norm(B).
178 *> \endverbatim
179 *>
180 *> \param[out] VL
181 *> \verbatim
182 *> VL is DOUBLE PRECISION array, dimension (LDVL,N)
183 *> If JOBVL = 'V', the left eigenvectors u(j) are stored one
184 *> after another in the columns of VL, in the same order as
185 *> their eigenvalues. If the j-th eigenvalue is real, then
186 *> u(j) = VL(:,j), the j-th column of VL. If the j-th and
187 *> (j+1)-th eigenvalues form a complex conjugate pair, then
188 *> u(j) = VL(:,j)+i*VL(:,j+1) and u(j+1) = VL(:,j)-i*VL(:,j+1).
189 *> Each eigenvector will be scaled so the largest component have
190 *> abs(real part) + abs(imag. part) = 1.
191 *> Not referenced if JOBVL = 'N'.
192 *> \endverbatim
193 *>
194 *> \param[in] LDVL
195 *> \verbatim
196 *> LDVL is INTEGER
197 *> The leading dimension of the matrix VL. LDVL >= 1, and
198 *> if JOBVL = 'V', LDVL >= N.
199 *> \endverbatim
200 *>
201 *> \param[out] VR
202 *> \verbatim
203 *> VR is DOUBLE PRECISION array, dimension (LDVR,N)
204 *> If JOBVR = 'V', the right eigenvectors v(j) are stored one
205 *> after another in the columns of VR, in the same order as
206 *> their eigenvalues. If the j-th eigenvalue is real, then
207 *> v(j) = VR(:,j), the j-th column of VR. If the j-th and
208 *> (j+1)-th eigenvalues form a complex conjugate pair, then
209 *> v(j) = VR(:,j)+i*VR(:,j+1) and v(j+1) = VR(:,j)-i*VR(:,j+1).
210 *> Each eigenvector will be scaled so the largest component have
211 *> abs(real part) + abs(imag. part) = 1.
212 *> Not referenced if JOBVR = 'N'.
213 *> \endverbatim
214 *>
215 *> \param[in] LDVR
216 *> \verbatim
217 *> LDVR is INTEGER
218 *> The leading dimension of the matrix VR. LDVR >= 1, and
219 *> if JOBVR = 'V', LDVR >= N.
220 *> \endverbatim
221 *>
222 *> \param[out] ILO
223 *> \verbatim
224 *> ILO is INTEGER
225 *> \endverbatim
226 *>
227 *> \param[out] IHI
228 *> \verbatim
229 *> IHI is INTEGER
230 *> ILO and IHI are integer values such that on exit
231 *> A(i,j) = 0 and B(i,j) = 0 if i > j and
232 *> j = 1,...,ILO-1 or i = IHI+1,...,N.
233 *> If BALANC = 'N' or 'S', ILO = 1 and IHI = N.
234 *> \endverbatim
235 *>
236 *> \param[out] LSCALE
237 *> \verbatim
238 *> LSCALE is DOUBLE PRECISION array, dimension (N)
239 *> Details of the permutations and scaling factors applied
240 *> to the left side of A and B. If PL(j) is the index of the
241 *> row interchanged with row j, and DL(j) is the scaling
242 *> factor applied to row j, then
243 *> LSCALE(j) = PL(j) for j = 1,...,ILO-1
244 *> = DL(j) for j = ILO,...,IHI
245 *> = PL(j) for j = IHI+1,...,N.
246 *> The order in which the interchanges are made is N to IHI+1,
247 *> then 1 to ILO-1.
248 *> \endverbatim
249 *>
250 *> \param[out] RSCALE
251 *> \verbatim
252 *> RSCALE is DOUBLE PRECISION array, dimension (N)
253 *> Details of the permutations and scaling factors applied
254 *> to the right side of A and B. If PR(j) is the index of the
255 *> column interchanged with column j, and DR(j) is the scaling
256 *> factor applied to column j, then
257 *> RSCALE(j) = PR(j) for j = 1,...,ILO-1
258 *> = DR(j) for j = ILO,...,IHI
259 *> = PR(j) for j = IHI+1,...,N
260 *> The order in which the interchanges are made is N to IHI+1,
261 *> then 1 to ILO-1.
262 *> \endverbatim
263 *>
264 *> \param[out] ABNRM
265 *> \verbatim
266 *> ABNRM is DOUBLE PRECISION
267 *> The one-norm of the balanced matrix A.
268 *> \endverbatim
269 *>
270 *> \param[out] BBNRM
271 *> \verbatim
272 *> BBNRM is DOUBLE PRECISION
273 *> The one-norm of the balanced matrix B.
274 *> \endverbatim
275 *>
276 *> \param[out] RCONDE
277 *> \verbatim
278 *> RCONDE is DOUBLE PRECISION array, dimension (N)
279 *> If SENSE = 'E' or 'B', the reciprocal condition numbers of
280 *> the eigenvalues, stored in consecutive elements of the array.
281 *> For a complex conjugate pair of eigenvalues two consecutive
282 *> elements of RCONDE are set to the same value. Thus RCONDE(j),
283 *> RCONDV(j), and the j-th columns of VL and VR all correspond
284 *> to the j-th eigenpair.
285 *> If SENSE = 'N or 'V', RCONDE is not referenced.
286 *> \endverbatim
287 *>
288 *> \param[out] RCONDV
289 *> \verbatim
290 *> RCONDV is DOUBLE PRECISION array, dimension (N)
291 *> If SENSE = 'V' or 'B', the estimated reciprocal condition
292 *> numbers of the eigenvectors, stored in consecutive elements
293 *> of the array. For a complex eigenvector two consecutive
294 *> elements of RCONDV are set to the same value. If the
295 *> eigenvalues cannot be reordered to compute RCONDV(j),
296 *> RCONDV(j) is set to 0; this can only occur when the true
297 *> value would be very small anyway.
298 *> If SENSE = 'N' or 'E', RCONDV is not referenced.
299 *> \endverbatim
300 *>
301 *> \param[out] WORK
302 *> \verbatim
303 *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
304 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
305 *> \endverbatim
306 *>
307 *> \param[in] LWORK
308 *> \verbatim
309 *> LWORK is INTEGER
310 *> The dimension of the array WORK. LWORK >= max(1,2*N).
311 *> If BALANC = 'S' or 'B', or JOBVL = 'V', or JOBVR = 'V',
312 *> LWORK >= max(1,6*N).
313 *> If SENSE = 'E' or 'B', LWORK >= max(1,10*N).
314 *> If SENSE = 'V' or 'B', LWORK >= 2*N*N+8*N+16.
315 *>
316 *> If LWORK = -1, then a workspace query is assumed; the routine
317 *> only calculates the optimal size of the WORK array, returns
318 *> this value as the first entry of the WORK array, and no error
319 *> message related to LWORK is issued by XERBLA.
320 *> \endverbatim
321 *>
322 *> \param[out] IWORK
323 *> \verbatim
324 *> IWORK is INTEGER array, dimension (N+6)
325 *> If SENSE = 'E', IWORK is not referenced.
326 *> \endverbatim
327 *>
328 *> \param[out] BWORK
329 *> \verbatim
330 *> BWORK is LOGICAL array, dimension (N)
331 *> If SENSE = 'N', BWORK is not referenced.
332 *> \endverbatim
333 *>
334 *> \param[out] INFO
335 *> \verbatim
336 *> INFO is INTEGER
337 *> = 0: successful exit
338 *> < 0: if INFO = -i, the i-th argument had an illegal value.
339 *> = 1,...,N:
340 *> The QZ iteration failed. No eigenvectors have been
341 *> calculated, but ALPHAR(j), ALPHAI(j), and BETA(j)
342 *> should be correct for j=INFO+1,...,N.
343 *> > N: =N+1: other than QZ iteration failed in DHGEQZ.
344 *> =N+2: error return from DTGEVC.
345 *> \endverbatim
346 *
347 * Authors:
348 * ========
349 *
350 *> \author Univ. of Tennessee
351 *> \author Univ. of California Berkeley
352 *> \author Univ. of Colorado Denver
353 *> \author NAG Ltd.
354 *
355 *> \date April 2012
356 *
357 *> \ingroup doubleGEeigen
358 *
359 *> \par Further Details:
360 * =====================
361 *>
362 *> \verbatim
363 *>
364 *> Balancing a matrix pair (A,B) includes, first, permuting rows and
365 *> columns to isolate eigenvalues, second, applying diagonal similarity
366 *> transformation to the rows and columns to make the rows and columns
367 *> as close in norm as possible. The computed reciprocal condition
368 *> numbers correspond to the balanced matrix. Permuting rows and columns
369 *> will not change the condition numbers (in exact arithmetic) but
370 *> diagonal scaling will. For further explanation of balancing, see
371 *> section 4.11.1.2 of LAPACK Users' Guide.
372 *>
373 *> An approximate error bound on the chordal distance between the i-th
374 *> computed generalized eigenvalue w and the corresponding exact
375 *> eigenvalue lambda is
376 *>
377 *> chord(w, lambda) <= EPS * norm(ABNRM, BBNRM) / RCONDE(I)
378 *>
379 *> An approximate error bound for the angle between the i-th computed
380 *> eigenvector VL(i) or VR(i) is given by
381 *>
382 *> EPS * norm(ABNRM, BBNRM) / DIF(i).
383 *>
384 *> For further explanation of the reciprocal condition numbers RCONDE
385 *> and RCONDV, see section 4.11 of LAPACK User's Guide.
386 *> \endverbatim
387 *>
388 * =====================================================================
389  SUBROUTINE dggevx( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, B, LDB,
390  $ alphar, alphai, beta, vl, ldvl, vr, ldvr, ilo,
391  $ ihi, lscale, rscale, abnrm, bbnrm, rconde,
392  $ rcondv, work, lwork, iwork, bwork, info )
393 *
394 * -- LAPACK driver routine (version 3.4.1) --
395 * -- LAPACK is a software package provided by Univ. of Tennessee, --
396 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
397 * April 2012
398 *
399 * .. Scalar Arguments ..
400  CHARACTER balanc, jobvl, jobvr, sense
401  INTEGER ihi, ilo, info, lda, ldb, ldvl, ldvr, lwork, n
402  DOUBLE PRECISION abnrm, bbnrm
403 * ..
404 * .. Array Arguments ..
405  LOGICAL bwork( * )
406  INTEGER iwork( * )
407  DOUBLE PRECISION a( lda, * ), alphai( * ), alphar( * ),
408  $ b( ldb, * ), beta( * ), lscale( * ),
409  $ rconde( * ), rcondv( * ), rscale( * ),
410  $ vl( ldvl, * ), vr( ldvr, * ), work( * )
411 * ..
412 *
413 * =====================================================================
414 *
415 * .. Parameters ..
416  DOUBLE PRECISION zero, one
417  parameter( zero = 0.0d+0, one = 1.0d+0 )
418 * ..
419 * .. Local Scalars ..
420  LOGICAL ilascl, ilbscl, ilv, ilvl, ilvr, lquery, noscl,
421  $ pair, wantsb, wantse, wantsn, wantsv
422  CHARACTER chtemp
423  INTEGER i, icols, ierr, ijobvl, ijobvr, in, irows,
424  $ itau, iwrk, iwrk1, j, jc, jr, m, maxwrk,
425  $ minwrk, mm
426  DOUBLE PRECISION anrm, anrmto, bignum, bnrm, bnrmto, eps,
427  $ smlnum, temp
428 * ..
429 * .. Local Arrays ..
430  LOGICAL ldumma( 1 )
431 * ..
432 * .. External Subroutines ..
433  EXTERNAL dgeqrf, dggbak, dggbal, dgghrd, dhgeqz, dlabad,
435  $ dtgsna, xerbla
436 * ..
437 * .. External Functions ..
438  LOGICAL lsame
439  INTEGER ilaenv
440  DOUBLE PRECISION dlamch, dlange
441  EXTERNAL lsame, ilaenv, dlamch, dlange
442 * ..
443 * .. Intrinsic Functions ..
444  INTRINSIC abs, max, sqrt
445 * ..
446 * .. Executable Statements ..
447 *
448 * Decode the input arguments
449 *
450  IF( lsame( jobvl, 'N' ) ) THEN
451  ijobvl = 1
452  ilvl = .false.
453  ELSE IF( lsame( jobvl, 'V' ) ) THEN
454  ijobvl = 2
455  ilvl = .true.
456  ELSE
457  ijobvl = -1
458  ilvl = .false.
459  END IF
460 *
461  IF( lsame( jobvr, 'N' ) ) THEN
462  ijobvr = 1
463  ilvr = .false.
464  ELSE IF( lsame( jobvr, 'V' ) ) THEN
465  ijobvr = 2
466  ilvr = .true.
467  ELSE
468  ijobvr = -1
469  ilvr = .false.
470  END IF
471  ilv = ilvl .OR. ilvr
472 *
473  noscl = lsame( balanc, 'N' ) .OR. lsame( balanc, 'P' )
474  wantsn = lsame( sense, 'N' )
475  wantse = lsame( sense, 'E' )
476  wantsv = lsame( sense, 'V' )
477  wantsb = lsame( sense, 'B' )
478 *
479 * Test the input arguments
480 *
481  info = 0
482  lquery = ( lwork.EQ.-1 )
483  IF( .NOT.( lsame( balanc, 'N' ) .OR. lsame( balanc,
484  $ 'S' ) .OR. lsame( balanc, 'P' ) .OR. lsame( balanc, 'B' ) ) )
485  $ THEN
486  info = -1
487  ELSE IF( ijobvl.LE.0 ) THEN
488  info = -2
489  ELSE IF( ijobvr.LE.0 ) THEN
490  info = -3
491  ELSE IF( .NOT.( wantsn .OR. wantse .OR. wantsb .OR. wantsv ) )
492  $ THEN
493  info = -4
494  ELSE IF( n.LT.0 ) THEN
495  info = -5
496  ELSE IF( lda.LT.max( 1, n ) ) THEN
497  info = -7
498  ELSE IF( ldb.LT.max( 1, n ) ) THEN
499  info = -9
500  ELSE IF( ldvl.LT.1 .OR. ( ilvl .AND. ldvl.LT.n ) ) THEN
501  info = -14
502  ELSE IF( ldvr.LT.1 .OR. ( ilvr .AND. ldvr.LT.n ) ) THEN
503  info = -16
504  END IF
505 *
506 * Compute workspace
507 * (Note: Comments in the code beginning "Workspace:" describe the
508 * minimal amount of workspace needed at that point in the code,
509 * as well as the preferred amount for good performance.
510 * NB refers to the optimal block size for the immediately
511 * following subroutine, as returned by ILAENV. The workspace is
512 * computed assuming ILO = 1 and IHI = N, the worst case.)
513 *
514  IF( info.EQ.0 ) THEN
515  IF( n.EQ.0 ) THEN
516  minwrk = 1
517  maxwrk = 1
518  ELSE
519  IF( noscl .AND. .NOT.ilv ) THEN
520  minwrk = 2*n
521  ELSE
522  minwrk = 6*n
523  END IF
524  IF( wantse .OR. wantsb ) THEN
525  minwrk = 10*n
526  END IF
527  IF( wantsv .OR. wantsb ) THEN
528  minwrk = max( minwrk, 2*n*( n + 4 ) + 16 )
529  END IF
530  maxwrk = minwrk
531  maxwrk = max( maxwrk,
532  $ n + n*ilaenv( 1, 'DGEQRF', ' ', n, 1, n, 0 ) )
533  maxwrk = max( maxwrk,
534  $ n + n*ilaenv( 1, 'DORMQR', ' ', n, 1, n, 0 ) )
535  IF( ilvl ) THEN
536  maxwrk = max( maxwrk, n +
537  $ n*ilaenv( 1, 'DORGQR', ' ', n, 1, n, 0 ) )
538  END IF
539  END IF
540  work( 1 ) = maxwrk
541 *
542  IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
543  info = -26
544  END IF
545  END IF
546 *
547  IF( info.NE.0 ) THEN
548  CALL xerbla( 'DGGEVX', -info )
549  return
550  ELSE IF( lquery ) THEN
551  return
552  END IF
553 *
554 * Quick return if possible
555 *
556  IF( n.EQ.0 )
557  $ return
558 *
559 *
560 * Get machine constants
561 *
562  eps = dlamch( 'P' )
563  smlnum = dlamch( 'S' )
564  bignum = one / smlnum
565  CALL dlabad( smlnum, bignum )
566  smlnum = sqrt( smlnum ) / eps
567  bignum = one / smlnum
568 *
569 * Scale A if max element outside range [SMLNUM,BIGNUM]
570 *
571  anrm = dlange( 'M', n, n, a, lda, work )
572  ilascl = .false.
573  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
574  anrmto = smlnum
575  ilascl = .true.
576  ELSE IF( anrm.GT.bignum ) THEN
577  anrmto = bignum
578  ilascl = .true.
579  END IF
580  IF( ilascl )
581  $ CALL dlascl( 'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
582 *
583 * Scale B if max element outside range [SMLNUM,BIGNUM]
584 *
585  bnrm = dlange( 'M', n, n, b, ldb, work )
586  ilbscl = .false.
587  IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
588  bnrmto = smlnum
589  ilbscl = .true.
590  ELSE IF( bnrm.GT.bignum ) THEN
591  bnrmto = bignum
592  ilbscl = .true.
593  END IF
594  IF( ilbscl )
595  $ CALL dlascl( 'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
596 *
597 * Permute and/or balance the matrix pair (A,B)
598 * (Workspace: need 6*N if BALANC = 'S' or 'B', 1 otherwise)
599 *
600  CALL dggbal( balanc, n, a, lda, b, ldb, ilo, ihi, lscale, rscale,
601  $ work, ierr )
602 *
603 * Compute ABNRM and BBNRM
604 *
605  abnrm = dlange( '1', n, n, a, lda, work( 1 ) )
606  IF( ilascl ) THEN
607  work( 1 ) = abnrm
608  CALL dlascl( 'G', 0, 0, anrmto, anrm, 1, 1, work( 1 ), 1,
609  $ ierr )
610  abnrm = work( 1 )
611  END IF
612 *
613  bbnrm = dlange( '1', n, n, b, ldb, work( 1 ) )
614  IF( ilbscl ) THEN
615  work( 1 ) = bbnrm
616  CALL dlascl( 'G', 0, 0, bnrmto, bnrm, 1, 1, work( 1 ), 1,
617  $ ierr )
618  bbnrm = work( 1 )
619  END IF
620 *
621 * Reduce B to triangular form (QR decomposition of B)
622 * (Workspace: need N, prefer N*NB )
623 *
624  irows = ihi + 1 - ilo
625  IF( ilv .OR. .NOT.wantsn ) THEN
626  icols = n + 1 - ilo
627  ELSE
628  icols = irows
629  END IF
630  itau = 1
631  iwrk = itau + irows
632  CALL dgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
633  $ work( iwrk ), lwork+1-iwrk, ierr )
634 *
635 * Apply the orthogonal transformation to A
636 * (Workspace: need N, prefer N*NB)
637 *
638  CALL dormqr( 'L', 'T', irows, icols, irows, b( ilo, ilo ), ldb,
639  $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
640  $ lwork+1-iwrk, ierr )
641 *
642 * Initialize VL and/or VR
643 * (Workspace: need N, prefer N*NB)
644 *
645  IF( ilvl ) THEN
646  CALL dlaset( 'Full', n, n, zero, one, vl, ldvl )
647  IF( irows.GT.1 ) THEN
648  CALL dlacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
649  $ vl( ilo+1, ilo ), ldvl )
650  END IF
651  CALL dorgqr( irows, irows, irows, vl( ilo, ilo ), ldvl,
652  $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
653  END IF
654 *
655  IF( ilvr )
656  $ CALL dlaset( 'Full', n, n, zero, one, vr, ldvr )
657 *
658 * Reduce to generalized Hessenberg form
659 * (Workspace: none needed)
660 *
661  IF( ilv .OR. .NOT.wantsn ) THEN
662 *
663 * Eigenvectors requested -- work on whole matrix.
664 *
665  CALL dgghrd( jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb, vl,
666  $ ldvl, vr, ldvr, ierr )
667  ELSE
668  CALL dgghrd( 'N', 'N', irows, 1, irows, a( ilo, ilo ), lda,
669  $ b( ilo, ilo ), ldb, vl, ldvl, vr, ldvr, ierr )
670  END IF
671 *
672 * Perform QZ algorithm (Compute eigenvalues, and optionally, the
673 * Schur forms and Schur vectors)
674 * (Workspace: need N)
675 *
676  IF( ilv .OR. .NOT.wantsn ) THEN
677  chtemp = 'S'
678  ELSE
679  chtemp = 'E'
680  END IF
681 *
682  CALL dhgeqz( chtemp, jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb,
683  $ alphar, alphai, beta, vl, ldvl, vr, ldvr, work,
684  $ lwork, ierr )
685  IF( ierr.NE.0 ) THEN
686  IF( ierr.GT.0 .AND. ierr.LE.n ) THEN
687  info = ierr
688  ELSE IF( ierr.GT.n .AND. ierr.LE.2*n ) THEN
689  info = ierr - n
690  ELSE
691  info = n + 1
692  END IF
693  go to 130
694  END IF
695 *
696 * Compute Eigenvectors and estimate condition numbers if desired
697 * (Workspace: DTGEVC: need 6*N
698 * DTGSNA: need 2*N*(N+2)+16 if SENSE = 'V' or 'B',
699 * need N otherwise )
700 *
701  IF( ilv .OR. .NOT.wantsn ) THEN
702  IF( ilv ) THEN
703  IF( ilvl ) THEN
704  IF( ilvr ) THEN
705  chtemp = 'B'
706  ELSE
707  chtemp = 'L'
708  END IF
709  ELSE
710  chtemp = 'R'
711  END IF
712 *
713  CALL dtgevc( chtemp, 'B', ldumma, n, a, lda, b, ldb, vl,
714  $ ldvl, vr, ldvr, n, in, work, ierr )
715  IF( ierr.NE.0 ) THEN
716  info = n + 2
717  go to 130
718  END IF
719  END IF
720 *
721  IF( .NOT.wantsn ) THEN
722 *
723 * compute eigenvectors (DTGEVC) and estimate condition
724 * numbers (DTGSNA). Note that the definition of the condition
725 * number is not invariant under transformation (u,v) to
726 * (Q*u, Z*v), where (u,v) are eigenvectors of the generalized
727 * Schur form (S,T), Q and Z are orthogonal matrices. In order
728 * to avoid using extra 2*N*N workspace, we have to recalculate
729 * eigenvectors and estimate one condition numbers at a time.
730 *
731  pair = .false.
732  DO 20 i = 1, n
733 *
734  IF( pair ) THEN
735  pair = .false.
736  go to 20
737  END IF
738  mm = 1
739  IF( i.LT.n ) THEN
740  IF( a( i+1, i ).NE.zero ) THEN
741  pair = .true.
742  mm = 2
743  END IF
744  END IF
745 *
746  DO 10 j = 1, n
747  bwork( j ) = .false.
748  10 continue
749  IF( mm.EQ.1 ) THEN
750  bwork( i ) = .true.
751  ELSE IF( mm.EQ.2 ) THEN
752  bwork( i ) = .true.
753  bwork( i+1 ) = .true.
754  END IF
755 *
756  iwrk = mm*n + 1
757  iwrk1 = iwrk + mm*n
758 *
759 * Compute a pair of left and right eigenvectors.
760 * (compute workspace: need up to 4*N + 6*N)
761 *
762  IF( wantse .OR. wantsb ) THEN
763  CALL dtgevc( 'B', 'S', bwork, n, a, lda, b, ldb,
764  $ work( 1 ), n, work( iwrk ), n, mm, m,
765  $ work( iwrk1 ), ierr )
766  IF( ierr.NE.0 ) THEN
767  info = n + 2
768  go to 130
769  END IF
770  END IF
771 *
772  CALL dtgsna( sense, 'S', bwork, n, a, lda, b, ldb,
773  $ work( 1 ), n, work( iwrk ), n, rconde( i ),
774  $ rcondv( i ), mm, m, work( iwrk1 ),
775  $ lwork-iwrk1+1, iwork, ierr )
776 *
777  20 continue
778  END IF
779  END IF
780 *
781 * Undo balancing on VL and VR and normalization
782 * (Workspace: none needed)
783 *
784  IF( ilvl ) THEN
785  CALL dggbak( balanc, 'L', n, ilo, ihi, lscale, rscale, n, vl,
786  $ ldvl, ierr )
787 *
788  DO 70 jc = 1, n
789  IF( alphai( jc ).LT.zero )
790  $ go to 70
791  temp = zero
792  IF( alphai( jc ).EQ.zero ) THEN
793  DO 30 jr = 1, n
794  temp = max( temp, abs( vl( jr, jc ) ) )
795  30 continue
796  ELSE
797  DO 40 jr = 1, n
798  temp = max( temp, abs( vl( jr, jc ) )+
799  $ abs( vl( jr, jc+1 ) ) )
800  40 continue
801  END IF
802  IF( temp.LT.smlnum )
803  $ go to 70
804  temp = one / temp
805  IF( alphai( jc ).EQ.zero ) THEN
806  DO 50 jr = 1, n
807  vl( jr, jc ) = vl( jr, jc )*temp
808  50 continue
809  ELSE
810  DO 60 jr = 1, n
811  vl( jr, jc ) = vl( jr, jc )*temp
812  vl( jr, jc+1 ) = vl( jr, jc+1 )*temp
813  60 continue
814  END IF
815  70 continue
816  END IF
817  IF( ilvr ) THEN
818  CALL dggbak( balanc, 'R', n, ilo, ihi, lscale, rscale, n, vr,
819  $ ldvr, ierr )
820  DO 120 jc = 1, n
821  IF( alphai( jc ).LT.zero )
822  $ go to 120
823  temp = zero
824  IF( alphai( jc ).EQ.zero ) THEN
825  DO 80 jr = 1, n
826  temp = max( temp, abs( vr( jr, jc ) ) )
827  80 continue
828  ELSE
829  DO 90 jr = 1, n
830  temp = max( temp, abs( vr( jr, jc ) )+
831  $ abs( vr( jr, jc+1 ) ) )
832  90 continue
833  END IF
834  IF( temp.LT.smlnum )
835  $ go to 120
836  temp = one / temp
837  IF( alphai( jc ).EQ.zero ) THEN
838  DO 100 jr = 1, n
839  vr( jr, jc ) = vr( jr, jc )*temp
840  100 continue
841  ELSE
842  DO 110 jr = 1, n
843  vr( jr, jc ) = vr( jr, jc )*temp
844  vr( jr, jc+1 ) = vr( jr, jc+1 )*temp
845  110 continue
846  END IF
847  120 continue
848  END IF
849 *
850 * Undo scaling if necessary
851 *
852  130 continue
853 *
854  IF( ilascl ) THEN
855  CALL dlascl( 'G', 0, 0, anrmto, anrm, n, 1, alphar, n, ierr )
856  CALL dlascl( 'G', 0, 0, anrmto, anrm, n, 1, alphai, n, ierr )
857  END IF
858 *
859  IF( ilbscl ) THEN
860  CALL dlascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
861  END IF
862 *
863  work( 1 ) = maxwrk
864  return
865 *
866 * End of DGGEVX
867 *
868  END