LAPACK  3.6.1
LAPACK: Linear Algebra PACKage
dgeesx.f
Go to the documentation of this file.
1 *> \brief <b> DGEESX 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 DGEESX + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgeesx.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgeesx.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgeesx.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DGEESX( JOBVS, SORT, SELECT, SENSE, N, A, LDA, SDIM,
22 * WR, WI, VS, LDVS, RCONDE, RCONDV, WORK, LWORK,
23 * IWORK, LIWORK, BWORK, INFO )
24 *
25 * .. Scalar Arguments ..
26 * CHARACTER JOBVS, SENSE, SORT
27 * INTEGER INFO, LDA, LDVS, LIWORK, LWORK, N, SDIM
28 * DOUBLE PRECISION RCONDE, RCONDV
29 * ..
30 * .. Array Arguments ..
31 * LOGICAL BWORK( * )
32 * INTEGER IWORK( * )
33 * DOUBLE PRECISION A( LDA, * ), VS( LDVS, * ), WI( * ), WORK( * ),
34 * $ WR( * )
35 * ..
36 * .. Function Arguments ..
37 * LOGICAL SELECT
38 * EXTERNAL SELECT
39 * ..
40 *
41 *
42 *> \par Purpose:
43 * =============
44 *>
45 *> \verbatim
46 *>
47 *> DGEESX computes for an N-by-N real nonsymmetric matrix A, the
48 *> eigenvalues, the real Schur form T, and, optionally, the matrix of
49 *> Schur vectors Z. This gives the Schur factorization A = Z*T*(Z**T).
50 *>
51 *> Optionally, it also orders the eigenvalues on the diagonal of the
52 *> real Schur form so that selected eigenvalues are at the top left;
53 *> computes a reciprocal condition number for the average of the
54 *> selected eigenvalues (RCONDE); and computes a reciprocal condition
55 *> number for the right invariant subspace corresponding to the
56 *> selected eigenvalues (RCONDV). The leading columns of Z form an
57 *> orthonormal basis for this invariant subspace.
58 *>
59 *> For further explanation of the reciprocal condition numbers RCONDE
60 *> and RCONDV, see Section 4.10 of the LAPACK Users' Guide (where
61 *> these quantities are called s and sep respectively).
62 *>
63 *> A real matrix is in real Schur form if it is upper quasi-triangular
64 *> with 1-by-1 and 2-by-2 blocks. 2-by-2 blocks will be standardized in
65 *> the form
66 *> [ a b ]
67 *> [ c a ]
68 *>
69 *> where b*c < 0. The eigenvalues of such a block are a +- sqrt(bc).
70 *> \endverbatim
71 *
72 * Arguments:
73 * ==========
74 *
75 *> \param[in] JOBVS
76 *> \verbatim
77 *> JOBVS is CHARACTER*1
78 *> = 'N': Schur vectors are not computed;
79 *> = 'V': Schur vectors are computed.
80 *> \endverbatim
81 *>
82 *> \param[in] SORT
83 *> \verbatim
84 *> SORT is CHARACTER*1
85 *> Specifies whether or not to order the eigenvalues on the
86 *> diagonal of the Schur form.
87 *> = 'N': Eigenvalues are not ordered;
88 *> = 'S': Eigenvalues are ordered (see SELECT).
89 *> \endverbatim
90 *>
91 *> \param[in] SELECT
92 *> \verbatim
93 *> SELECT is a LOGICAL FUNCTION of two DOUBLE PRECISION arguments
94 *> SELECT must be declared EXTERNAL in the calling subroutine.
95 *> If SORT = 'S', SELECT is used to select eigenvalues to sort
96 *> to the top left of the Schur form.
97 *> If SORT = 'N', SELECT is not referenced.
98 *> An eigenvalue WR(j)+sqrt(-1)*WI(j) is selected if
99 *> SELECT(WR(j),WI(j)) is true; i.e., if either one of a
100 *> complex conjugate pair of eigenvalues is selected, then both
101 *> are. Note that a selected complex eigenvalue may no longer
102 *> satisfy SELECT(WR(j),WI(j)) = .TRUE. after ordering, since
103 *> ordering may change the value of complex eigenvalues
104 *> (especially if the eigenvalue is ill-conditioned); in this
105 *> case INFO may be set to N+3 (see INFO below).
106 *> \endverbatim
107 *>
108 *> \param[in] SENSE
109 *> \verbatim
110 *> SENSE is CHARACTER*1
111 *> Determines which reciprocal condition numbers are computed.
112 *> = 'N': None are computed;
113 *> = 'E': Computed for average of selected eigenvalues only;
114 *> = 'V': Computed for selected right invariant subspace only;
115 *> = 'B': Computed for both.
116 *> If SENSE = 'E', 'V' or 'B', SORT must equal 'S'.
117 *> \endverbatim
118 *>
119 *> \param[in] N
120 *> \verbatim
121 *> N is INTEGER
122 *> The order of the matrix A. N >= 0.
123 *> \endverbatim
124 *>
125 *> \param[in,out] A
126 *> \verbatim
127 *> A is DOUBLE PRECISION array, dimension (LDA, N)
128 *> On entry, the N-by-N matrix A.
129 *> On exit, A is overwritten by its real Schur form T.
130 *> \endverbatim
131 *>
132 *> \param[in] LDA
133 *> \verbatim
134 *> LDA is INTEGER
135 *> The leading dimension of the array A. LDA >= max(1,N).
136 *> \endverbatim
137 *>
138 *> \param[out] SDIM
139 *> \verbatim
140 *> SDIM is INTEGER
141 *> If SORT = 'N', SDIM = 0.
142 *> If SORT = 'S', SDIM = number of eigenvalues (after sorting)
143 *> for which SELECT is true. (Complex conjugate
144 *> pairs for which SELECT is true for either
145 *> eigenvalue count as 2.)
146 *> \endverbatim
147 *>
148 *> \param[out] WR
149 *> \verbatim
150 *> WR is DOUBLE PRECISION array, dimension (N)
151 *> \endverbatim
152 *>
153 *> \param[out] WI
154 *> \verbatim
155 *> WI is DOUBLE PRECISION array, dimension (N)
156 *> WR and WI contain the real and imaginary parts, respectively,
157 *> of the computed eigenvalues, in the same order that they
158 *> appear on the diagonal of the output Schur form T. Complex
159 *> conjugate pairs of eigenvalues appear consecutively with the
160 *> eigenvalue having the positive imaginary part first.
161 *> \endverbatim
162 *>
163 *> \param[out] VS
164 *> \verbatim
165 *> VS is DOUBLE PRECISION array, dimension (LDVS,N)
166 *> If JOBVS = 'V', VS contains the orthogonal matrix Z of Schur
167 *> vectors.
168 *> If JOBVS = 'N', VS is not referenced.
169 *> \endverbatim
170 *>
171 *> \param[in] LDVS
172 *> \verbatim
173 *> LDVS is INTEGER
174 *> The leading dimension of the array VS. LDVS >= 1, and if
175 *> JOBVS = 'V', LDVS >= N.
176 *> \endverbatim
177 *>
178 *> \param[out] RCONDE
179 *> \verbatim
180 *> RCONDE is DOUBLE PRECISION
181 *> If SENSE = 'E' or 'B', RCONDE contains the reciprocal
182 *> condition number for the average of the selected eigenvalues.
183 *> Not referenced if SENSE = 'N' or 'V'.
184 *> \endverbatim
185 *>
186 *> \param[out] RCONDV
187 *> \verbatim
188 *> RCONDV is DOUBLE PRECISION
189 *> If SENSE = 'V' or 'B', RCONDV contains the reciprocal
190 *> condition number for the selected right invariant subspace.
191 *> Not referenced if SENSE = 'N' or 'E'.
192 *> \endverbatim
193 *>
194 *> \param[out] WORK
195 *> \verbatim
196 *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
197 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
198 *> \endverbatim
199 *>
200 *> \param[in] LWORK
201 *> \verbatim
202 *> LWORK is INTEGER
203 *> The dimension of the array WORK. LWORK >= max(1,3*N).
204 *> Also, if SENSE = 'E' or 'V' or 'B',
205 *> LWORK >= N+2*SDIM*(N-SDIM), where SDIM is the number of
206 *> selected eigenvalues computed by this routine. Note that
207 *> N+2*SDIM*(N-SDIM) <= N+N*N/2. Note also that an error is only
208 *> returned if LWORK < max(1,3*N), but if SENSE = 'E' or 'V' or
209 *> 'B' this may not be large enough.
210 *> For good performance, LWORK must generally be larger.
211 *>
212 *> If LWORK = -1, then a workspace query is assumed; the routine
213 *> only calculates upper bounds on the optimal sizes of the
214 *> arrays WORK and IWORK, returns these values as the first
215 *> entries of the WORK and IWORK arrays, and no error messages
216 *> related to LWORK or LIWORK are issued by XERBLA.
217 *> \endverbatim
218 *>
219 *> \param[out] IWORK
220 *> \verbatim
221 *> IWORK is INTEGER array, dimension (MAX(1,LIWORK))
222 *> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
223 *> \endverbatim
224 *>
225 *> \param[in] LIWORK
226 *> \verbatim
227 *> LIWORK is INTEGER
228 *> The dimension of the array IWORK.
229 *> LIWORK >= 1; if SENSE = 'V' or 'B', LIWORK >= SDIM*(N-SDIM).
230 *> Note that SDIM*(N-SDIM) <= N*N/4. Note also that an error is
231 *> only returned if LIWORK < 1, but if SENSE = 'V' or 'B' this
232 *> may not be large enough.
233 *>
234 *> If LIWORK = -1, then a workspace query is assumed; the
235 *> routine only calculates upper bounds on the optimal sizes of
236 *> the arrays WORK and IWORK, returns these values as the first
237 *> entries of the WORK and IWORK arrays, and no error messages
238 *> related to LWORK or LIWORK are issued by XERBLA.
239 *> \endverbatim
240 *>
241 *> \param[out] BWORK
242 *> \verbatim
243 *> BWORK is LOGICAL array, dimension (N)
244 *> Not referenced if SORT = 'N'.
245 *> \endverbatim
246 *>
247 *> \param[out] INFO
248 *> \verbatim
249 *> INFO is INTEGER
250 *> = 0: successful exit
251 *> < 0: if INFO = -i, the i-th argument had an illegal value.
252 *> > 0: if INFO = i, and i is
253 *> <= N: the QR algorithm failed to compute all the
254 *> eigenvalues; elements 1:ILO-1 and i+1:N of WR and WI
255 *> contain those eigenvalues which have converged; if
256 *> JOBVS = 'V', VS contains the transformation which
257 *> reduces A to its partially converged Schur form.
258 *> = N+1: the eigenvalues could not be reordered because some
259 *> eigenvalues were too close to separate (the problem
260 *> is very ill-conditioned);
261 *> = N+2: after reordering, roundoff changed values of some
262 *> complex eigenvalues so that leading eigenvalues in
263 *> the Schur form no longer satisfy SELECT=.TRUE. This
264 *> could also be caused by underflow due to scaling.
265 *> \endverbatim
266 *
267 * Authors:
268 * ========
269 *
270 *> \author Univ. of Tennessee
271 *> \author Univ. of California Berkeley
272 *> \author Univ. of Colorado Denver
273 *> \author NAG Ltd.
274 *
275 *> \date June 2016
276 *
277 *> \ingroup doubleGEeigen
278 *
279 * =====================================================================
280  SUBROUTINE dgeesx( JOBVS, SORT, SELECT, SENSE, N, A, LDA, SDIM,
281  $ wr, wi, vs, ldvs, rconde, rcondv, work, lwork,
282  $ iwork, liwork, bwork, info )
283 *
284 * -- LAPACK driver routine (version 3.6.1) --
285 * -- LAPACK is a software package provided by Univ. of Tennessee, --
286 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
287 * June 2016
288 *
289 * .. Scalar Arguments ..
290  CHARACTER JOBVS, SENSE, SORT
291  INTEGER INFO, LDA, LDVS, LIWORK, LWORK, N, SDIM
292  DOUBLE PRECISION RCONDE, RCONDV
293 * ..
294 * .. Array Arguments ..
295  LOGICAL BWORK( * )
296  INTEGER IWORK( * )
297  DOUBLE PRECISION A( lda, * ), VS( ldvs, * ), WI( * ), WORK( * ),
298  $ wr( * )
299 * ..
300 * .. Function Arguments ..
301  LOGICAL SELECT
302  EXTERNAL SELECT
303 * ..
304 *
305 * =====================================================================
306 *
307 * .. Parameters ..
308  DOUBLE PRECISION ZERO, ONE
309  parameter ( zero = 0.0d0, one = 1.0d0 )
310 * ..
311 * .. Local Scalars ..
312  LOGICAL CURSL, LASTSL, LQUERY, LST2SL, SCALEA, WANTSB,
313  $ wantse, wantsn, wantst, wantsv, wantvs
314  INTEGER HSWORK, I, I1, I2, IBAL, ICOND, IERR, IEVAL,
315  $ ihi, ilo, inxt, ip, itau, iwrk, liwrk, lwrk,
316  $ maxwrk, minwrk
317  DOUBLE PRECISION ANRM, BIGNUM, CSCALE, EPS, SMLNUM
318 * ..
319 * .. Local Arrays ..
320  DOUBLE PRECISION DUM( 1 )
321 * ..
322 * .. External Subroutines ..
323  EXTERNAL dcopy, dgebak, dgebal, dgehrd, dhseqr, dlacpy,
325 * ..
326 * .. External Functions ..
327  LOGICAL LSAME
328  INTEGER ILAENV
329  DOUBLE PRECISION DLAMCH, DLANGE
330  EXTERNAL lsame, ilaenv, dlabad, dlamch, dlange
331 * ..
332 * .. Intrinsic Functions ..
333  INTRINSIC max, sqrt
334 * ..
335 * .. Executable Statements ..
336 *
337 * Test the input arguments
338 *
339  info = 0
340  wantvs = lsame( jobvs, 'V' )
341  wantst = lsame( sort, 'S' )
342  wantsn = lsame( sense, 'N' )
343  wantse = lsame( sense, 'E' )
344  wantsv = lsame( sense, 'V' )
345  wantsb = lsame( sense, 'B' )
346  lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
347 *
348  IF( ( .NOT.wantvs ) .AND. ( .NOT.lsame( jobvs, 'N' ) ) ) THEN
349  info = -1
350  ELSE IF( ( .NOT.wantst ) .AND. ( .NOT.lsame( sort, 'N' ) ) ) THEN
351  info = -2
352  ELSE IF( .NOT.( wantsn .OR. wantse .OR. wantsv .OR. wantsb ) .OR.
353  $ ( .NOT.wantst .AND. .NOT.wantsn ) ) THEN
354  info = -4
355  ELSE IF( n.LT.0 ) THEN
356  info = -5
357  ELSE IF( lda.LT.max( 1, n ) ) THEN
358  info = -7
359  ELSE IF( ldvs.LT.1 .OR. ( wantvs .AND. ldvs.LT.n ) ) THEN
360  info = -12
361  END IF
362 *
363 * Compute workspace
364 * (Note: Comments in the code beginning "RWorkspace:" describe the
365 * minimal amount of real workspace needed at that point in the
366 * code, as well as the preferred amount for good performance.
367 * IWorkspace refers to integer workspace.
368 * NB refers to the optimal block size for the immediately
369 * following subroutine, as returned by ILAENV.
370 * HSWORK refers to the workspace preferred by DHSEQR, as
371 * calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
372 * the worst case.
373 * If SENSE = 'E', 'V' or 'B', then the amount of workspace needed
374 * depends on SDIM, which is computed by the routine DTRSEN later
375 * in the code.)
376 *
377  IF( info.EQ.0 ) THEN
378  liwrk = 1
379  IF( n.EQ.0 ) THEN
380  minwrk = 1
381  lwrk = 1
382  ELSE
383  maxwrk = 2*n + n*ilaenv( 1, 'DGEHRD', ' ', n, 1, n, 0 )
384  minwrk = 3*n
385 *
386  CALL dhseqr( 'S', jobvs, n, 1, n, a, lda, wr, wi, vs, ldvs,
387  $ work, -1, ieval )
388  hswork = work( 1 )
389 *
390  IF( .NOT.wantvs ) THEN
391  maxwrk = max( maxwrk, n + hswork )
392  ELSE
393  maxwrk = max( maxwrk, 2*n + ( n - 1 )*ilaenv( 1,
394  $ 'DORGHR', ' ', n, 1, n, -1 ) )
395  maxwrk = max( maxwrk, n + hswork )
396  END IF
397  lwrk = maxwrk
398  IF( .NOT.wantsn )
399  $ lwrk = max( lwrk, n + ( n*n )/2 )
400  IF( wantsv .OR. wantsb )
401  $ liwrk = ( n*n )/4
402  END IF
403  iwork( 1 ) = liwrk
404  work( 1 ) = lwrk
405 *
406  IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
407  info = -16
408  ELSE IF( liwork.LT.1 .AND. .NOT.lquery ) THEN
409  info = -18
410  END IF
411  END IF
412 *
413  IF( info.NE.0 ) THEN
414  CALL xerbla( 'DGEESX', -info )
415  RETURN
416  ELSE IF( lquery ) THEN
417  RETURN
418  END IF
419 *
420 * Quick return if possible
421 *
422  IF( n.EQ.0 ) THEN
423  sdim = 0
424  RETURN
425  END IF
426 *
427 * Get machine constants
428 *
429  eps = dlamch( 'P' )
430  smlnum = dlamch( 'S' )
431  bignum = one / smlnum
432  CALL dlabad( smlnum, bignum )
433  smlnum = sqrt( smlnum ) / eps
434  bignum = one / smlnum
435 *
436 * Scale A if max element outside range [SMLNUM,BIGNUM]
437 *
438  anrm = dlange( 'M', n, n, a, lda, dum )
439  scalea = .false.
440  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
441  scalea = .true.
442  cscale = smlnum
443  ELSE IF( anrm.GT.bignum ) THEN
444  scalea = .true.
445  cscale = bignum
446  END IF
447  IF( scalea )
448  $ CALL dlascl( 'G', 0, 0, anrm, cscale, n, n, a, lda, ierr )
449 *
450 * Permute the matrix to make it more nearly triangular
451 * (RWorkspace: need N)
452 *
453  ibal = 1
454  CALL dgebal( 'P', n, a, lda, ilo, ihi, work( ibal ), ierr )
455 *
456 * Reduce to upper Hessenberg form
457 * (RWorkspace: need 3*N, prefer 2*N+N*NB)
458 *
459  itau = n + ibal
460  iwrk = n + itau
461  CALL dgehrd( n, ilo, ihi, a, lda, work( itau ), work( iwrk ),
462  $ lwork-iwrk+1, ierr )
463 *
464  IF( wantvs ) THEN
465 *
466 * Copy Householder vectors to VS
467 *
468  CALL dlacpy( 'L', n, n, a, lda, vs, ldvs )
469 *
470 * Generate orthogonal matrix in VS
471 * (RWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
472 *
473  CALL dorghr( n, ilo, ihi, vs, ldvs, work( itau ), work( iwrk ),
474  $ lwork-iwrk+1, ierr )
475  END IF
476 *
477  sdim = 0
478 *
479 * Perform QR iteration, accumulating Schur vectors in VS if desired
480 * (RWorkspace: need N+1, prefer N+HSWORK (see comments) )
481 *
482  iwrk = itau
483  CALL dhseqr( 'S', jobvs, n, ilo, ihi, a, lda, wr, wi, vs, ldvs,
484  $ work( iwrk ), lwork-iwrk+1, ieval )
485  IF( ieval.GT.0 )
486  $ info = ieval
487 *
488 * Sort eigenvalues if desired
489 *
490  IF( wantst .AND. info.EQ.0 ) THEN
491  IF( scalea ) THEN
492  CALL dlascl( 'G', 0, 0, cscale, anrm, n, 1, wr, n, ierr )
493  CALL dlascl( 'G', 0, 0, cscale, anrm, n, 1, wi, n, ierr )
494  END IF
495  DO 10 i = 1, n
496  bwork( i ) = SELECT( wr( i ), wi( i ) )
497  10 CONTINUE
498 *
499 * Reorder eigenvalues, transform Schur vectors, and compute
500 * reciprocal condition numbers
501 * (RWorkspace: if SENSE is not 'N', need N+2*SDIM*(N-SDIM)
502 * otherwise, need N )
503 * (IWorkspace: if SENSE is 'V' or 'B', need SDIM*(N-SDIM)
504 * otherwise, need 0 )
505 *
506  CALL dtrsen( sense, jobvs, bwork, n, a, lda, vs, ldvs, wr, wi,
507  $ sdim, rconde, rcondv, work( iwrk ), lwork-iwrk+1,
508  $ iwork, liwork, icond )
509  IF( .NOT.wantsn )
510  $ maxwrk = max( maxwrk, n+2*sdim*( n-sdim ) )
511  IF( icond.EQ.-15 ) THEN
512 *
513 * Not enough real workspace
514 *
515  info = -16
516  ELSE IF( icond.EQ.-17 ) THEN
517 *
518 * Not enough integer workspace
519 *
520  info = -18
521  ELSE IF( icond.GT.0 ) THEN
522 *
523 * DTRSEN failed to reorder or to restore standard Schur form
524 *
525  info = icond + n
526  END IF
527  END IF
528 *
529  IF( wantvs ) THEN
530 *
531 * Undo balancing
532 * (RWorkspace: need N)
533 *
534  CALL dgebak( 'P', 'R', n, ilo, ihi, work( ibal ), n, vs, ldvs,
535  $ ierr )
536  END IF
537 *
538  IF( scalea ) THEN
539 *
540 * Undo scaling for the Schur form of A
541 *
542  CALL dlascl( 'H', 0, 0, cscale, anrm, n, n, a, lda, ierr )
543  CALL dcopy( n, a, lda+1, wr, 1 )
544  IF( ( wantsv .OR. wantsb ) .AND. info.EQ.0 ) THEN
545  dum( 1 ) = rcondv
546  CALL dlascl( 'G', 0, 0, cscale, anrm, 1, 1, dum, 1, ierr )
547  rcondv = dum( 1 )
548  END IF
549  IF( cscale.EQ.smlnum ) THEN
550 *
551 * If scaling back towards underflow, adjust WI if an
552 * offdiagonal element of a 2-by-2 block in the Schur form
553 * underflows.
554 *
555  IF( ieval.GT.0 ) THEN
556  i1 = ieval + 1
557  i2 = ihi - 1
558  CALL dlascl( 'G', 0, 0, cscale, anrm, ilo-1, 1, wi, n,
559  $ ierr )
560  ELSE IF( wantst ) THEN
561  i1 = 1
562  i2 = n - 1
563  ELSE
564  i1 = ilo
565  i2 = ihi - 1
566  END IF
567  inxt = i1 - 1
568  DO 20 i = i1, i2
569  IF( i.LT.inxt )
570  $ GO TO 20
571  IF( wi( i ).EQ.zero ) THEN
572  inxt = i + 1
573  ELSE
574  IF( a( i+1, i ).EQ.zero ) THEN
575  wi( i ) = zero
576  wi( i+1 ) = zero
577  ELSE IF( a( i+1, i ).NE.zero .AND. a( i, i+1 ).EQ.
578  $ zero ) THEN
579  wi( i ) = zero
580  wi( i+1 ) = zero
581  IF( i.GT.1 )
582  $ CALL dswap( i-1, a( 1, i ), 1, a( 1, i+1 ), 1 )
583  IF( n.GT.i+1 )
584  $ CALL dswap( n-i-1, a( i, i+2 ), lda,
585  $ a( i+1, i+2 ), lda )
586  CALL dswap( n, vs( 1, i ), 1, vs( 1, i+1 ), 1 )
587  a( i, i+1 ) = a( i+1, i )
588  a( i+1, i ) = zero
589  END IF
590  inxt = i + 2
591  END IF
592  20 CONTINUE
593  END IF
594  CALL dlascl( 'G', 0, 0, cscale, anrm, n-ieval, 1,
595  $ wi( ieval+1 ), max( n-ieval, 1 ), ierr )
596  END IF
597 *
598  IF( wantst .AND. info.EQ.0 ) THEN
599 *
600 * Check if reordering successful
601 *
602  lastsl = .true.
603  lst2sl = .true.
604  sdim = 0
605  ip = 0
606  DO 30 i = 1, n
607  cursl = SELECT( wr( i ), wi( i ) )
608  IF( wi( i ).EQ.zero ) THEN
609  IF( cursl )
610  $ sdim = sdim + 1
611  ip = 0
612  IF( cursl .AND. .NOT.lastsl )
613  $ info = n + 2
614  ELSE
615  IF( ip.EQ.1 ) THEN
616 *
617 * Last eigenvalue of conjugate pair
618 *
619  cursl = cursl .OR. lastsl
620  lastsl = cursl
621  IF( cursl )
622  $ sdim = sdim + 2
623  ip = -1
624  IF( cursl .AND. .NOT.lst2sl )
625  $ info = n + 2
626  ELSE
627 *
628 * First eigenvalue of conjugate pair
629 *
630  ip = 1
631  END IF
632  END IF
633  lst2sl = lastsl
634  lastsl = cursl
635  30 CONTINUE
636  END IF
637 *
638  work( 1 ) = maxwrk
639  IF( wantsv .OR. wantsb ) THEN
640  iwork( 1 ) = max( 1, sdim*( n-sdim ) )
641  ELSE
642  iwork( 1 ) = 1
643  END IF
644 *
645  RETURN
646 *
647 * End of DGEESX
648 *
649  END
subroutine dgeesx(JOBVS, SORT, SELECT, SENSE, N, A, LDA, SDIM, WR, WI, VS, LDVS, RCONDE, RCONDV, WORK, LWORK, IWORK, LIWORK, BWORK, INFO)
DGEESX computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE...
Definition: dgeesx.f:283
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:53
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 dgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
DGEHRD
Definition: dgehrd.f:169
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 dgebak(JOB, SIDE, N, ILO, IHI, SCALE, M, V, LDV, INFO)
DGEBAK
Definition: dgebak.f:132
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:62
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:76
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
Definition: dswap.f:53
subroutine dgebal(JOB, N, A, LDA, ILO, IHI, SCALE, INFO)
DGEBAL
Definition: dgebal.f:162
subroutine dorghr(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
DORGHR
Definition: dorghr.f:128
subroutine dtrsen(JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, WR, WI, M, S, SEP, WORK, LWORK, IWORK, LIWORK, INFO)
DTRSEN
Definition: dtrsen.f:315
subroutine dhseqr(JOB, COMPZ, N, ILO, IHI, H, LDH, WR, WI, Z, LDZ, WORK, LWORK, INFO)
DHSEQR
Definition: dhseqr.f:318