1*> \brief <b> CGEES 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 CGEES + dependencies
10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgees.f">
11*> [TGZ]</a>
12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgees.f">
13*> [ZIP]</a>
14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgees.f">
15*> [TXT]</a>
16*> \endhtmlonly
17*
18*  Definition:
19*  ===========
20*
21*       SUBROUTINE CGEES( JOBVS, SORT, SELECT, N, A, LDA, SDIM, W, VS,
22*                         LDVS, WORK, LWORK, RWORK, BWORK, INFO )
23*
24*       .. Scalar Arguments ..
25*       CHARACTER          JOBVS, SORT
26*       INTEGER            INFO, LDA, LDVS, LWORK, N, SDIM
27*       ..
28*       .. Array Arguments ..
29*       LOGICAL            BWORK( * )
30*       REAL               RWORK( * )
31*       COMPLEX            A( LDA, * ), VS( LDVS, * ), W( * ), WORK( * )
32*       ..
33*       .. Function Arguments ..
34*       LOGICAL            SELECT
35*       EXTERNAL           SELECT
36*       ..
37*
38*
39*> \par Purpose:
40*  =============
41*>
42*> \verbatim
43*>
44*> CGEES computes for an N-by-N complex nonsymmetric matrix A, the
45*> eigenvalues, the Schur form T, and, optionally, the matrix of Schur
46*> vectors Z.  This gives the Schur factorization A = Z*T*(Z**H).
47*>
48*> Optionally, it also orders the eigenvalues on the diagonal of the
49*> Schur form so that selected eigenvalues are at the top left.
50*> The leading columns of Z then form an orthonormal basis for the
51*> invariant subspace corresponding to the selected eigenvalues.
52*>
53*> A complex matrix is in Schur form if it is upper triangular.
54*> \endverbatim
55*
56*  Arguments:
57*  ==========
58*
59*> \param[in] JOBVS
60*> \verbatim
61*>          JOBVS is CHARACTER*1
62*>          = 'N': Schur vectors are not computed;
63*>          = 'V': Schur vectors are computed.
64*> \endverbatim
65*>
66*> \param[in] SORT
67*> \verbatim
68*>          SORT is CHARACTER*1
69*>          Specifies whether or not to order the eigenvalues on the
70*>          diagonal of the Schur form.
71*>          = 'N': Eigenvalues are not ordered:
72*>          = 'S': Eigenvalues are ordered (see SELECT).
73*> \endverbatim
74*>
75*> \param[in] SELECT
76*> \verbatim
77*>          SELECT is a LOGICAL FUNCTION of one COMPLEX argument
78*>          SELECT must be declared EXTERNAL in the calling subroutine.
79*>          If SORT = 'S', SELECT is used to select eigenvalues to order
80*>          to the top left of the Schur form.
81*>          IF SORT = 'N', SELECT is not referenced.
82*>          The eigenvalue W(j) is selected if SELECT(W(j)) is true.
83*> \endverbatim
84*>
85*> \param[in] N
86*> \verbatim
87*>          N is INTEGER
88*>          The order of the matrix A. N >= 0.
89*> \endverbatim
90*>
91*> \param[in,out] A
92*> \verbatim
93*>          A is COMPLEX array, dimension (LDA,N)
94*>          On entry, the N-by-N matrix A.
95*>          On exit, A has been overwritten by its Schur form T.
96*> \endverbatim
97*>
98*> \param[in] LDA
99*> \verbatim
100*>          LDA is INTEGER
101*>          The leading dimension of the array A.  LDA >= max(1,N).
102*> \endverbatim
103*>
104*> \param[out] SDIM
105*> \verbatim
106*>          SDIM is INTEGER
107*>          If SORT = 'N', SDIM = 0.
108*>          If SORT = 'S', SDIM = number of eigenvalues for which
109*>                         SELECT is true.
110*> \endverbatim
111*>
112*> \param[out] W
113*> \verbatim
114*>          W is COMPLEX array, dimension (N)
115*>          W contains the computed eigenvalues, in the same order that
116*>          they appear on the diagonal of the output Schur form T.
117*> \endverbatim
118*>
119*> \param[out] VS
120*> \verbatim
121*>          VS is COMPLEX array, dimension (LDVS,N)
122*>          If JOBVS = 'V', VS contains the unitary matrix Z of Schur
123*>          vectors.
124*>          If JOBVS = 'N', VS is not referenced.
125*> \endverbatim
126*>
127*> \param[in] LDVS
128*> \verbatim
129*>          LDVS is INTEGER
130*>          The leading dimension of the array VS.  LDVS >= 1; if
131*>          JOBVS = 'V', LDVS >= N.
132*> \endverbatim
133*>
134*> \param[out] WORK
135*> \verbatim
136*>          WORK is COMPLEX array, dimension (MAX(1,LWORK))
137*>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
138*> \endverbatim
139*>
140*> \param[in] LWORK
141*> \verbatim
142*>          LWORK is INTEGER
143*>          The dimension of the array WORK.  LWORK >= max(1,2*N).
144*>          For good performance, LWORK must generally be larger.
145*>
146*>          If LWORK = -1, then a workspace query is assumed; the routine
147*>          only calculates the optimal size of the WORK array, returns
148*>          this value as the first entry of the WORK array, and no error
149*>          message related to LWORK is issued by XERBLA.
150*> \endverbatim
151*>
152*> \param[out] RWORK
153*> \verbatim
154*>          RWORK is REAL array, dimension (N)
155*> \endverbatim
156*>
157*> \param[out] BWORK
158*> \verbatim
159*>          BWORK is LOGICAL array, dimension (N)
160*>          Not referenced if SORT = 'N'.
161*> \endverbatim
162*>
163*> \param[out] INFO
164*> \verbatim
165*>          INFO is INTEGER
166*>          = 0: successful exit
167*>          < 0: if INFO = -i, the i-th argument had an illegal value.
168*>          > 0: if INFO = i, and i is
169*>               <= N:  the QR algorithm failed to compute all the
170*>                      eigenvalues; elements 1:ILO-1 and i+1:N of W
171*>                      contain those eigenvalues which have converged;
172*>                      if JOBVS = 'V', VS contains the matrix which
173*>                      reduces A to its partially converged Schur form.
174*>               = N+1: the eigenvalues could not be reordered because
175*>                      some eigenvalues were too close to separate (the
176*>                      problem is very ill-conditioned);
177*>               = N+2: after reordering, roundoff changed values of
178*>                      some complex eigenvalues so that leading
179*>                      eigenvalues in the Schur form no longer satisfy
180*>                      SELECT = .TRUE..  This could also be caused by
181*>                      underflow due to scaling.
182*> \endverbatim
183*
184*  Authors:
185*  ========
186*
187*> \author Univ. of Tennessee
188*> \author Univ. of California Berkeley
189*> \author Univ. of Colorado Denver
190*> \author NAG Ltd.
191*
192*> \date November 2011
193*
194*> \ingroup complexGEeigen
195*
196*  =====================================================================
197      SUBROUTINE CGEES( JOBVS, SORT, SELECT, N, A, LDA, SDIM, W, VS,
198     $                  LDVS, WORK, LWORK, RWORK, BWORK, INFO )
199*
200*  -- LAPACK driver routine (version 3.4.0) --
201*  -- LAPACK is a software package provided by Univ. of Tennessee,    --
202*  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
203*     November 2011
204*
205*     .. Scalar Arguments ..
206      CHARACTER          JOBVS, SORT
207      INTEGER            INFO, LDA, LDVS, LWORK, N, SDIM
208*     ..
209*     .. Array Arguments ..
210      LOGICAL            BWORK( * )
211      REAL               RWORK( * )
212      COMPLEX            A( LDA, * ), VS( LDVS, * ), W( * ), WORK( * )
213*     ..
214*     .. Function Arguments ..
215      LOGICAL            SELECT
216      EXTERNAL           SELECT
217*     ..
218*
219*  =====================================================================
220*
221*     .. Parameters ..
222      REAL               ZERO, ONE
223      PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
224*     ..
225*     .. Local Scalars ..
226      LOGICAL            LQUERY, SCALEA, WANTST, WANTVS
227      INTEGER            HSWORK, I, IBAL, ICOND, IERR, IEVAL, IHI, ILO,
228     $                   ITAU, IWRK, MAXWRK, MINWRK
229      REAL               ANRM, BIGNUM, CSCALE, EPS, S, SEP, SMLNUM
230*     ..
231*     .. Local Arrays ..
232      REAL               DUM( 1 )
233*     ..
234*     .. External Subroutines ..
235      EXTERNAL           CCOPY, CGEBAK, CGEBAL, CGEHRD, CHSEQR, CLACPY,
236     $                   CLASCL, CTRSEN, CUNGHR, SLABAD, XERBLA
237*     ..
238*     .. External Functions ..
239      LOGICAL            LSAME
240      INTEGER            ILAENV
241      REAL               CLANGE, SLAMCH
242      EXTERNAL           LSAME, ILAENV, CLANGE, SLAMCH
243*     ..
244*     .. Intrinsic Functions ..
245      INTRINSIC          MAX, SQRT
246*     ..
247*     .. Executable Statements ..
248*
249*     Test the input arguments
250*
251      INFO = 0
252      LQUERY = ( LWORK.EQ.-1 )
253      WANTVS = LSAME( JOBVS, 'V' )
254      WANTST = LSAME( SORT, 'S' )
255      IF( ( .NOT.WANTVS ) .AND. ( .NOT.LSAME( JOBVS, 'N' ) ) ) THEN
256         INFO = -1
257      ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN
258         INFO = -2
259      ELSE IF( N.LT.0 ) THEN
260         INFO = -4
261      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
262         INFO = -6
263      ELSE IF( LDVS.LT.1 .OR. ( WANTVS .AND. LDVS.LT.N ) ) THEN
264         INFO = -10
265      END IF
266*
267*     Compute workspace
268*      (Note: Comments in the code beginning "Workspace:" describe the
269*       minimal amount of workspace needed at that point in the code,
270*       as well as the preferred amount for good performance.
271*       CWorkspace refers to complex workspace, and RWorkspace to real
272*       workspace. NB refers to the optimal block size for the
273*       immediately following subroutine, as returned by ILAENV.
274*       HSWORK refers to the workspace preferred by CHSEQR, as
275*       calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
276*       the worst case.)
277*
278      IF( INFO.EQ.0 ) THEN
279         IF( N.EQ.0 ) THEN
280            MINWRK = 1
281            MAXWRK = 1
282         ELSE
283            MAXWRK = N + N*ILAENV( 1, 'CGEHRD', ' ', N, 1, N, 0 )
284            MINWRK = 2*N
285*
286            CALL CHSEQR( 'S', JOBVS, N, 1, N, A, LDA, W, VS, LDVS,
287     $             WORK, -1, IEVAL )
288            HSWORK = WORK( 1 )
289*
290            IF( .NOT.WANTVS ) THEN
291               MAXWRK = MAX( MAXWRK, HSWORK )
292            ELSE
293               MAXWRK = MAX( MAXWRK, N + ( N - 1 )*ILAENV( 1, 'CUNGHR',
294     $                       ' ', N, 1, N, -1 ) )
295               MAXWRK = MAX( MAXWRK, HSWORK )
296            END IF
297         END IF
298         WORK( 1 ) = MAXWRK
299*
300         IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
301            INFO = -12
302         END IF
303      END IF
304*
305      IF( INFO.NE.0 ) THEN
306         CALL XERBLA( 'CGEES ', -INFO )
307         RETURN
308      ELSE IF( LQUERY ) THEN
309         RETURN
310      END IF
311*
312*     Quick return if possible
313*
314      IF( N.EQ.0 ) THEN
315         SDIM = 0
316         RETURN
317      END IF
318*
319*     Get machine constants
320*
321      EPS = SLAMCH( 'P' )
322      SMLNUM = SLAMCH( 'S' )
323      BIGNUM = ONE / SMLNUM
324      CALL SLABAD( SMLNUM, BIGNUM )
325      SMLNUM = SQRT( SMLNUM ) / EPS
326      BIGNUM = ONE / SMLNUM
327*
328*     Scale A if max element outside range [SMLNUM,BIGNUM]
329*
330      ANRM = CLANGE( 'M', N, N, A, LDA, DUM )
331      SCALEA = .FALSE.
332      IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
333         SCALEA = .TRUE.
334         CSCALE = SMLNUM
335      ELSE IF( ANRM.GT.BIGNUM ) THEN
336         SCALEA = .TRUE.
337         CSCALE = BIGNUM
338      END IF
339      IF( SCALEA )
340     $   CALL CLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR )
341*
342*     Permute the matrix to make it more nearly triangular
343*     (CWorkspace: none)
344*     (RWorkspace: need N)
345*
346      IBAL = 1
347      CALL CGEBAL( 'P', N, A, LDA, ILO, IHI, RWORK( IBAL ), IERR )
348*
349*     Reduce to upper Hessenberg form
350*     (CWorkspace: need 2*N, prefer N+N*NB)
351*     (RWorkspace: none)
352*
353      ITAU = 1
354      IWRK = N + ITAU
355      CALL CGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ),
356     $             LWORK-IWRK+1, IERR )
357*
358      IF( WANTVS ) THEN
359*
360*        Copy Householder vectors to VS
361*
362         CALL CLACPY( 'L', N, N, A, LDA, VS, LDVS )
363*
364*        Generate unitary matrix in VS
365*        (CWorkspace: need 2*N-1, prefer N+(N-1)*NB)
366*        (RWorkspace: none)
367*
368         CALL CUNGHR( N, ILO, IHI, VS, LDVS, WORK( ITAU ), WORK( IWRK ),
369     $                LWORK-IWRK+1, IERR )
370      END IF
371*
372      SDIM = 0
373*
374*     Perform QR iteration, accumulating Schur vectors in VS if desired
375*     (CWorkspace: need 1, prefer HSWORK (see comments) )
376*     (RWorkspace: none)
377*
378      IWRK = ITAU
379      CALL CHSEQR( 'S', JOBVS, N, ILO, IHI, A, LDA, W, VS, LDVS,
380     $             WORK( IWRK ), LWORK-IWRK+1, IEVAL )
381      IF( IEVAL.GT.0 )
382     $   INFO = IEVAL
383*
384*     Sort eigenvalues if desired
385*
386      IF( WANTST .AND. INFO.EQ.0 ) THEN
387         IF( SCALEA )
388     $      CALL CLASCL( 'G', 0, 0, CSCALE, ANRM, N, 1, W, N, IERR )
389         DO 10 I = 1, N
390            BWORK( I ) = SELECT( W( I ) )
391   10    CONTINUE
392*
393*        Reorder eigenvalues and transform Schur vectors
394*        (CWorkspace: none)
395*        (RWorkspace: none)
396*
397         CALL CTRSEN( 'N', JOBVS, BWORK, N, A, LDA, VS, LDVS, W, SDIM,
398     $                S, SEP, WORK( IWRK ), LWORK-IWRK+1, ICOND )
399      END IF
400*
401      IF( WANTVS ) THEN
402*
403*        Undo balancing
404*        (CWorkspace: none)
405*        (RWorkspace: need N)
406*
407         CALL CGEBAK( 'P', 'R', N, ILO, IHI, RWORK( IBAL ), N, VS, LDVS,
408     $                IERR )
409      END IF
410*
411      IF( SCALEA ) THEN
412*
413*        Undo scaling for the Schur form of A
414*
415         CALL CLASCL( 'U', 0, 0, CSCALE, ANRM, N, N, A, LDA, IERR )
416         CALL CCOPY( N, A, LDA+1, W, 1 )
417      END IF
418*
419      WORK( 1 ) = MAXWRK
420      RETURN
421*
422*     End of CGEES
423*
424      END
425