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