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