LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
zggesx.f
Go to the documentation of this file.
1 *> \brief <b> ZGGESX computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors 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 ZGGESX + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zggesx.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zggesx.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zggesx.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE ZGGESX( JOBVSL, JOBVSR, SORT, SELCTG, SENSE, N, A, LDA,
22 * B, LDB, SDIM, ALPHA, BETA, VSL, LDVSL, VSR,
23 * LDVSR, RCONDE, RCONDV, WORK, LWORK, RWORK,
24 * IWORK, LIWORK, BWORK, INFO )
25 *
26 * .. Scalar Arguments ..
27 * CHARACTER JOBVSL, JOBVSR, SENSE, SORT
28 * INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LIWORK, LWORK, N,
29 * $ SDIM
30 * ..
31 * .. Array Arguments ..
32 * LOGICAL BWORK( * )
33 * INTEGER IWORK( * )
34 * DOUBLE PRECISION RCONDE( 2 ), RCONDV( 2 ), RWORK( * )
35 * COMPLEX*16 A( LDA, * ), ALPHA( * ), B( LDB, * ),
36 * $ BETA( * ), VSL( LDVSL, * ), VSR( LDVSR, * ),
37 * $ WORK( * )
38 * ..
39 * .. Function Arguments ..
40 * LOGICAL SELCTG
41 * EXTERNAL SELCTG
42 * ..
43 *
44 *
45 *> \par Purpose:
46 * =============
47 *>
48 *> \verbatim
49 *>
50 *> ZGGESX computes for a pair of N-by-N complex nonsymmetric matrices
51 *> (A,B), the generalized eigenvalues, the complex Schur form (S,T),
52 *> and, optionally, the left and/or right matrices of Schur vectors (VSL
53 *> and VSR). This gives the generalized Schur factorization
54 *>
55 *> (A,B) = ( (VSL) S (VSR)**H, (VSL) T (VSR)**H )
56 *>
57 *> where (VSR)**H is the conjugate-transpose of VSR.
58 *>
59 *> Optionally, it also orders the eigenvalues so that a selected cluster
60 *> of eigenvalues appears in the leading diagonal blocks of the upper
61 *> triangular matrix S and the upper triangular matrix T; computes
62 *> a reciprocal condition number for the average of the selected
63 *> eigenvalues (RCONDE); and computes a reciprocal condition number for
64 *> the right and left deflating subspaces corresponding to the selected
65 *> eigenvalues (RCONDV). The leading columns of VSL and VSR then form
66 *> an orthonormal basis for the corresponding left and right eigenspaces
67 *> (deflating subspaces).
68 *>
69 *> A generalized eigenvalue for a pair of matrices (A,B) is a scalar w
70 *> or a ratio alpha/beta = w, such that A - w*B is singular. It is
71 *> usually represented as the pair (alpha,beta), as there is a
72 *> reasonable interpretation for beta=0 or for both being zero.
73 *>
74 *> A pair of matrices (S,T) is in generalized complex Schur form if T is
75 *> upper triangular with non-negative diagonal and S is upper
76 *> triangular.
77 *> \endverbatim
78 *
79 * Arguments:
80 * ==========
81 *
82 *> \param[in] JOBVSL
83 *> \verbatim
84 *> JOBVSL is CHARACTER*1
85 *> = 'N': do not compute the left Schur vectors;
86 *> = 'V': compute the left Schur vectors.
87 *> \endverbatim
88 *>
89 *> \param[in] JOBVSR
90 *> \verbatim
91 *> JOBVSR is CHARACTER*1
92 *> = 'N': do not compute the right Schur vectors;
93 *> = 'V': compute the right Schur vectors.
94 *> \endverbatim
95 *>
96 *> \param[in] SORT
97 *> \verbatim
98 *> SORT is CHARACTER*1
99 *> Specifies whether or not to order the eigenvalues on the
100 *> diagonal of the generalized Schur form.
101 *> = 'N': Eigenvalues are not ordered;
102 *> = 'S': Eigenvalues are ordered (see SELCTG).
103 *> \endverbatim
104 *>
105 *> \param[in] SELCTG
106 *> \verbatim
107 *> SELCTG is procedure) LOGICAL FUNCTION of two COMPLEX*16 arguments
108 *> SELCTG must be declared EXTERNAL in the calling subroutine.
109 *> If SORT = 'N', SELCTG is not referenced.
110 *> If SORT = 'S', SELCTG is used to select eigenvalues to sort
111 *> to the top left of the Schur form.
112 *> Note that a selected complex eigenvalue may no longer satisfy
113 *> SELCTG(ALPHA(j),BETA(j)) = .TRUE. after ordering, since
114 *> ordering may change the value of complex eigenvalues
115 *> (especially if the eigenvalue is ill-conditioned), in this
116 *> case INFO is set to N+3 see INFO below).
117 *> \endverbatim
118 *>
119 *> \param[in] SENSE
120 *> \verbatim
121 *> SENSE is CHARACTER*1
122 *> Determines which reciprocal condition numbers are computed.
123 *> = 'N' : None are computed;
124 *> = 'E' : Computed for average of selected eigenvalues only;
125 *> = 'V' : Computed for selected deflating subspaces only;
126 *> = 'B' : Computed for both.
127 *> If SENSE = 'E', 'V', or 'B', SORT must equal 'S'.
128 *> \endverbatim
129 *>
130 *> \param[in] N
131 *> \verbatim
132 *> N is INTEGER
133 *> The order of the matrices A, B, VSL, and VSR. N >= 0.
134 *> \endverbatim
135 *>
136 *> \param[in,out] A
137 *> \verbatim
138 *> A is COMPLEX*16 array, dimension (LDA, N)
139 *> On entry, the first of the pair of matrices.
140 *> On exit, A has been overwritten by its generalized Schur
141 *> form S.
142 *> \endverbatim
143 *>
144 *> \param[in] LDA
145 *> \verbatim
146 *> LDA is INTEGER
147 *> The leading dimension of A. LDA >= max(1,N).
148 *> \endverbatim
149 *>
150 *> \param[in,out] B
151 *> \verbatim
152 *> B is COMPLEX*16 array, dimension (LDB, N)
153 *> On entry, the second of the pair of matrices.
154 *> On exit, B has been overwritten by its generalized Schur
155 *> form T.
156 *> \endverbatim
157 *>
158 *> \param[in] LDB
159 *> \verbatim
160 *> LDB is INTEGER
161 *> The leading dimension of B. LDB >= max(1,N).
162 *> \endverbatim
163 *>
164 *> \param[out] SDIM
165 *> \verbatim
166 *> SDIM is INTEGER
167 *> If SORT = 'N', SDIM = 0.
168 *> If SORT = 'S', SDIM = number of eigenvalues (after sorting)
169 *> for which SELCTG is true.
170 *> \endverbatim
171 *>
172 *> \param[out] ALPHA
173 *> \verbatim
174 *> ALPHA is COMPLEX*16 array, dimension (N)
175 *> \endverbatim
176 *>
177 *> \param[out] BETA
178 *> \verbatim
179 *> BETA is COMPLEX*16 array, dimension (N)
180 *> On exit, ALPHA(j)/BETA(j), j=1,...,N, will be the
181 *> generalized eigenvalues. ALPHA(j) and BETA(j),j=1,...,N are
182 *> the diagonals of the complex Schur form (S,T). BETA(j) will
183 *> be non-negative real.
184 *>
185 *> Note: the quotients ALPHA(j)/BETA(j) may easily over- or
186 *> underflow, and BETA(j) may even be zero. Thus, the user
187 *> should avoid naively computing the ratio alpha/beta.
188 *> However, ALPHA will be always less than and usually
189 *> comparable with norm(A) in magnitude, and BETA always less
190 *> than and usually comparable with norm(B).
191 *> \endverbatim
192 *>
193 *> \param[out] VSL
194 *> \verbatim
195 *> VSL is COMPLEX*16 array, dimension (LDVSL,N)
196 *> If JOBVSL = 'V', VSL will contain the left Schur vectors.
197 *> Not referenced if JOBVSL = 'N'.
198 *> \endverbatim
199 *>
200 *> \param[in] LDVSL
201 *> \verbatim
202 *> LDVSL is INTEGER
203 *> The leading dimension of the matrix VSL. LDVSL >=1, and
204 *> if JOBVSL = 'V', LDVSL >= N.
205 *> \endverbatim
206 *>
207 *> \param[out] VSR
208 *> \verbatim
209 *> VSR is COMPLEX*16 array, dimension (LDVSR,N)
210 *> If JOBVSR = 'V', VSR will contain the right Schur vectors.
211 *> Not referenced if JOBVSR = 'N'.
212 *> \endverbatim
213 *>
214 *> \param[in] LDVSR
215 *> \verbatim
216 *> LDVSR is INTEGER
217 *> The leading dimension of the matrix VSR. LDVSR >= 1, and
218 *> if JOBVSR = 'V', LDVSR >= N.
219 *> \endverbatim
220 *>
221 *> \param[out] RCONDE
222 *> \verbatim
223 *> RCONDE is DOUBLE PRECISION array, dimension ( 2 )
224 *> If SENSE = 'E' or 'B', RCONDE(1) and RCONDE(2) contain the
225 *> reciprocal condition numbers for the average of the selected
226 *> eigenvalues.
227 *> Not referenced if SENSE = 'N' or 'V'.
228 *> \endverbatim
229 *>
230 *> \param[out] RCONDV
231 *> \verbatim
232 *> RCONDV is DOUBLE PRECISION array, dimension ( 2 )
233 *> If SENSE = 'V' or 'B', RCONDV(1) and RCONDV(2) contain the
234 *> reciprocal condition number for the selected deflating
235 *> subspaces.
236 *> Not referenced if SENSE = 'N' or 'E'.
237 *> \endverbatim
238 *>
239 *> \param[out] WORK
240 *> \verbatim
241 *> WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
242 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
243 *> \endverbatim
244 *>
245 *> \param[in] LWORK
246 *> \verbatim
247 *> LWORK is INTEGER
248 *> The dimension of the array WORK.
249 *> If N = 0, LWORK >= 1, else if SENSE = 'E', 'V', or 'B',
250 *> LWORK >= MAX(1,2*N,2*SDIM*(N-SDIM)), else
251 *> LWORK >= MAX(1,2*N). Note that 2*SDIM*(N-SDIM) <= N*N/2.
252 *> Note also that an error is only returned if
253 *> LWORK < MAX(1,2*N), but if SENSE = 'E' or 'V' or 'B' this may
254 *> not be large enough.
255 *>
256 *> If LWORK = -1, then a workspace query is assumed; the routine
257 *> only calculates the bound on the optimal size of the WORK
258 *> array and the minimum size of the IWORK array, returns these
259 *> values as the first entries of the WORK and IWORK arrays, and
260 *> no error message related to LWORK or LIWORK is issued by
261 *> XERBLA.
262 *> \endverbatim
263 *>
264 *> \param[out] RWORK
265 *> \verbatim
266 *> RWORK is DOUBLE PRECISION array, dimension ( 8*N )
267 *> Real workspace.
268 *> \endverbatim
269 *>
270 *> \param[out] IWORK
271 *> \verbatim
272 *> IWORK is INTEGER array, dimension (MAX(1,LIWORK))
273 *> On exit, if INFO = 0, IWORK(1) returns the minimum LIWORK.
274 *> \endverbatim
275 *>
276 *> \param[in] LIWORK
277 *> \verbatim
278 *> LIWORK is INTEGER
279 *> The dimension of the array IWORK.
280 *> If SENSE = 'N' or N = 0, LIWORK >= 1, otherwise
281 *> LIWORK >= N+2.
282 *>
283 *> If LIWORK = -1, then a workspace query is assumed; the
284 *> routine only calculates the bound on the optimal size of the
285 *> WORK array and the minimum size of the IWORK array, returns
286 *> these values as the first entries of the WORK and IWORK
287 *> arrays, and no error message related to LWORK or LIWORK is
288 *> issued by XERBLA.
289 *> \endverbatim
290 *>
291 *> \param[out] BWORK
292 *> \verbatim
293 *> BWORK is LOGICAL array, dimension (N)
294 *> Not referenced if SORT = 'N'.
295 *> \endverbatim
296 *>
297 *> \param[out] INFO
298 *> \verbatim
299 *> INFO is INTEGER
300 *> = 0: successful exit
301 *> < 0: if INFO = -i, the i-th argument had an illegal value.
302 *> = 1,...,N:
303 *> The QZ iteration failed. (A,B) are not in Schur
304 *> form, but ALPHA(j) and BETA(j) should be correct for
305 *> j=INFO+1,...,N.
306 *> > N: =N+1: other than QZ iteration failed in ZHGEQZ
307 *> =N+2: after reordering, roundoff changed values of
308 *> some complex eigenvalues so that leading
309 *> eigenvalues in the Generalized Schur form no
310 *> longer satisfy SELCTG=.TRUE. This could also
311 *> be caused due to scaling.
312 *> =N+3: reordering failed in ZTGSEN.
313 *> \endverbatim
314 *
315 * Authors:
316 * ========
317 *
318 *> \author Univ. of Tennessee
319 *> \author Univ. of California Berkeley
320 *> \author Univ. of Colorado Denver
321 *> \author NAG Ltd.
322 *
323 *> \date November 2011
324 *
325 *> \ingroup complex16GEeigen
326 *
327 * =====================================================================
328  SUBROUTINE zggesx( JOBVSL, JOBVSR, SORT, SELCTG, SENSE, N, A, LDA,
329  $ b, ldb, sdim, alpha, beta, vsl, ldvsl, vsr,
330  $ ldvsr, rconde, rcondv, work, lwork, rwork,
331  $ iwork, liwork, bwork, info )
332 *
333 * -- LAPACK driver routine (version 3.4.0) --
334 * -- LAPACK is a software package provided by Univ. of Tennessee, --
335 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
336 * November 2011
337 *
338 * .. Scalar Arguments ..
339  CHARACTER jobvsl, jobvsr, sense, sort
340  INTEGER info, lda, ldb, ldvsl, ldvsr, liwork, lwork, n,
341  $ sdim
342 * ..
343 * .. Array Arguments ..
344  LOGICAL bwork( * )
345  INTEGER iwork( * )
346  DOUBLE PRECISION rconde( 2 ), rcondv( 2 ), rwork( * )
347  COMPLEX*16 a( lda, * ), alpha( * ), b( ldb, * ),
348  $ beta( * ), vsl( ldvsl, * ), vsr( ldvsr, * ),
349  $ work( * )
350 * ..
351 * .. Function Arguments ..
352  LOGICAL selctg
353  EXTERNAL selctg
354 * ..
355 *
356 * =====================================================================
357 *
358 * .. Parameters ..
359  DOUBLE PRECISION zero, one
360  parameter( zero = 0.0d+0, one = 1.0d+0 )
361  COMPLEX*16 czero, cone
362  parameter( czero = ( 0.0d+0, 0.0d+0 ),
363  $ cone = ( 1.0d+0, 0.0d+0 ) )
364 * ..
365 * .. Local Scalars ..
366  LOGICAL cursl, ilascl, ilbscl, ilvsl, ilvsr, lastsl,
367  $ lquery, wantsb, wantse, wantsn, wantst, wantsv
368  INTEGER i, icols, ierr, ihi, ijob, ijobvl, ijobvr,
369  $ ileft, ilo, iright, irows, irwrk, itau, iwrk,
370  $ liwmin, lwrk, maxwrk, minwrk
371  DOUBLE PRECISION anrm, anrmto, bignum, bnrm, bnrmto, eps, pl,
372  $ pr, smlnum
373 * ..
374 * .. Local Arrays ..
375  DOUBLE PRECISION dif( 2 )
376 * ..
377 * .. External Subroutines ..
378  EXTERNAL dlabad, xerbla, zgeqrf, zggbak, zggbal, zgghrd,
380  $ zunmqr
381 * ..
382 * .. External Functions ..
383  LOGICAL lsame
384  INTEGER ilaenv
385  DOUBLE PRECISION dlamch, zlange
386  EXTERNAL lsame, ilaenv, dlamch, zlange
387 * ..
388 * .. Intrinsic Functions ..
389  INTRINSIC max, sqrt
390 * ..
391 * .. Executable Statements ..
392 *
393 * Decode the input arguments
394 *
395  IF( lsame( jobvsl, 'N' ) ) THEN
396  ijobvl = 1
397  ilvsl = .false.
398  ELSE IF( lsame( jobvsl, 'V' ) ) THEN
399  ijobvl = 2
400  ilvsl = .true.
401  ELSE
402  ijobvl = -1
403  ilvsl = .false.
404  END IF
405 *
406  IF( lsame( jobvsr, 'N' ) ) THEN
407  ijobvr = 1
408  ilvsr = .false.
409  ELSE IF( lsame( jobvsr, 'V' ) ) THEN
410  ijobvr = 2
411  ilvsr = .true.
412  ELSE
413  ijobvr = -1
414  ilvsr = .false.
415  END IF
416 *
417  wantst = lsame( sort, 'S' )
418  wantsn = lsame( sense, 'N' )
419  wantse = lsame( sense, 'E' )
420  wantsv = lsame( sense, 'V' )
421  wantsb = lsame( sense, 'B' )
422  lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
423  IF( wantsn ) THEN
424  ijob = 0
425  ELSE IF( wantse ) THEN
426  ijob = 1
427  ELSE IF( wantsv ) THEN
428  ijob = 2
429  ELSE IF( wantsb ) THEN
430  ijob = 4
431  END IF
432 *
433 * Test the input arguments
434 *
435  info = 0
436  IF( ijobvl.LE.0 ) THEN
437  info = -1
438  ELSE IF( ijobvr.LE.0 ) THEN
439  info = -2
440  ELSE IF( ( .NOT.wantst ) .AND. ( .NOT.lsame( sort, 'N' ) ) ) THEN
441  info = -3
442  ELSE IF( .NOT.( wantsn .OR. wantse .OR. wantsv .OR. wantsb ) .OR.
443  $ ( .NOT.wantst .AND. .NOT.wantsn ) ) THEN
444  info = -5
445  ELSE IF( n.LT.0 ) THEN
446  info = -6
447  ELSE IF( lda.LT.max( 1, n ) ) THEN
448  info = -8
449  ELSE IF( ldb.LT.max( 1, n ) ) THEN
450  info = -10
451  ELSE IF( ldvsl.LT.1 .OR. ( ilvsl .AND. ldvsl.LT.n ) ) THEN
452  info = -15
453  ELSE IF( ldvsr.LT.1 .OR. ( ilvsr .AND. ldvsr.LT.n ) ) THEN
454  info = -17
455  END IF
456 *
457 * Compute workspace
458 * (Note: Comments in the code beginning "Workspace:" describe the
459 * minimal amount of workspace needed at that point in the code,
460 * as well as the preferred amount for good performance.
461 * NB refers to the optimal block size for the immediately
462 * following subroutine, as returned by ILAENV.)
463 *
464  IF( info.EQ.0 ) THEN
465  IF( n.GT.0) THEN
466  minwrk = 2*n
467  maxwrk = n*(1 + ilaenv( 1, 'ZGEQRF', ' ', n, 1, n, 0 ) )
468  maxwrk = max( maxwrk, n*( 1 +
469  $ ilaenv( 1, 'ZUNMQR', ' ', n, 1, n, -1 ) ) )
470  IF( ilvsl ) THEN
471  maxwrk = max( maxwrk, n*( 1 +
472  $ ilaenv( 1, 'ZUNGQR', ' ', n, 1, n, -1 ) ) )
473  END IF
474  lwrk = maxwrk
475  IF( ijob.GE.1 )
476  $ lwrk = max( lwrk, n*n/2 )
477  ELSE
478  minwrk = 1
479  maxwrk = 1
480  lwrk = 1
481  END IF
482  work( 1 ) = lwrk
483  IF( wantsn .OR. n.EQ.0 ) THEN
484  liwmin = 1
485  ELSE
486  liwmin = n + 2
487  END IF
488  iwork( 1 ) = liwmin
489 *
490  IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
491  info = -21
492  ELSE IF( liwork.LT.liwmin .AND. .NOT.lquery) THEN
493  info = -24
494  END IF
495  END IF
496 *
497  IF( info.NE.0 ) THEN
498  CALL xerbla( 'ZGGESX', -info )
499  return
500  ELSE IF (lquery) THEN
501  return
502  END IF
503 *
504 * Quick return if possible
505 *
506  IF( n.EQ.0 ) THEN
507  sdim = 0
508  return
509  END IF
510 *
511 * Get machine constants
512 *
513  eps = dlamch( 'P' )
514  smlnum = dlamch( 'S' )
515  bignum = one / smlnum
516  CALL dlabad( smlnum, bignum )
517  smlnum = sqrt( smlnum ) / eps
518  bignum = one / smlnum
519 *
520 * Scale A if max element outside range [SMLNUM,BIGNUM]
521 *
522  anrm = zlange( 'M', n, n, a, lda, rwork )
523  ilascl = .false.
524  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
525  anrmto = smlnum
526  ilascl = .true.
527  ELSE IF( anrm.GT.bignum ) THEN
528  anrmto = bignum
529  ilascl = .true.
530  END IF
531  IF( ilascl )
532  $ CALL zlascl( 'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
533 *
534 * Scale B if max element outside range [SMLNUM,BIGNUM]
535 *
536  bnrm = zlange( 'M', n, n, b, ldb, rwork )
537  ilbscl = .false.
538  IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
539  bnrmto = smlnum
540  ilbscl = .true.
541  ELSE IF( bnrm.GT.bignum ) THEN
542  bnrmto = bignum
543  ilbscl = .true.
544  END IF
545  IF( ilbscl )
546  $ CALL zlascl( 'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
547 *
548 * Permute the matrix to make it more nearly triangular
549 * (Real Workspace: need 6*N)
550 *
551  ileft = 1
552  iright = n + 1
553  irwrk = iright + n
554  CALL zggbal( 'P', n, a, lda, b, ldb, ilo, ihi, rwork( ileft ),
555  $ rwork( iright ), rwork( irwrk ), ierr )
556 *
557 * Reduce B to triangular form (QR decomposition of B)
558 * (Complex Workspace: need N, prefer N*NB)
559 *
560  irows = ihi + 1 - ilo
561  icols = n + 1 - ilo
562  itau = 1
563  iwrk = itau + irows
564  CALL zgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
565  $ work( iwrk ), lwork+1-iwrk, ierr )
566 *
567 * Apply the unitary transformation to matrix A
568 * (Complex Workspace: need N, prefer N*NB)
569 *
570  CALL zunmqr( 'L', 'C', irows, icols, irows, b( ilo, ilo ), ldb,
571  $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
572  $ lwork+1-iwrk, ierr )
573 *
574 * Initialize VSL
575 * (Complex Workspace: need N, prefer N*NB)
576 *
577  IF( ilvsl ) THEN
578  CALL zlaset( 'Full', n, n, czero, cone, vsl, ldvsl )
579  IF( irows.GT.1 ) THEN
580  CALL zlacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
581  $ vsl( ilo+1, ilo ), ldvsl )
582  END IF
583  CALL zungqr( irows, irows, irows, vsl( ilo, ilo ), ldvsl,
584  $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
585  END IF
586 *
587 * Initialize VSR
588 *
589  IF( ilvsr )
590  $ CALL zlaset( 'Full', n, n, czero, cone, vsr, ldvsr )
591 *
592 * Reduce to generalized Hessenberg form
593 * (Workspace: none needed)
594 *
595  CALL zgghrd( jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb, vsl,
596  $ ldvsl, vsr, ldvsr, ierr )
597 *
598  sdim = 0
599 *
600 * Perform QZ algorithm, computing Schur vectors if desired
601 * (Complex Workspace: need N)
602 * (Real Workspace: need N)
603 *
604  iwrk = itau
605  CALL zhgeqz( 'S', jobvsl, jobvsr, n, ilo, ihi, a, lda, b, ldb,
606  $ alpha, beta, vsl, ldvsl, vsr, ldvsr, work( iwrk ),
607  $ lwork+1-iwrk, rwork( irwrk ), ierr )
608  IF( ierr.NE.0 ) THEN
609  IF( ierr.GT.0 .AND. ierr.LE.n ) THEN
610  info = ierr
611  ELSE IF( ierr.GT.n .AND. ierr.LE.2*n ) THEN
612  info = ierr - n
613  ELSE
614  info = n + 1
615  END IF
616  go to 40
617  END IF
618 *
619 * Sort eigenvalues ALPHA/BETA and compute the reciprocal of
620 * condition number(s)
621 *
622  IF( wantst ) THEN
623 *
624 * Undo scaling on eigenvalues before SELCTGing
625 *
626  IF( ilascl )
627  $ CALL zlascl( 'G', 0, 0, anrmto, anrm, n, 1, alpha, n, ierr )
628  IF( ilbscl )
629  $ CALL zlascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
630 *
631 * Select eigenvalues
632 *
633  DO 10 i = 1, n
634  bwork( i ) = selctg( alpha( i ), beta( i ) )
635  10 continue
636 *
637 * Reorder eigenvalues, transform Generalized Schur vectors, and
638 * compute reciprocal condition numbers
639 * (Complex Workspace: If IJOB >= 1, need MAX(1, 2*SDIM*(N-SDIM))
640 * otherwise, need 1 )
641 *
642  CALL ztgsen( ijob, ilvsl, ilvsr, bwork, n, a, lda, b, ldb,
643  $ alpha, beta, vsl, ldvsl, vsr, ldvsr, sdim, pl, pr,
644  $ dif, work( iwrk ), lwork-iwrk+1, iwork, liwork,
645  $ ierr )
646 *
647  IF( ijob.GE.1 )
648  $ maxwrk = max( maxwrk, 2*sdim*( n-sdim ) )
649  IF( ierr.EQ.-21 ) THEN
650 *
651 * not enough complex workspace
652 *
653  info = -21
654  ELSE
655  IF( ijob.EQ.1 .OR. ijob.EQ.4 ) THEN
656  rconde( 1 ) = pl
657  rconde( 2 ) = pr
658  END IF
659  IF( ijob.EQ.2 .OR. ijob.EQ.4 ) THEN
660  rcondv( 1 ) = dif( 1 )
661  rcondv( 2 ) = dif( 2 )
662  END IF
663  IF( ierr.EQ.1 )
664  $ info = n + 3
665  END IF
666 *
667  END IF
668 *
669 * Apply permutation to VSL and VSR
670 * (Workspace: none needed)
671 *
672  IF( ilvsl )
673  $ CALL zggbak( 'P', 'L', n, ilo, ihi, rwork( ileft ),
674  $ rwork( iright ), n, vsl, ldvsl, ierr )
675 *
676  IF( ilvsr )
677  $ CALL zggbak( 'P', 'R', n, ilo, ihi, rwork( ileft ),
678  $ rwork( iright ), n, vsr, ldvsr, ierr )
679 *
680 * Undo scaling
681 *
682  IF( ilascl ) THEN
683  CALL zlascl( 'U', 0, 0, anrmto, anrm, n, n, a, lda, ierr )
684  CALL zlascl( 'G', 0, 0, anrmto, anrm, n, 1, alpha, n, ierr )
685  END IF
686 *
687  IF( ilbscl ) THEN
688  CALL zlascl( 'U', 0, 0, bnrmto, bnrm, n, n, b, ldb, ierr )
689  CALL zlascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
690  END IF
691 *
692  IF( wantst ) THEN
693 *
694 * Check if reordering is correct
695 *
696  lastsl = .true.
697  sdim = 0
698  DO 30 i = 1, n
699  cursl = selctg( alpha( i ), beta( i ) )
700  IF( cursl )
701  $ sdim = sdim + 1
702  IF( cursl .AND. .NOT.lastsl )
703  $ info = n + 2
704  lastsl = cursl
705  30 continue
706 *
707  END IF
708 *
709  40 continue
710 *
711  work( 1 ) = maxwrk
712  iwork( 1 ) = liwmin
713 *
714  return
715 *
716 * End of ZGGESX
717 *
718  END