LAPACK  3.6.0
LAPACK: Linear Algebra PACKage
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
subroutine stgevc(SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL, LDVL, VR, LDVR, MM, M, WORK, INFO)
STGEVC
Definition: stgevc.f:297
subroutine sormqr(SIDE, TRANS, M, N, K, A, LDA, TAU, C, LDC, WORK, LWORK, INFO)
SORMQR
Definition: sormqr.f:170
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine sgghrd(COMPQ, COMPZ, N, ILO, IHI, A, LDA, B, LDB, Q, LDQ, Z, LDZ, INFO)
SGGHRD
Definition: sgghrd.f:209
subroutine slabad(SMALL, LARGE)
SLABAD
Definition: slabad.f:76
subroutine sgeqrf(M, N, A, LDA, TAU, WORK, LWORK, INFO)
SGEQRF
Definition: sgeqrf.f:138
subroutine sggbal(JOB, N, A, LDA, B, LDB, ILO, IHI, LSCALE, RSCALE, WORK, INFO)
SGGBAL
Definition: sggbal.f:179
subroutine sggevx(BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, B, LDB, ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR, ILO, IHI, LSCALE, RSCALE, ABNRM, BBNRM, RCONDE, RCONDV, WORK, LWORK, IWORK, BWORK, INFO)
SGGEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices ...
Definition: sggevx.f:393
subroutine shgeqz(JOB, COMPQ, COMPZ, N, ILO, IHI, H, LDH, T, LDT, ALPHAR, ALPHAI, BETA, Q, LDQ, Z, LDZ, WORK, LWORK, INFO)
SHGEQZ
Definition: shgeqz.f:306
subroutine slacpy(UPLO, M, N, A, LDA, B, LDB)
SLACPY copies all or part of one two-dimensional array to another.
Definition: slacpy.f:105
subroutine stgsna(JOB, HOWMNY, SELECT, N, A, LDA, B, LDB, VL, LDVL, VR, LDVR, S, DIF, MM, M, WORK, LWORK, IWORK, INFO)
STGSNA
Definition: stgsna.f:383
subroutine slascl(TYPE, KL, KU, CFROM, CTO, M, N, A, LDA, INFO)
SLASCL multiplies a general rectangular matrix by a real scalar defined as cto/cfrom.
Definition: slascl.f:141
subroutine sorgqr(M, N, K, A, LDA, TAU, WORK, LWORK, INFO)
SORGQR
Definition: sorgqr.f:130
subroutine sggbak(JOB, SIDE, N, ILO, IHI, LSCALE, RSCALE, M, V, LDV, INFO)
SGGBAK
Definition: sggbak.f:149
subroutine slaset(UPLO, M, N, ALPHA, BETA, A, LDA)
SLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: slaset.f:112