 LAPACK  3.9.1 LAPACK: Linear Algebra PACKage

## ◆ dgeesx()

 subroutine dgeesx ( character JOBVS, character SORT, external SELECT, character SENSE, integer N, double precision, dimension( lda, * ) A, integer LDA, integer SDIM, double precision, dimension( * ) WR, double precision, dimension( * ) WI, double precision, dimension( ldvs, * ) VS, integer LDVS, double precision RCONDE, double precision RCONDV, double precision, dimension( * ) WORK, integer LWORK, integer, dimension( * ) IWORK, integer LIWORK, logical, dimension( * ) BWORK, integer INFO )

DGEESX computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE matrices

Download DGEESX + dependencies [TGZ] [ZIP] [TXT]

Purpose:
``` DGEESX computes for an N-by-N real nonsymmetric matrix A, the
eigenvalues, the real Schur form T, and, optionally, the matrix of
Schur vectors Z.  This gives the Schur factorization A = Z*T*(Z**T).

Optionally, it also orders the eigenvalues on the diagonal of the
real Schur form so that selected eigenvalues are at the top left;
computes a reciprocal condition number for the average of the
selected eigenvalues (RCONDE); and computes a reciprocal condition
number for the right invariant subspace corresponding to the
selected eigenvalues (RCONDV).  The leading columns of Z form an
orthonormal basis for this invariant subspace.

For further explanation of the reciprocal condition numbers RCONDE
and RCONDV, see Section 4.10 of the LAPACK Users' Guide (where
these quantities are called s and sep respectively).

A real matrix is in real Schur form if it is upper quasi-triangular
with 1-by-1 and 2-by-2 blocks. 2-by-2 blocks will be standardized in
the form
[  a  b  ]
[  c  a  ]

where b*c < 0. The eigenvalues of such a block are a +- sqrt(bc).```
Parameters
 [in] JOBVS ``` JOBVS is CHARACTER*1 = 'N': Schur vectors are not computed; = 'V': Schur vectors are computed.``` [in] SORT ``` SORT is CHARACTER*1 Specifies whether or not to order the eigenvalues on the diagonal of the Schur form. = 'N': Eigenvalues are not ordered; = 'S': Eigenvalues are ordered (see SELECT).``` [in] SELECT ``` SELECT is a LOGICAL FUNCTION of two DOUBLE PRECISION arguments SELECT must be declared EXTERNAL in the calling subroutine. If SORT = 'S', SELECT is used to select eigenvalues to sort to the top left of the Schur form. If SORT = 'N', SELECT is not referenced. An eigenvalue WR(j)+sqrt(-1)*WI(j) is selected if SELECT(WR(j),WI(j)) is true; i.e., if either one of a complex conjugate pair of eigenvalues is selected, then both are. Note that a selected complex eigenvalue may no longer satisfy SELECT(WR(j),WI(j)) = .TRUE. after ordering, since ordering may change the value of complex eigenvalues (especially if the eigenvalue is ill-conditioned); in this case INFO may be set to N+3 (see INFO below).``` [in] SENSE ``` SENSE is CHARACTER*1 Determines which reciprocal condition numbers are computed. = 'N': None are computed; = 'E': Computed for average of selected eigenvalues only; = 'V': Computed for selected right invariant subspace only; = 'B': Computed for both. If SENSE = 'E', 'V' or 'B', SORT must equal 'S'.``` [in] N ``` N is INTEGER The order of the matrix A. N >= 0.``` [in,out] A ``` A is DOUBLE PRECISION array, dimension (LDA, N) On entry, the N-by-N matrix A. On exit, A is overwritten by its real Schur form T.``` [in] LDA ``` LDA is INTEGER The leading dimension of the array A. LDA >= max(1,N).``` [out] SDIM ``` SDIM is INTEGER If SORT = 'N', SDIM = 0. If SORT = 'S', SDIM = number of eigenvalues (after sorting) for which SELECT is true. (Complex conjugate pairs for which SELECT is true for either eigenvalue count as 2.)``` [out] WR ` WR is DOUBLE PRECISION array, dimension (N)` [out] WI ``` WI is DOUBLE PRECISION array, dimension (N) WR and WI contain the real and imaginary parts, respectively, of the computed eigenvalues, in the same order that they appear on the diagonal of the output Schur form T. Complex conjugate pairs of eigenvalues appear consecutively with the eigenvalue having the positive imaginary part first.``` [out] VS ``` VS is DOUBLE PRECISION array, dimension (LDVS,N) If JOBVS = 'V', VS contains the orthogonal matrix Z of Schur vectors. If JOBVS = 'N', VS is not referenced.``` [in] LDVS ``` LDVS is INTEGER The leading dimension of the array VS. LDVS >= 1, and if JOBVS = 'V', LDVS >= N.``` [out] RCONDE ``` RCONDE is DOUBLE PRECISION If SENSE = 'E' or 'B', RCONDE contains the reciprocal condition number for the average of the selected eigenvalues. Not referenced if SENSE = 'N' or 'V'.``` [out] RCONDV ``` RCONDV is DOUBLE PRECISION If SENSE = 'V' or 'B', RCONDV contains the reciprocal condition number for the selected right invariant subspace. Not referenced if SENSE = 'N' or 'E'.``` [out] WORK ``` WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK)) On exit, if INFO = 0, WORK(1) returns the optimal LWORK.``` [in] LWORK ``` LWORK is INTEGER The dimension of the array WORK. LWORK >= max(1,3*N). Also, if SENSE = 'E' or 'V' or 'B', LWORK >= N+2*SDIM*(N-SDIM), where SDIM is the number of selected eigenvalues computed by this routine. Note that N+2*SDIM*(N-SDIM) <= N+N*N/2. Note also that an error is only returned if LWORK < max(1,3*N), but if SENSE = 'E' or 'V' or 'B' this may not be large enough. For good performance, LWORK must generally be larger. If LWORK = -1, then a workspace query is assumed; the routine only calculates upper bounds on the optimal sizes of the arrays WORK and IWORK, returns these values as the first entries of the WORK and IWORK arrays, and no error messages related to LWORK or LIWORK are issued by XERBLA.``` [out] IWORK ``` IWORK is INTEGER array, dimension (MAX(1,LIWORK)) On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.``` [in] LIWORK ``` LIWORK is INTEGER The dimension of the array IWORK. LIWORK >= 1; if SENSE = 'V' or 'B', LIWORK >= SDIM*(N-SDIM). Note that SDIM*(N-SDIM) <= N*N/4. Note also that an error is only returned if LIWORK < 1, but if SENSE = 'V' or 'B' this may not be large enough. If LIWORK = -1, then a workspace query is assumed; the routine only calculates upper bounds on the optimal sizes of the arrays WORK and IWORK, returns these values as the first entries of the WORK and IWORK arrays, and no error messages related to LWORK or LIWORK are issued by XERBLA.``` [out] BWORK ``` BWORK is LOGICAL array, dimension (N) Not referenced if SORT = 'N'.``` [out] INFO ``` INFO is INTEGER = 0: successful exit < 0: if INFO = -i, the i-th argument had an illegal value. > 0: if INFO = i, and i is <= N: the QR algorithm failed to compute all the eigenvalues; elements 1:ILO-1 and i+1:N of WR and WI contain those eigenvalues which have converged; if JOBVS = 'V', VS contains the transformation which reduces A to its partially converged Schur form. = N+1: the eigenvalues could not be reordered because some eigenvalues were too close to separate (the problem is very ill-conditioned); = N+2: after reordering, roundoff changed values of some complex eigenvalues so that leading eigenvalues in the Schur form no longer satisfy SELECT=.TRUE. This could also be caused by underflow due to scaling.```

Definition at line 278 of file dgeesx.f.

281 *
282 * -- LAPACK driver routine --
283 * -- LAPACK is a software package provided by Univ. of Tennessee, --
284 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
285 *
286 * .. Scalar Arguments ..
287  CHARACTER JOBVS, SENSE, SORT
288  INTEGER INFO, LDA, LDVS, LIWORK, LWORK, N, SDIM
289  DOUBLE PRECISION RCONDE, RCONDV
290 * ..
291 * .. Array Arguments ..
292  LOGICAL BWORK( * )
293  INTEGER IWORK( * )
294  DOUBLE PRECISION A( LDA, * ), VS( LDVS, * ), WI( * ), WORK( * ),
295  \$ WR( * )
296 * ..
297 * .. Function Arguments ..
298  LOGICAL SELECT
299  EXTERNAL SELECT
300 * ..
301 *
302 * =====================================================================
303 *
304 * .. Parameters ..
305  DOUBLE PRECISION ZERO, ONE
306  parameter( zero = 0.0d0, one = 1.0d0 )
307 * ..
308 * .. Local Scalars ..
309  LOGICAL CURSL, LASTSL, LQUERY, LST2SL, SCALEA, WANTSB,
310  \$ WANTSE, WANTSN, WANTST, WANTSV, WANTVS
311  INTEGER HSWORK, I, I1, I2, IBAL, ICOND, IERR, IEVAL,
312  \$ IHI, ILO, INXT, IP, ITAU, IWRK, LIWRK, LWRK,
313  \$ MAXWRK, MINWRK
314  DOUBLE PRECISION ANRM, BIGNUM, CSCALE, EPS, SMLNUM
315 * ..
316 * .. Local Arrays ..
317  DOUBLE PRECISION DUM( 1 )
318 * ..
319 * .. External Subroutines ..
320  EXTERNAL dcopy, dgebak, dgebal, dgehrd, dhseqr, dlacpy,
322 * ..
323 * .. External Functions ..
324  LOGICAL LSAME
325  INTEGER ILAENV
326  DOUBLE PRECISION DLAMCH, DLANGE
327  EXTERNAL lsame, ilaenv, dlabad, dlamch, dlange
328 * ..
329 * .. Intrinsic Functions ..
330  INTRINSIC max, sqrt
331 * ..
332 * .. Executable Statements ..
333 *
334 * Test the input arguments
335 *
336  info = 0
337  wantvs = lsame( jobvs, 'V' )
338  wantst = lsame( sort, 'S' )
339  wantsn = lsame( sense, 'N' )
340  wantse = lsame( sense, 'E' )
341  wantsv = lsame( sense, 'V' )
342  wantsb = lsame( sense, 'B' )
343  lquery = ( lwork.EQ.-1 .OR. liwork.EQ.-1 )
344 *
345  IF( ( .NOT.wantvs ) .AND. ( .NOT.lsame( jobvs, 'N' ) ) ) THEN
346  info = -1
347  ELSE IF( ( .NOT.wantst ) .AND. ( .NOT.lsame( sort, 'N' ) ) ) THEN
348  info = -2
349  ELSE IF( .NOT.( wantsn .OR. wantse .OR. wantsv .OR. wantsb ) .OR.
350  \$ ( .NOT.wantst .AND. .NOT.wantsn ) ) THEN
351  info = -4
352  ELSE IF( n.LT.0 ) THEN
353  info = -5
354  ELSE IF( lda.LT.max( 1, n ) ) THEN
355  info = -7
356  ELSE IF( ldvs.LT.1 .OR. ( wantvs .AND. ldvs.LT.n ) ) THEN
357  info = -12
358  END IF
359 *
360 * Compute workspace
361 * (Note: Comments in the code beginning "RWorkspace:" describe the
362 * minimal amount of real workspace needed at that point in the
363 * code, as well as the preferred amount for good performance.
364 * IWorkspace refers to integer workspace.
365 * NB refers to the optimal block size for the immediately
366 * following subroutine, as returned by ILAENV.
367 * HSWORK refers to the workspace preferred by DHSEQR, as
368 * calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
369 * the worst case.
370 * If SENSE = 'E', 'V' or 'B', then the amount of workspace needed
371 * depends on SDIM, which is computed by the routine DTRSEN later
372 * in the code.)
373 *
374  IF( info.EQ.0 ) THEN
375  liwrk = 1
376  IF( n.EQ.0 ) THEN
377  minwrk = 1
378  lwrk = 1
379  ELSE
380  maxwrk = 2*n + n*ilaenv( 1, 'DGEHRD', ' ', n, 1, n, 0 )
381  minwrk = 3*n
382 *
383  CALL dhseqr( 'S', jobvs, n, 1, n, a, lda, wr, wi, vs, ldvs,
384  \$ work, -1, ieval )
385  hswork = work( 1 )
386 *
387  IF( .NOT.wantvs ) THEN
388  maxwrk = max( maxwrk, n + hswork )
389  ELSE
390  maxwrk = max( maxwrk, 2*n + ( n - 1 )*ilaenv( 1,
391  \$ 'DORGHR', ' ', n, 1, n, -1 ) )
392  maxwrk = max( maxwrk, n + hswork )
393  END IF
394  lwrk = maxwrk
395  IF( .NOT.wantsn )
396  \$ lwrk = max( lwrk, n + ( n*n )/2 )
397  IF( wantsv .OR. wantsb )
398  \$ liwrk = ( n*n )/4
399  END IF
400  iwork( 1 ) = liwrk
401  work( 1 ) = lwrk
402 *
403  IF( lwork.LT.minwrk .AND. .NOT.lquery ) THEN
404  info = -16
405  ELSE IF( liwork.LT.1 .AND. .NOT.lquery ) THEN
406  info = -18
407  END IF
408  END IF
409 *
410  IF( info.NE.0 ) THEN
411  CALL xerbla( 'DGEESX', -info )
412  RETURN
413  ELSE IF( lquery ) THEN
414  RETURN
415  END IF
416 *
417 * Quick return if possible
418 *
419  IF( n.EQ.0 ) THEN
420  sdim = 0
421  RETURN
422  END IF
423 *
424 * Get machine constants
425 *
426  eps = dlamch( 'P' )
427  smlnum = dlamch( 'S' )
428  bignum = one / smlnum
429  CALL dlabad( smlnum, bignum )
430  smlnum = sqrt( smlnum ) / eps
431  bignum = one / smlnum
432 *
433 * Scale A if max element outside range [SMLNUM,BIGNUM]
434 *
435  anrm = dlange( 'M', n, n, a, lda, dum )
436  scalea = .false.
437  IF( anrm.GT.zero .AND. anrm.LT.smlnum ) THEN
438  scalea = .true.
439  cscale = smlnum
440  ELSE IF( anrm.GT.bignum ) THEN
441  scalea = .true.
442  cscale = bignum
443  END IF
444  IF( scalea )
445  \$ CALL dlascl( 'G', 0, 0, anrm, cscale, n, n, a, lda, ierr )
446 *
447 * Permute the matrix to make it more nearly triangular
448 * (RWorkspace: need N)
449 *
450  ibal = 1
451  CALL dgebal( 'P', n, a, lda, ilo, ihi, work( ibal ), ierr )
452 *
453 * Reduce to upper Hessenberg form
454 * (RWorkspace: need 3*N, prefer 2*N+N*NB)
455 *
456  itau = n + ibal
457  iwrk = n + itau
458  CALL dgehrd( n, ilo, ihi, a, lda, work( itau ), work( iwrk ),
459  \$ lwork-iwrk+1, ierr )
460 *
461  IF( wantvs ) THEN
462 *
463 * Copy Householder vectors to VS
464 *
465  CALL dlacpy( 'L', n, n, a, lda, vs, ldvs )
466 *
467 * Generate orthogonal matrix in VS
468 * (RWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
469 *
470  CALL dorghr( n, ilo, ihi, vs, ldvs, work( itau ), work( iwrk ),
471  \$ lwork-iwrk+1, ierr )
472  END IF
473 *
474  sdim = 0
475 *
476 * Perform QR iteration, accumulating Schur vectors in VS if desired
477 * (RWorkspace: need N+1, prefer N+HSWORK (see comments) )
478 *
479  iwrk = itau
480  CALL dhseqr( 'S', jobvs, n, ilo, ihi, a, lda, wr, wi, vs, ldvs,
481  \$ work( iwrk ), lwork-iwrk+1, ieval )
482  IF( ieval.GT.0 )
483  \$ info = ieval
484 *
485 * Sort eigenvalues if desired
486 *
487  IF( wantst .AND. info.EQ.0 ) THEN
488  IF( scalea ) THEN
489  CALL dlascl( 'G', 0, 0, cscale, anrm, n, 1, wr, n, ierr )
490  CALL dlascl( 'G', 0, 0, cscale, anrm, n, 1, wi, n, ierr )
491  END IF
492  DO 10 i = 1, n
493  bwork( i ) = SELECT( wr( i ), wi( i ) )
494  10 CONTINUE
495 *
496 * Reorder eigenvalues, transform Schur vectors, and compute
497 * reciprocal condition numbers
498 * (RWorkspace: if SENSE is not 'N', need N+2*SDIM*(N-SDIM)
499 * otherwise, need N )
500 * (IWorkspace: if SENSE is 'V' or 'B', need SDIM*(N-SDIM)
501 * otherwise, need 0 )
502 *
503  CALL dtrsen( sense, jobvs, bwork, n, a, lda, vs, ldvs, wr, wi,
504  \$ sdim, rconde, rcondv, work( iwrk ), lwork-iwrk+1,
505  \$ iwork, liwork, icond )
506  IF( .NOT.wantsn )
507  \$ maxwrk = max( maxwrk, n+2*sdim*( n-sdim ) )
508  IF( icond.EQ.-15 ) THEN
509 *
510 * Not enough real workspace
511 *
512  info = -16
513  ELSE IF( icond.EQ.-17 ) THEN
514 *
515 * Not enough integer workspace
516 *
517  info = -18
518  ELSE IF( icond.GT.0 ) THEN
519 *
520 * DTRSEN failed to reorder or to restore standard Schur form
521 *
522  info = icond + n
523  END IF
524  END IF
525 *
526  IF( wantvs ) THEN
527 *
528 * Undo balancing
529 * (RWorkspace: need N)
530 *
531  CALL dgebak( 'P', 'R', n, ilo, ihi, work( ibal ), n, vs, ldvs,
532  \$ ierr )
533  END IF
534 *
535  IF( scalea ) THEN
536 *
537 * Undo scaling for the Schur form of A
538 *
539  CALL dlascl( 'H', 0, 0, cscale, anrm, n, n, a, lda, ierr )
540  CALL dcopy( n, a, lda+1, wr, 1 )
541  IF( ( wantsv .OR. wantsb ) .AND. info.EQ.0 ) THEN
542  dum( 1 ) = rcondv
543  CALL dlascl( 'G', 0, 0, cscale, anrm, 1, 1, dum, 1, ierr )
544  rcondv = dum( 1 )
545  END IF
546  IF( cscale.EQ.smlnum ) THEN
547 *
548 * If scaling back towards underflow, adjust WI if an
549 * offdiagonal element of a 2-by-2 block in the Schur form
550 * underflows.
551 *
552  IF( ieval.GT.0 ) THEN
553  i1 = ieval + 1
554  i2 = ihi - 1
555  CALL dlascl( 'G', 0, 0, cscale, anrm, ilo-1, 1, wi, n,
556  \$ ierr )
557  ELSE IF( wantst ) THEN
558  i1 = 1
559  i2 = n - 1
560  ELSE
561  i1 = ilo
562  i2 = ihi - 1
563  END IF
564  inxt = i1 - 1
565  DO 20 i = i1, i2
566  IF( i.LT.inxt )
567  \$ GO TO 20
568  IF( wi( i ).EQ.zero ) THEN
569  inxt = i + 1
570  ELSE
571  IF( a( i+1, i ).EQ.zero ) THEN
572  wi( i ) = zero
573  wi( i+1 ) = zero
574  ELSE IF( a( i+1, i ).NE.zero .AND. a( i, i+1 ).EQ.
575  \$ zero ) THEN
576  wi( i ) = zero
577  wi( i+1 ) = zero
578  IF( i.GT.1 )
579  \$ CALL dswap( i-1, a( 1, i ), 1, a( 1, i+1 ), 1 )
580  IF( n.GT.i+1 )
581  \$ CALL dswap( n-i-1, a( i, i+2 ), lda,
582  \$ a( i+1, i+2 ), lda )
583  IF( wantvs ) THEN
584  CALL dswap( n, vs( 1, i ), 1, vs( 1, i+1 ), 1 )
585  END IF
586  a( i, i+1 ) = a( i+1, i )
587  a( i+1, i ) = zero
588  END IF
589  inxt = i + 2
590  END IF
591  20 CONTINUE
592  END IF
593  CALL dlascl( 'G', 0, 0, cscale, anrm, n-ieval, 1,
594  \$ wi( ieval+1 ), max( n-ieval, 1 ), ierr )
595  END IF
596 *
597  IF( wantst .AND. info.EQ.0 ) THEN
598 *
599 * Check if reordering successful
600 *
601  lastsl = .true.
602  lst2sl = .true.
603  sdim = 0
604  ip = 0
605  DO 30 i = 1, n
606  cursl = SELECT( wr( i ), wi( i ) )
607  IF( wi( i ).EQ.zero ) THEN
608  IF( cursl )
609  \$ sdim = sdim + 1
610  ip = 0
611  IF( cursl .AND. .NOT.lastsl )
612  \$ info = n + 2
613  ELSE
614  IF( ip.EQ.1 ) THEN
615 *
616 * Last eigenvalue of conjugate pair
617 *
618  cursl = cursl .OR. lastsl
619  lastsl = cursl
620  IF( cursl )
621  \$ sdim = sdim + 2
622  ip = -1
623  IF( cursl .AND. .NOT.lst2sl )
624  \$ info = n + 2
625  ELSE
626 *
627 * First eigenvalue of conjugate pair
628 *
629  ip = 1
630  END IF
631  END IF
632  lst2sl = lastsl
633  lastsl = cursl
634  30 CONTINUE
635  END IF
636 *
637  work( 1 ) = maxwrk
638  IF( wantsv .OR. wantsb ) THEN
639  iwork( 1 ) = max( 1, sdim*( n-sdim ) )
640  ELSE
641  iwork( 1 ) = 1
642  END IF
643 *
644  RETURN
645 *
646 * End of DGEESX
647 *
double precision function dlamch(CMACH)
DLAMCH
Definition: dlamch.f:69
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:74
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:143
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:103
integer function ilaenv(ISPEC, NAME, OPTS, N1, N2, N3, N4)
ILAENV
Definition: ilaenv.f:162
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
logical function lsame(CA, CB)
LSAME
Definition: lsame.f:53
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82
subroutine dswap(N, DX, INCX, DY, INCY)
DSWAP
Definition: dswap.f:82
double precision function dlange(NORM, M, N, A, LDA, WORK)
DLANGE returns the value of the 1-norm, Frobenius norm, infinity-norm, or the largest absolute value ...
Definition: dlange.f:114
subroutine dgehrd(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
DGEHRD
Definition: dgehrd.f:167
subroutine dgebal(JOB, N, A, LDA, ILO, IHI, SCALE, INFO)
DGEBAL
Definition: dgebal.f:160
subroutine dgebak(JOB, SIDE, N, ILO, IHI, SCALE, M, V, LDV, INFO)
DGEBAK
Definition: dgebak.f:130
subroutine dhseqr(JOB, COMPZ, N, ILO, IHI, H, LDH, WR, WI, Z, LDZ, WORK, LWORK, INFO)
DHSEQR
Definition: dhseqr.f:316
subroutine dorghr(N, ILO, IHI, A, LDA, TAU, WORK, LWORK, INFO)
DORGHR
Definition: dorghr.f:126
subroutine dtrsen(JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, WR, WI, M, S, SEP, WORK, LWORK, IWORK, LIWORK, INFO)
DTRSEN
Definition: dtrsen.f:313
Here is the call graph for this function:
Here is the caller graph for this function: