LAPACK  3.5.0
LAPACK: Linear Algebra PACKage
 All Classes Files Functions Variables Typedefs Macros
sggevx.f
Go to the documentation of this file.
1 *> \brief <b> SGGEVX 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 SGGEVX + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sggevx.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sggevx.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sggevx.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE SGGEVX( 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 * REAL ABNRM, BBNRM
30 * ..
31 * .. Array Arguments ..
32 * LOGICAL BWORK( * )
33 * INTEGER IWORK( * )
34 * REAL 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 *> SGGEVX 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 REAL 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 REAL 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 REAL array, dimension (N)
156 *> \endverbatim
157 *>
158 *> \param[out] ALPHAI
159 *> \verbatim
160 *> ALPHAI is REAL array, dimension (N)
161 *> \endverbatim
162 *>
163 *> \param[out] BETA
164 *> \verbatim
165 *> BETA is REAL 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 REAL 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 REAL 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 REAL 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 REAL 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 REAL
267 *> The one-norm of the balanced matrix A.
268 *> \endverbatim
269 *>
270 *> \param[out] BBNRM
271 *> \verbatim
272 *> BBNRM is REAL
273 *> The one-norm of the balanced matrix B.
274 *> \endverbatim
275 *>
276 *> \param[out] RCONDE
277 *> \verbatim
278 *> RCONDE is REAL 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 REAL 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 REAL 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', 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 SHGEQZ.
344 *> =N+2: error return from STGEVC.
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 realGEeigen
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 sggevx( 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  REAL abnrm, bbnrm
403 * ..
404 * .. Array Arguments ..
405  LOGICAL bwork( * )
406  INTEGER iwork( * )
407  REAL 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  REAL zero, one
417  parameter( zero = 0.0e+0, one = 1.0e+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  REAL anrm, anrmto, bignum, bnrm, bnrmto, eps,
427  $ smlnum, temp
428 * ..
429 * .. Local Arrays ..
430  LOGICAL ldumma( 1 )
431 * ..
432 * .. External Subroutines ..
433  EXTERNAL sgeqrf, sggbak, sggbal, sgghrd, shgeqz, slabad,
435  $ stgsna, xerbla
436 * ..
437 * .. External Functions ..
438  LOGICAL lsame
439  INTEGER ilaenv
440  REAL slamch, slange
441  EXTERNAL lsame, ilaenv, slamch, slange
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.( noscl .OR. lsame( balanc, 'S' ) .OR.
484  $ lsame( balanc, 'B' ) ) ) THEN
485  info = -1
486  ELSE IF( ijobvl.LE.0 ) THEN
487  info = -2
488  ELSE IF( ijobvr.LE.0 ) THEN
489  info = -3
490  ELSE IF( .NOT.( wantsn .OR. wantse .OR. wantsb .OR. wantsv ) )
491  $ THEN
492  info = -4
493  ELSE IF( n.LT.0 ) THEN
494  info = -5
495  ELSE IF( lda.LT.max( 1, n ) ) THEN
496  info = -7
497  ELSE IF( ldb.LT.max( 1, n ) ) THEN
498  info = -9
499  ELSE IF( ldvl.LT.1 .OR. ( ilvl .AND. ldvl.LT.n ) ) THEN
500  info = -14
501  ELSE IF( ldvr.LT.1 .OR. ( ilvr .AND. ldvr.LT.n ) ) THEN
502  info = -16
503  END IF
504 *
505 * Compute workspace
506 * (Note: Comments in the code beginning "Workspace:" describe the
507 * minimal amount of workspace needed at that point in the code,
508 * as well as the preferred amount for good performance.
509 * NB refers to the optimal block size for the immediately
510 * following subroutine, as returned by ILAENV. The workspace is
511 * computed assuming ILO = 1 and IHI = N, the worst case.)
512 *
513  IF( info.EQ.0 ) THEN
514  IF( n.EQ.0 ) THEN
515  minwrk = 1
516  maxwrk = 1
517  ELSE
518  IF( noscl .AND. .NOT.ilv ) THEN
519  minwrk = 2*n
520  ELSE
521  minwrk = 6*n
522  END IF
523  IF( wantse ) THEN
524  minwrk = 10*n
525  ELSE IF( wantsv .OR. wantsb ) THEN
526  minwrk = 2*n*( n + 4 ) + 16
527  END IF
528  maxwrk = minwrk
529  maxwrk = max( maxwrk,
530  $ n + n*ilaenv( 1, 'SGEQRF', ' ', n, 1, n, 0 ) )
531  maxwrk = max( maxwrk,
532  $ n + n*ilaenv( 1, 'SORMQR', ' ', n, 1, n, 0 ) )
533  IF( ilvl ) THEN
534  maxwrk = max( maxwrk, n +
535  $ n*ilaenv( 1, 'SORGQR', ' ', n, 1, n, 0 ) )
536  END IF
537  END IF
538  work( 1 ) = maxwrk
539 *
540  IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
541  info = -26
542  END IF
543  END IF
544 *
545  IF( info.NE.0 ) THEN
546  CALL xerbla( 'SGGEVX', -info )
547  RETURN
548  ELSE IF( lquery ) THEN
549  RETURN
550  END IF
551 *
552 * Quick return if possible
553 *
554  IF( n.EQ.0 )
555  $ RETURN
556 *
557 *
558 * Get machine constants
559 *
560  eps = slamch( 'P' )
561  smlnum = slamch( 'S' )
562  bignum = one / smlnum
563  CALL slabad( smlnum, bignum )
564  smlnum = sqrt( smlnum ) / eps
565  bignum = one / smlnum
566 *
567 * Scale A if max element outside range [SMLNUM,BIGNUM]
568 *
569  anrm = slange( 'M', n, n, a, lda, work )
570  ilascl = .false.
571  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
572  anrmto = smlnum
573  ilascl = .true.
574  ELSE IF( anrm.GT.bignum ) THEN
575  anrmto = bignum
576  ilascl = .true.
577  END IF
578  IF( ilascl )
579  $ CALL slascl( 'G', 0, 0, anrm, anrmto, n, n, a, lda, ierr )
580 *
581 * Scale B if max element outside range [SMLNUM,BIGNUM]
582 *
583  bnrm = slange( 'M', n, n, b, ldb, work )
584  ilbscl = .false.
585  IF( bnrm.GT.zero .AND. bnrm.LT.smlnum ) THEN
586  bnrmto = smlnum
587  ilbscl = .true.
588  ELSE IF( bnrm.GT.bignum ) THEN
589  bnrmto = bignum
590  ilbscl = .true.
591  END IF
592  IF( ilbscl )
593  $ CALL slascl( 'G', 0, 0, bnrm, bnrmto, n, n, b, ldb, ierr )
594 *
595 * Permute and/or balance the matrix pair (A,B)
596 * (Workspace: need 6*N if BALANC = 'S' or 'B', 1 otherwise)
597 *
598  CALL sggbal( balanc, n, a, lda, b, ldb, ilo, ihi, lscale, rscale,
599  $ work, ierr )
600 *
601 * Compute ABNRM and BBNRM
602 *
603  abnrm = slange( '1', n, n, a, lda, work( 1 ) )
604  IF( ilascl ) THEN
605  work( 1 ) = abnrm
606  CALL slascl( 'G', 0, 0, anrmto, anrm, 1, 1, work( 1 ), 1,
607  $ ierr )
608  abnrm = work( 1 )
609  END IF
610 *
611  bbnrm = slange( '1', n, n, b, ldb, work( 1 ) )
612  IF( ilbscl ) THEN
613  work( 1 ) = bbnrm
614  CALL slascl( 'G', 0, 0, bnrmto, bnrm, 1, 1, work( 1 ), 1,
615  $ ierr )
616  bbnrm = work( 1 )
617  END IF
618 *
619 * Reduce B to triangular form (QR decomposition of B)
620 * (Workspace: need N, prefer N*NB )
621 *
622  irows = ihi + 1 - ilo
623  IF( ilv .OR. .NOT.wantsn ) THEN
624  icols = n + 1 - ilo
625  ELSE
626  icols = irows
627  END IF
628  itau = 1
629  iwrk = itau + irows
630  CALL sgeqrf( irows, icols, b( ilo, ilo ), ldb, work( itau ),
631  $ work( iwrk ), lwork+1-iwrk, ierr )
632 *
633 * Apply the orthogonal transformation to A
634 * (Workspace: need N, prefer N*NB)
635 *
636  CALL sormqr( 'L', 'T', irows, icols, irows, b( ilo, ilo ), ldb,
637  $ work( itau ), a( ilo, ilo ), lda, work( iwrk ),
638  $ lwork+1-iwrk, ierr )
639 *
640 * Initialize VL and/or VR
641 * (Workspace: need N, prefer N*NB)
642 *
643  IF( ilvl ) THEN
644  CALL slaset( 'Full', n, n, zero, one, vl, ldvl )
645  IF( irows.GT.1 ) THEN
646  CALL slacpy( 'L', irows-1, irows-1, b( ilo+1, ilo ), ldb,
647  $ vl( ilo+1, ilo ), ldvl )
648  END IF
649  CALL sorgqr( irows, irows, irows, vl( ilo, ilo ), ldvl,
650  $ work( itau ), work( iwrk ), lwork+1-iwrk, ierr )
651  END IF
652 *
653  IF( ilvr )
654  $ CALL slaset( 'Full', n, n, zero, one, vr, ldvr )
655 *
656 * Reduce to generalized Hessenberg form
657 * (Workspace: none needed)
658 *
659  IF( ilv .OR. .NOT.wantsn ) THEN
660 *
661 * Eigenvectors requested -- work on whole matrix.
662 *
663  CALL sgghrd( jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb, vl,
664  $ ldvl, vr, ldvr, ierr )
665  ELSE
666  CALL sgghrd( 'N', 'N', irows, 1, irows, a( ilo, ilo ), lda,
667  $ b( ilo, ilo ), ldb, vl, ldvl, vr, ldvr, ierr )
668  END IF
669 *
670 * Perform QZ algorithm (Compute eigenvalues, and optionally, the
671 * Schur forms and Schur vectors)
672 * (Workspace: need N)
673 *
674  IF( ilv .OR. .NOT.wantsn ) THEN
675  chtemp = 'S'
676  ELSE
677  chtemp = 'E'
678  END IF
679 *
680  CALL shgeqz( chtemp, jobvl, jobvr, n, ilo, ihi, a, lda, b, ldb,
681  $ alphar, alphai, beta, vl, ldvl, vr, ldvr, work,
682  $ lwork, ierr )
683  IF( ierr.NE.0 ) THEN
684  IF( ierr.GT.0 .AND. ierr.LE.n ) THEN
685  info = ierr
686  ELSE IF( ierr.GT.n .AND. ierr.LE.2*n ) THEN
687  info = ierr - n
688  ELSE
689  info = n + 1
690  END IF
691  go to 130
692  END IF
693 *
694 * Compute Eigenvectors and estimate condition numbers if desired
695 * (Workspace: STGEVC: need 6*N
696 * STGSNA: need 2*N*(N+2)+16 if SENSE = 'V' or 'B',
697 * need N otherwise )
698 *
699  IF( ilv .OR. .NOT.wantsn ) THEN
700  IF( ilv ) THEN
701  IF( ilvl ) THEN
702  IF( ilvr ) THEN
703  chtemp = 'B'
704  ELSE
705  chtemp = 'L'
706  END IF
707  ELSE
708  chtemp = 'R'
709  END IF
710 *
711  CALL stgevc( chtemp, 'B', ldumma, n, a, lda, b, ldb, vl,
712  $ ldvl, vr, ldvr, n, in, work, ierr )
713  IF( ierr.NE.0 ) THEN
714  info = n + 2
715  go to 130
716  END IF
717  END IF
718 *
719  IF( .NOT.wantsn ) THEN
720 *
721 * compute eigenvectors (STGEVC) and estimate condition
722 * numbers (STGSNA). Note that the definition of the condition
723 * number is not invariant under transformation (u,v) to
724 * (Q*u, Z*v), where (u,v) are eigenvectors of the generalized
725 * Schur form (S,T), Q and Z are orthogonal matrices. In order
726 * to avoid using extra 2*N*N workspace, we have to recalculate
727 * eigenvectors and estimate one condition numbers at a time.
728 *
729  pair = .false.
730  DO 20 i = 1, n
731 *
732  IF( pair ) THEN
733  pair = .false.
734  go to 20
735  END IF
736  mm = 1
737  IF( i.LT.n ) THEN
738  IF( a( i+1, i ).NE.zero ) THEN
739  pair = .true.
740  mm = 2
741  END IF
742  END IF
743 *
744  DO 10 j = 1, n
745  bwork( j ) = .false.
746  10 CONTINUE
747  IF( mm.EQ.1 ) THEN
748  bwork( i ) = .true.
749  ELSE IF( mm.EQ.2 ) THEN
750  bwork( i ) = .true.
751  bwork( i+1 ) = .true.
752  END IF
753 *
754  iwrk = mm*n + 1
755  iwrk1 = iwrk + mm*n
756 *
757 * Compute a pair of left and right eigenvectors.
758 * (compute workspace: need up to 4*N + 6*N)
759 *
760  IF( wantse .OR. wantsb ) THEN
761  CALL stgevc( 'B', 'S', bwork, n, a, lda, b, ldb,
762  $ work( 1 ), n, work( iwrk ), n, mm, m,
763  $ work( iwrk1 ), ierr )
764  IF( ierr.NE.0 ) THEN
765  info = n + 2
766  go to 130
767  END IF
768  END IF
769 *
770  CALL stgsna( sense, 'S', bwork, n, a, lda, b, ldb,
771  $ work( 1 ), n, work( iwrk ), n, rconde( i ),
772  $ rcondv( i ), mm, m, work( iwrk1 ),
773  $ lwork-iwrk1+1, iwork, ierr )
774 *
775  20 CONTINUE
776  END IF
777  END IF
778 *
779 * Undo balancing on VL and VR and normalization
780 * (Workspace: none needed)
781 *
782  IF( ilvl ) THEN
783  CALL sggbak( balanc, 'L', n, ilo, ihi, lscale, rscale, n, vl,
784  $ ldvl, ierr )
785 *
786  DO 70 jc = 1, n
787  IF( alphai( jc ).LT.zero )
788  $ go to 70
789  temp = zero
790  IF( alphai( jc ).EQ.zero ) THEN
791  DO 30 jr = 1, n
792  temp = max( temp, abs( vl( jr, jc ) ) )
793  30 CONTINUE
794  ELSE
795  DO 40 jr = 1, n
796  temp = max( temp, abs( vl( jr, jc ) )+
797  $ abs( vl( jr, jc+1 ) ) )
798  40 CONTINUE
799  END IF
800  IF( temp.LT.smlnum )
801  $ go to 70
802  temp = one / temp
803  IF( alphai( jc ).EQ.zero ) THEN
804  DO 50 jr = 1, n
805  vl( jr, jc ) = vl( jr, jc )*temp
806  50 CONTINUE
807  ELSE
808  DO 60 jr = 1, n
809  vl( jr, jc ) = vl( jr, jc )*temp
810  vl( jr, jc+1 ) = vl( jr, jc+1 )*temp
811  60 CONTINUE
812  END IF
813  70 CONTINUE
814  END IF
815  IF( ilvr ) THEN
816  CALL sggbak( balanc, 'R', n, ilo, ihi, lscale, rscale, n, vr,
817  $ ldvr, ierr )
818  DO 120 jc = 1, n
819  IF( alphai( jc ).LT.zero )
820  $ go to 120
821  temp = zero
822  IF( alphai( jc ).EQ.zero ) THEN
823  DO 80 jr = 1, n
824  temp = max( temp, abs( vr( jr, jc ) ) )
825  80 CONTINUE
826  ELSE
827  DO 90 jr = 1, n
828  temp = max( temp, abs( vr( jr, jc ) )+
829  $ abs( vr( jr, jc+1 ) ) )
830  90 CONTINUE
831  END IF
832  IF( temp.LT.smlnum )
833  $ go to 120
834  temp = one / temp
835  IF( alphai( jc ).EQ.zero ) THEN
836  DO 100 jr = 1, n
837  vr( jr, jc ) = vr( jr, jc )*temp
838  100 CONTINUE
839  ELSE
840  DO 110 jr = 1, n
841  vr( jr, jc ) = vr( jr, jc )*temp
842  vr( jr, jc+1 ) = vr( jr, jc+1 )*temp
843  110 CONTINUE
844  END IF
845  120 CONTINUE
846  END IF
847 *
848 * Undo scaling if necessary
849 *
850  130 CONTINUE
851 *
852  IF( ilascl ) THEN
853  CALL slascl( 'G', 0, 0, anrmto, anrm, n, 1, alphar, n, ierr )
854  CALL slascl( 'G', 0, 0, anrmto, anrm, n, 1, alphai, n, ierr )
855  END IF
856 *
857  IF( ilbscl ) THEN
858  CALL slascl( 'G', 0, 0, bnrmto, bnrm, n, 1, beta, n, ierr )
859  END IF
860 *
861  work( 1 ) = maxwrk
862  RETURN
863 *
864 * End of SGGEVX
865 *
866  END