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