Lots of trailing whitespaces in the files of Syd. Cleaning this. No big deal.
[platform/upstream/lapack.git] / SRC / sgges3.f
1 *> \brief <b> SGGES3 computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE matrices (blocked algorithm)</b>
2 *
3 *  =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 *            http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download SGGES3 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgges3.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgges3.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgges3.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 *  Definition:
19 *  ===========
20 *
21 *       SUBROUTINE SGGES3( JOBVSL, JOBVSR, SORT, SELCTG, N, A, LDA, B,
22 *      $                   LDB, SDIM, ALPHAR, ALPHAI, BETA, VSL, LDVSL,
23 *      $                   VSR, 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 *       REAL               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 *> SGGES3 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 *> SGGEV 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 REAL 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 REAL 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 REAL 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 REAL array, dimension (N)
171 *> \endverbatim
172 *>
173 *> \param[out] ALPHAI
174 *> \verbatim
175 *>          ALPHAI is REAL array, dimension (N)
176 *> \endverbatim
177 *>
178 *> \param[out] BETA
179 *> \verbatim
180 *>          BETA is REAL 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 REAL 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 REAL 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 REAL 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 *>
238 *>          If LWORK = -1, then a workspace query is assumed; the routine
239 *>          only calculates the optimal size of the WORK array, returns
240 *>          this value as the first entry of the WORK array, and no error
241 *>          message related to LWORK is issued by XERBLA.
242 *> \endverbatim
243 *>
244 *> \param[out] BWORK
245 *> \verbatim
246 *>          BWORK is LOGICAL array, dimension (N)
247 *>          Not referenced if SORT = 'N'.
248 *> \endverbatim
249 *>
250 *> \param[out] INFO
251 *> \verbatim
252 *>          INFO is INTEGER
253 *>          = 0:  successful exit
254 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
255 *>          = 1,...,N:
256 *>                The QZ iteration failed.  (A,B) are not in Schur
257 *>                form, but ALPHAR(j), ALPHAI(j), and BETA(j) should
258 *>                be correct for j=INFO+1,...,N.
259 *>          > N:  =N+1: other than QZ iteration failed in SHGEQZ.
260 *>                =N+2: after reordering, roundoff changed values of
261 *>                      some complex eigenvalues so that leading
262 *>                      eigenvalues in the Generalized Schur form no
263 *>                      longer satisfy SELCTG=.TRUE.  This could also
264 *>                      be caused due to scaling.
265 *>                =N+3: reordering failed in STGSEN.
266 *> \endverbatim
267 *
268 *  Authors:
269 *  ========
270 *
271 *> \author Univ. of Tennessee
272 *> \author Univ. of California Berkeley
273 *> \author Univ. of Colorado Denver
274 *> \author NAG Ltd.
275 *
276 *> \date January 2015
277 *
278 *> \ingroup realGEeigen
279 *
280 *  =====================================================================
281       SUBROUTINE SGGES3( JOBVSL, JOBVSR, SORT, SELCTG, N, A, LDA, B,
282      $                   LDB, SDIM, ALPHAR, ALPHAI, BETA, VSL, LDVSL,
283      $                   VSR, LDVSR, WORK, LWORK, BWORK, INFO )
284 *
285 *  -- LAPACK driver routine (version 3.6.0) --
286 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
287 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
288 *     January 2015
289 *
290 *     .. Scalar Arguments ..
291       CHARACTER          JOBVSL, JOBVSR, SORT
292       INTEGER            INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N, SDIM
293 *     ..
294 *     .. Array Arguments ..
295       LOGICAL            BWORK( * )
296       REAL               A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
297      $                   B( LDB, * ), BETA( * ), VSL( LDVSL, * ),
298      $                   VSR( LDVSR, * ), WORK( * )
299 *     ..
300 *     .. Function Arguments ..
301       LOGICAL            SELCTG
302       EXTERNAL           SELCTG
303 *     ..
304 *
305 *  =====================================================================
306 *
307 *     .. Parameters ..
308       REAL               ZERO, ONE
309       PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
310 *     ..
311 *     .. Local Scalars ..
312       LOGICAL            CURSL, ILASCL, ILBSCL, ILVSL, ILVSR, LASTSL,
313      $                   LQUERY, LST2SL, WANTST
314       INTEGER            I, ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT,
315      $                   ILO, IP, IRIGHT, IROWS, ITAU, IWRK, LWKOPT
316       REAL               ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, PVSL,
317      $                   PVSR, SAFMAX, SAFMIN, SMLNUM
318 *     ..
319 *     .. Local Arrays ..
320       INTEGER            IDUM( 1 )
321       REAL               DIF( 2 )
322 *     ..
323 *     .. External Subroutines ..
324       EXTERNAL           SGEQRF, SGGBAK, SGGBAL, SGGHD3, SHGEQZ, SLABAD,
325      $                   SLACPY, SLASCL, SLASET, SORGQR, SORMQR, STGSEN,
326      $                   XERBLA
327 *     ..
328 *     .. External Functions ..
329       LOGICAL            LSAME
330       REAL               SLAMCH, SLANGE
331       EXTERNAL           LSAME, SLAMCH, SLANGE
332 *     ..
333 *     .. Intrinsic Functions ..
334       INTRINSIC          ABS, MAX, SQRT
335 *     ..
336 *     .. Executable Statements ..
337 *
338 *     Decode the input arguments
339 *
340       IF( LSAME( JOBVSL, 'N' ) ) THEN
341          IJOBVL = 1
342          ILVSL = .FALSE.
343       ELSE IF( LSAME( JOBVSL, 'V' ) ) THEN
344          IJOBVL = 2
345          ILVSL = .TRUE.
346       ELSE
347          IJOBVL = -1
348          ILVSL = .FALSE.
349       END IF
350 *
351       IF( LSAME( JOBVSR, 'N' ) ) THEN
352          IJOBVR = 1
353          ILVSR = .FALSE.
354       ELSE IF( LSAME( JOBVSR, 'V' ) ) THEN
355          IJOBVR = 2
356          ILVSR = .TRUE.
357       ELSE
358          IJOBVR = -1
359          ILVSR = .FALSE.
360       END IF
361 *
362       WANTST = LSAME( SORT, 'S' )
363 *
364 *     Test the input arguments
365 *
366       INFO = 0
367       LQUERY = ( LWORK.EQ.-1 )
368       IF( IJOBVL.LE.0 ) THEN
369          INFO = -1
370       ELSE IF( IJOBVR.LE.0 ) THEN
371          INFO = -2
372       ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN
373          INFO = -3
374       ELSE IF( N.LT.0 ) THEN
375          INFO = -5
376       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
377          INFO = -7
378       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
379          INFO = -9
380       ELSE IF( LDVSL.LT.1 .OR. ( ILVSL .AND. LDVSL.LT.N ) ) THEN
381          INFO = -15
382       ELSE IF( LDVSR.LT.1 .OR. ( ILVSR .AND. LDVSR.LT.N ) ) THEN
383          INFO = -17
384       ELSE IF( LWORK.LT.6*N+16 .AND. .NOT.LQUERY ) THEN
385          INFO = -19
386       END IF
387 *
388 *     Compute workspace
389 *
390       IF( INFO.EQ.0 ) THEN
391          CALL SGEQRF( N, N, B, LDB, WORK, WORK, -1, IERR )
392          LWKOPT = MAX( 6*N+16, 3*N+INT( WORK( 1 ) ) )
393          CALL SORMQR( 'L', 'T', N, N, N, B, LDB, WORK, A, LDA, WORK,
394      $                -1, IERR )
395          LWKOPT = MAX( LWKOPT, 3*N+INT( WORK( 1 ) ) )
396          IF( ILVSL ) THEN
397             CALL SORGQR( N, N, N, VSL, LDVSL, WORK, WORK, -1, IERR )
398             LWKOPT = MAX( LWKOPT, 3*N+INT( WORK( 1 ) ) )
399          END IF
400          CALL SGGHD3( JOBVSL, JOBVSR, N, 1, N, A, LDA, B, LDB, VSL,
401      $                LDVSL, VSR, LDVSR, WORK, -1, IERR )
402          LWKOPT = MAX( LWKOPT, 3*N+INT( WORK( 1 ) ) )
403          CALL SHGEQZ( 'S', JOBVSL, JOBVSR, N, 1, N, A, LDA, B, LDB,
404      $                ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR,
405      $                WORK, -1, IERR )
406          LWKOPT = MAX( LWKOPT, 2*N+INT( WORK( 1 ) ) )
407          IF( WANTST ) THEN
408             CALL STGSEN( 0, ILVSL, ILVSR, BWORK, N, A, LDA, B, LDB,
409      $                   ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR,
410      $                   SDIM, PVSL, PVSR, DIF, WORK, -1, IDUM, 1,
411      $                   IERR )
412             LWKOPT = MAX( LWKOPT, 2*N+INT( WORK( 1 ) ) )
413          END IF
414          WORK( 1 ) = LWKOPT
415       END IF
416 *
417       IF( INFO.NE.0 ) THEN
418          CALL XERBLA( 'SGGES3 ', -INFO )
419          RETURN
420       ELSE IF( LQUERY ) THEN
421          RETURN
422       END IF
423 *
424 *     Quick return if possible
425 *
426       IF( N.EQ.0 ) THEN
427          SDIM = 0
428          RETURN
429       END IF
430 *
431 *     Get machine constants
432 *
433       EPS = SLAMCH( 'P' )
434       SAFMIN = SLAMCH( 'S' )
435       SAFMAX = ONE / SAFMIN
436       CALL SLABAD( SAFMIN, SAFMAX )
437       SMLNUM = SQRT( SAFMIN ) / EPS
438       BIGNUM = ONE / SMLNUM
439 *
440 *     Scale A if max element outside range [SMLNUM,BIGNUM]
441 *
442       ANRM = SLANGE( 'M', N, N, A, LDA, WORK )
443       ILASCL = .FALSE.
444       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
445          ANRMTO = SMLNUM
446          ILASCL = .TRUE.
447       ELSE IF( ANRM.GT.BIGNUM ) THEN
448          ANRMTO = BIGNUM
449          ILASCL = .TRUE.
450       END IF
451       IF( ILASCL )
452      $   CALL SLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
453 *
454 *     Scale B if max element outside range [SMLNUM,BIGNUM]
455 *
456       BNRM = SLANGE( 'M', N, N, B, LDB, WORK )
457       ILBSCL = .FALSE.
458       IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
459          BNRMTO = SMLNUM
460          ILBSCL = .TRUE.
461       ELSE IF( BNRM.GT.BIGNUM ) THEN
462          BNRMTO = BIGNUM
463          ILBSCL = .TRUE.
464       END IF
465       IF( ILBSCL )
466      $   CALL SLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
467 *
468 *     Permute the matrix to make it more nearly triangular
469 *
470       ILEFT = 1
471       IRIGHT = N + 1
472       IWRK = IRIGHT + N
473       CALL SGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, WORK( ILEFT ),
474      $             WORK( IRIGHT ), WORK( IWRK ), IERR )
475 *
476 *     Reduce B to triangular form (QR decomposition of B)
477 *
478       IROWS = IHI + 1 - ILO
479       ICOLS = N + 1 - ILO
480       ITAU = IWRK
481       IWRK = ITAU + IROWS
482       CALL SGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
483      $             WORK( IWRK ), LWORK+1-IWRK, IERR )
484 *
485 *     Apply the orthogonal transformation to matrix A
486 *
487       CALL SORMQR( 'L', 'T', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
488      $             WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ),
489      $             LWORK+1-IWRK, IERR )
490 *
491 *     Initialize VSL
492 *
493       IF( ILVSL ) THEN
494          CALL SLASET( 'Full', N, N, ZERO, ONE, VSL, LDVSL )
495          IF( IROWS.GT.1 ) THEN
496             CALL SLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
497      $                   VSL( ILO+1, ILO ), LDVSL )
498          END IF
499          CALL SORGQR( IROWS, IROWS, IROWS, VSL( ILO, ILO ), LDVSL,
500      $                WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
501       END IF
502 *
503 *     Initialize VSR
504 *
505       IF( ILVSR )
506      $   CALL SLASET( 'Full', N, N, ZERO, ONE, VSR, LDVSR )
507 *
508 *     Reduce to generalized Hessenberg form
509 *
510       CALL SGGHD3( JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB, VSL,
511      $             LDVSL, VSR, LDVSR, WORK( IWRK ), LWORK+1-IWRK, IERR )
512 *
513 *     Perform QZ algorithm, computing Schur vectors if desired
514 *
515       IWRK = ITAU
516       CALL SHGEQZ( 'S', JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB,
517      $             ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR,
518      $             WORK( IWRK ), LWORK+1-IWRK, IERR )
519       IF( IERR.NE.0 ) THEN
520          IF( IERR.GT.0 .AND. IERR.LE.N ) THEN
521             INFO = IERR
522          ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN
523             INFO = IERR - N
524          ELSE
525             INFO = N + 1
526          END IF
527          GO TO 40
528       END IF
529 *
530 *     Sort eigenvalues ALPHA/BETA if desired
531 *
532       SDIM = 0
533       IF( WANTST ) THEN
534 *
535 *        Undo scaling on eigenvalues before SELCTGing
536 *
537          IF( ILASCL ) THEN
538             CALL SLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAR, N,
539      $                   IERR )
540             CALL SLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAI, N,
541      $                   IERR )
542          END IF
543          IF( ILBSCL )
544      $      CALL SLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
545 *
546 *        Select eigenvalues
547 *
548          DO 10 I = 1, N
549             BWORK( I ) = SELCTG( ALPHAR( I ), ALPHAI( I ), BETA( I ) )
550    10    CONTINUE
551 *
552          CALL STGSEN( 0, ILVSL, ILVSR, BWORK, N, A, LDA, B, LDB, ALPHAR,
553      $                ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR, SDIM, PVSL,
554      $                PVSR, DIF, WORK( IWRK ), LWORK-IWRK+1, IDUM, 1,
555      $                IERR )
556          IF( IERR.EQ.1 )
557      $      INFO = N + 3
558 *
559       END IF
560 *
561 *     Apply back-permutation to VSL and VSR
562 *
563       IF( ILVSL )
564      $   CALL SGGBAK( 'P', 'L', N, ILO, IHI, WORK( ILEFT ),
565      $                WORK( IRIGHT ), N, VSL, LDVSL, IERR )
566 *
567       IF( ILVSR )
568      $   CALL SGGBAK( 'P', 'R', N, ILO, IHI, WORK( ILEFT ),
569      $                WORK( IRIGHT ), N, VSR, LDVSR, IERR )
570 *
571 *     Check if unscaling would cause over/underflow, if so, rescale
572 *     (ALPHAR(I),ALPHAI(I),BETA(I)) so BETA(I) is on the order of
573 *     B(I,I) and ALPHAR(I) and ALPHAI(I) are on the order of A(I,I)
574 *
575       IF( ILASCL )THEN
576          DO 50 I = 1, N
577             IF( ALPHAI( I ).NE.ZERO ) THEN
578                IF( ( ALPHAR( I )/SAFMAX ).GT.( ANRMTO/ANRM ) .OR.
579      $             ( SAFMIN/ALPHAR( I ) ).GT.( ANRM/ANRMTO ) ) THEN
580                   WORK( 1 ) = ABS( A( I, I )/ALPHAR( I ) )
581                   BETA( I ) = BETA( I )*WORK( 1 )
582                   ALPHAR( I ) = ALPHAR( I )*WORK( 1 )
583                   ALPHAI( I ) = ALPHAI( I )*WORK( 1 )
584                ELSE IF( ( ALPHAI( I )/SAFMAX ).GT.( ANRMTO/ANRM ) .OR.
585      $             ( SAFMIN/ALPHAI( I ) ).GT.( ANRM/ANRMTO ) ) THEN
586                   WORK( 1 ) = ABS( A( I, I+1 )/ALPHAI( I ) )
587                   BETA( I ) = BETA( I )*WORK( 1 )
588                   ALPHAR( I ) = ALPHAR( I )*WORK( 1 )
589                   ALPHAI( I ) = ALPHAI( I )*WORK( 1 )
590                END IF
591             END IF
592    50    CONTINUE
593       END IF
594 *
595       IF( ILBSCL )THEN
596          DO 60 I = 1, N
597             IF( ALPHAI( I ).NE.ZERO ) THEN
598                 IF( ( BETA( I )/SAFMAX ).GT.( BNRMTO/BNRM ) .OR.
599      $              ( SAFMIN/BETA( I ) ).GT.( BNRM/BNRMTO ) ) THEN
600                    WORK( 1 ) = ABS(B( I, I )/BETA( I ))
601                    BETA( I ) = BETA( I )*WORK( 1 )
602                    ALPHAR( I ) = ALPHAR( I )*WORK( 1 )
603                    ALPHAI( I ) = ALPHAI( I )*WORK( 1 )
604                 END IF
605              END IF
606    60    CONTINUE
607       END IF
608 *
609 *     Undo scaling
610 *
611       IF( ILASCL ) THEN
612          CALL SLASCL( 'H', 0, 0, ANRMTO, ANRM, N, N, A, LDA, IERR )
613          CALL SLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAR, N, IERR )
614          CALL SLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAI, N, IERR )
615       END IF
616 *
617       IF( ILBSCL ) THEN
618          CALL SLASCL( 'U', 0, 0, BNRMTO, BNRM, N, N, B, LDB, IERR )
619          CALL SLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
620       END IF
621 *
622       IF( WANTST ) THEN
623 *
624 *        Check if reordering is correct
625 *
626          LASTSL = .TRUE.
627          LST2SL = .TRUE.
628          SDIM = 0
629          IP = 0
630          DO 30 I = 1, N
631             CURSL = SELCTG( ALPHAR( I ), ALPHAI( I ), BETA( I ) )
632             IF( ALPHAI( I ).EQ.ZERO ) THEN
633                IF( CURSL )
634      $            SDIM = SDIM + 1
635                IP = 0
636                IF( CURSL .AND. .NOT.LASTSL )
637      $            INFO = N + 2
638             ELSE
639                IF( IP.EQ.1 ) THEN
640 *
641 *                 Last eigenvalue of conjugate pair
642 *
643                   CURSL = CURSL .OR. LASTSL
644                   LASTSL = CURSL
645                   IF( CURSL )
646      $               SDIM = SDIM + 2
647                   IP = -1
648                   IF( CURSL .AND. .NOT.LST2SL )
649      $               INFO = N + 2
650                ELSE
651 *
652 *                 First eigenvalue of conjugate pair
653 *
654                   IP = 1
655                END IF
656             END IF
657             LST2SL = LASTSL
658             LASTSL = CURSL
659    30    CONTINUE
660 *
661       END IF
662 *
663    40 CONTINUE
664 *
665       WORK( 1 ) = LWKOPT
666 *
667       RETURN
668 *
669 *     End of SGGES3
670 *
671       END