LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
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
subroutine zlacpy(UPLO, M, N, A, LDA, B, LDB)
ZLACPY copies all or part of one two-dimensional array to another.
Definition: zlacpy.f:105
subroutine zgghrd(COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, LDQ, Z, LDZ, INFO)
ZGGHRD
Definition: zgghrd.f:206
subroutine zgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
ZGEQRF VARIANT: left-looking Level 3 BLAS of the algorithm.
Definition: zgeqrf.f:151
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:145
subroutine zggbak(JOB, SIDE, N, ILO, IHI, LSCALE, RSCALE, M, V, LDV, INFO)
ZGGBAK
Definition: zggbak.f:150
subroutine zlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
ZLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: zlaset.f:108
subroutine zggbal(JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE, WORK, INFO)
ZGGBAL
Definition: zggbal.f:179
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:76
subroutine zggevx(BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, B, LDB, ALPHA, BETA, VL, LDVL, VR, LDVR, ILO, IHI, LSCALE, RSCALE, ABNRM, BBNRM, RCONDE, RCONDV, WORK, LWORK, RWORK, IWORK, BWORK, INFO)
ZGGEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices ...
Definition: zggevx.f:376
subroutine zunmqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
ZUNMQR
Definition: zunmqr.f:169
subroutine zungqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
ZUNGQR
Definition: zungqr.f:130
subroutine ztgevc(SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL, LDVL, VR, LDVR, MM, M, WORK, RWORK, INFO)
ZTGEVC
Definition: ztgevc.f:221
subroutine zhgeqz(JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT, ALPHA, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, RWORK, INFO)
ZHGEQZ
Definition: zhgeqz.f:286
subroutine zlascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
ZLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: zlascl.f:145
subroutine ztgsna(JOB, HOWMNY, SELECT, N, A, LDA, B, LDB, VL, LDVL, VR, LDVR, S, DIF, MM, M, WORK, LWORK, IWORK, INFO)
ZTGSNA
Definition: ztgsna.f:313