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