SRC/dgees.f

   1 *> \brief <b> DGEES 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 DGEES + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgees.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgees.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgees.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE DGEES( JOBVS, SORT, SELECT, N, A, LDA, SDIM, WR, WI,
  22 *                         VS, LDVS, WORK, LWORK, 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 *       DOUBLE PRECISION   A( LDA, * ), VS( LDVS, * ), WI( * ), WORK( * ),
  31 *      $                   WR( * )
  32 *       ..
  33 *       .. Function Arguments ..
  34 *       LOGICAL            SELECT
  35 *       EXTERNAL           SELECT
  36 *       ..
  37 *
  38 *
  39 *> \par Purpose:
  40 *  =============
  41 *>
  42 *> \verbatim
  43 *>
  44 *> DGEES computes for an N-by-N real nonsymmetric matrix A, the
  45 *> eigenvalues, the real Schur form T, and, optionally, the matrix of
  46 *> Schur vectors Z.  This gives the Schur factorization A = Z*T*(Z**T).
  47 *>
  48 *> Optionally, it also orders the eigenvalues on the diagonal of the
  49 *> real 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 matrix is in real Schur form if it is upper quasi-triangular with
  54 *> 1-by-1 and 2-by-2 blocks. 2-by-2 blocks will be standardized in the
  55 *> form
  56 *>         [  a  b  ]
  57 *>         [  c  a  ]
  58 *>
  59 *> where b*c < 0. The eigenvalues of such a block are a +- sqrt(bc).
  60 *> \endverbatim
  61 *
  62 *  Arguments:
  63 *  ==========
  64 *
  65 *> \param[in] JOBVS
  66 *> \verbatim
  67 *>          JOBVS is CHARACTER*1
  68 *>          = 'N': Schur vectors are not computed;
  69 *>          = 'V': Schur vectors are computed.
  70 *> \endverbatim
  71 *>
  72 *> \param[in] SORT
  73 *> \verbatim
  74 *>          SORT is CHARACTER*1
  75 *>          Specifies whether or not to order the eigenvalues on the
  76 *>          diagonal of the Schur form.
  77 *>          = 'N': Eigenvalues are not ordered;
  78 *>          = 'S': Eigenvalues are ordered (see SELECT).
  79 *> \endverbatim
  80 *>
  81 *> \param[in] SELECT
  82 *> \verbatim
  83 *>          SELECT is a LOGICAL FUNCTION of two DOUBLE PRECISION arguments
  84 *>          SELECT must be declared EXTERNAL in the calling subroutine.
  85 *>          If SORT = 'S', SELECT is used to select eigenvalues to sort
  86 *>          to the top left of the Schur form.
  87 *>          If SORT = 'N', SELECT is not referenced.
  88 *>          An eigenvalue WR(j)+sqrt(-1)*WI(j) is selected if
  89 *>          SELECT(WR(j),WI(j)) is true; i.e., if either one of a complex
  90 *>          conjugate pair of eigenvalues is selected, then both complex
  91 *>          eigenvalues are selected.
  92 *>          Note that a selected complex eigenvalue may no longer
  93 *>          satisfy SELECT(WR(j),WI(j)) = .TRUE. after ordering, since
  94 *>          ordering may change the value of complex eigenvalues
  95 *>          (especially if the eigenvalue is ill-conditioned); in this
  96 *>          case INFO is set to N+2 (see INFO below).
  97 *> \endverbatim
  98 *>
  99 *> \param[in] N
 100 *> \verbatim
 101 *>          N is INTEGER
 102 *>          The order of the matrix A. N >= 0.
 103 *> \endverbatim
 104 *>
 105 *> \param[in,out] A
 106 *> \verbatim
 107 *>          A is DOUBLE PRECISION array, dimension (LDA,N)
 108 *>          On entry, the N-by-N matrix A.
 109 *>          On exit, A has been overwritten by its real Schur form T.
 110 *> \endverbatim
 111 *>
 112 *> \param[in] LDA
 113 *> \verbatim
 114 *>          LDA is INTEGER
 115 *>          The leading dimension of the array A.  LDA >= max(1,N).
 116 *> \endverbatim
 117 *>
 118 *> \param[out] SDIM
 119 *> \verbatim
 120 *>          SDIM is INTEGER
 121 *>          If SORT = 'N', SDIM = 0.
 122 *>          If SORT = 'S', SDIM = number of eigenvalues (after sorting)
 123 *>                         for which SELECT is true. (Complex conjugate
 124 *>                         pairs for which SELECT is true for either
 125 *>                         eigenvalue count as 2.)
 126 *> \endverbatim
 127 *>
 128 *> \param[out] WR
 129 *> \verbatim
 130 *>          WR is DOUBLE PRECISION array, dimension (N)
 131 *> \endverbatim
 132 *>
 133 *> \param[out] WI
 134 *> \verbatim
 135 *>          WI is DOUBLE PRECISION array, dimension (N)
 136 *>          WR and WI contain the real and imaginary parts,
 137 *>          respectively, of the computed eigenvalues in the same order
 138 *>          that they appear on the diagonal of the output Schur form T.
 139 *>          Complex conjugate pairs of eigenvalues will appear
 140 *>          consecutively with the eigenvalue having the positive
 141 *>          imaginary part first.
 142 *> \endverbatim
 143 *>
 144 *> \param[out] VS
 145 *> \verbatim
 146 *>          VS is DOUBLE PRECISION array, dimension (LDVS,N)
 147 *>          If JOBVS = 'V', VS contains the orthogonal matrix Z of Schur
 148 *>          vectors.
 149 *>          If JOBVS = 'N', VS is not referenced.
 150 *> \endverbatim
 151 *>
 152 *> \param[in] LDVS
 153 *> \verbatim
 154 *>          LDVS is INTEGER
 155 *>          The leading dimension of the array VS.  LDVS >= 1; if
 156 *>          JOBVS = 'V', LDVS >= N.
 157 *> \endverbatim
 158 *>
 159 *> \param[out] WORK
 160 *> \verbatim
 161 *>          WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
 162 *>          On exit, if INFO = 0, WORK(1) contains the optimal LWORK.
 163 *> \endverbatim
 164 *>
 165 *> \param[in] LWORK
 166 *> \verbatim
 167 *>          LWORK is INTEGER
 168 *>          The dimension of the array WORK.  LWORK >= max(1,3*N).
 169 *>          For good performance, LWORK must generally be larger.
 170 *>
 171 *>          If LWORK = -1, then a workspace query is assumed; the routine
 172 *>          only calculates the optimal size of the WORK array, returns
 173 *>          this value as the first entry of the WORK array, and no error
 174 *>          message related to LWORK is issued by XERBLA.
 175 *> \endverbatim
 176 *>
 177 *> \param[out] BWORK
 178 *> \verbatim
 179 *>          BWORK is LOGICAL array, dimension (N)
 180 *>          Not referenced if SORT = 'N'.
 181 *> \endverbatim
 182 *>
 183 *> \param[out] INFO
 184 *> \verbatim
 185 *>          INFO is INTEGER
 186 *>          = 0: successful exit
 187 *>          < 0: if INFO = -i, the i-th argument had an illegal value.
 188 *>          > 0: if INFO = i, and i is
 189 *>             <= N: the QR algorithm failed to compute all the
 190 *>                   eigenvalues; elements 1:ILO-1 and i+1:N of WR and WI
 191 *>                   contain those eigenvalues which have converged; if
 192 *>                   JOBVS = 'V', VS contains the matrix which reduces A
 193 *>                   to its partially converged Schur form.
 194 *>             = N+1: the eigenvalues could not be reordered because some
 195 *>                   eigenvalues were too close to separate (the problem
 196 *>                   is very ill-conditioned);
 197 *>             = N+2: after reordering, roundoff changed values of some
 198 *>                   complex eigenvalues so that leading eigenvalues in
 199 *>                   the Schur form no longer satisfy SELECT=.TRUE.  This
 200 *>                   could also be caused by underflow due to scaling.
 201 *> \endverbatim
 202 *
 203 *  Authors:
 204 *  ========
 205 *
 206 *> \author Univ. of Tennessee
 207 *> \author Univ. of California Berkeley
 208 *> \author Univ. of Colorado Denver
 209 *> \author NAG Ltd.
 210 *
 211 *> \date November 2011
 212 *
 213 *> \ingroup doubleGEeigen
 214 *
 215 *  =====================================================================
 216       SUBROUTINE DGEES( JOBVS, SORT, SELECT, N, A, LDA, SDIM, WR, WI,
 217      $                  VS, LDVS, WORK, LWORK, BWORK, INFO )
 218 *
 219 *  -- LAPACK driver routine (version 3.4.0) --
 220 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 221 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
 222 *     November 2011
 223 *
 224 *     .. Scalar Arguments ..
 225       CHARACTER          JOBVS, SORT
 226       INTEGER            INFO, LDA, LDVS, LWORK, N, SDIM
 227 *     ..
 228 *     .. Array Arguments ..
 229       LOGICAL            BWORK( * )
 230       DOUBLE PRECISION   A( LDA, * ), VS( LDVS, * ), WI( * ), WORK( * ),
 231      $                   WR( * )
 232 *     ..
 233 *     .. Function Arguments ..
 234       LOGICAL            SELECT
 235       EXTERNAL           SELECT
 236 *     ..
 237 *
 238 *  =====================================================================
 239 *
 240 *     .. Parameters ..
 241       DOUBLE PRECISION   ZERO, ONE
 242       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
 243 *     ..
 244 *     .. Local Scalars ..
 245       LOGICAL            CURSL, LASTSL, LQUERY, LST2SL, SCALEA, WANTST,
 246      $                   WANTVS
 247       INTEGER            HSWORK, I, I1, I2, IBAL, ICOND, IERR, IEVAL,
 248      $                   IHI, ILO, INXT, IP, ITAU, IWRK, MAXWRK, MINWRK
 249       DOUBLE PRECISION   ANRM, BIGNUM, CSCALE, EPS, S, SEP, SMLNUM
 250 *     ..
 251 *     .. Local Arrays ..
 252       INTEGER            IDUM( 1 )
 253       DOUBLE PRECISION   DUM( 1 )
 254 *     ..
 255 *     .. External Subroutines ..
 256       EXTERNAL           DCOPY, DGEBAK, DGEBAL, DGEHRD, DHSEQR, DLACPY,
 257      $                   DLABAD, DLASCL, DORGHR, DSWAP, DTRSEN, XERBLA
 258 *     ..
 259 *     .. External Functions ..
 260       LOGICAL            LSAME
 261       INTEGER            ILAENV
 262       DOUBLE PRECISION   DLAMCH, DLANGE
 263       EXTERNAL           LSAME, ILAENV, DLAMCH, DLANGE
 264 *     ..
 265 *     .. Intrinsic Functions ..
 266       INTRINSIC          MAX, SQRT
 267 *     ..
 268 *     .. Executable Statements ..
 269 *
 270 *     Test the input arguments
 271 *
 272       INFO = 0
 273       LQUERY = ( LWORK.EQ.-1 )
 274       WANTVS = LSAME( JOBVS, 'V' )
 275       WANTST = LSAME( SORT, 'S' )
 276       IF( ( .NOT.WANTVS ) .AND. ( .NOT.LSAME( JOBVS, 'N' ) ) ) THEN
 277          INFO = -1
 278       ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN
 279          INFO = -2
 280       ELSE IF( N.LT.0 ) THEN
 281          INFO = -4
 282       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
 283          INFO = -6
 284       ELSE IF( LDVS.LT.1 .OR. ( WANTVS .AND. LDVS.LT.N ) ) THEN
 285          INFO = -11
 286       END IF
 287 *
 288 *     Compute workspace
 289 *      (Note: Comments in the code beginning "Workspace:" describe the
 290 *       minimal amount of workspace needed at that point in the code,
 291 *       as well as the preferred amount for good performance.
 292 *       NB refers to the optimal block size for the immediately
 293 *       following subroutine, as returned by ILAENV.
 294 *       HSWORK refers to the workspace preferred by DHSEQR, as
 295 *       calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
 296 *       the worst case.)
 297 *
 298       IF( INFO.EQ.0 ) THEN
 299          IF( N.EQ.0 ) THEN
 300             MINWRK = 1
 301             MAXWRK = 1
 302          ELSE
 303             MAXWRK = 2*N + N*ILAENV( 1, 'DGEHRD', ' ', N, 1, N, 0 )
 304             MINWRK = 3*N
 305 *
 306             CALL DHSEQR( 'S', JOBVS, N, 1, N, A, LDA, WR, WI, VS, LDVS,
 307      $             WORK, -1, IEVAL )
 308             HSWORK = WORK( 1 )
 309 *
 310             IF( .NOT.WANTVS ) THEN
 311                MAXWRK = MAX( MAXWRK, N + HSWORK )
 312             ELSE
 313                MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1,
 314      $                       'DORGHR', ' ', N, 1, N, -1 ) )
 315                MAXWRK = MAX( MAXWRK, N + HSWORK )
 316             END IF
 317          END IF
 318          WORK( 1 ) = MAXWRK
 319 *
 320          IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
 321             INFO = -13
 322          END IF
 323       END IF
 324 *
 325       IF( INFO.NE.0 ) THEN
 326          CALL XERBLA( 'DGEES ', -INFO )
 327          RETURN
 328       ELSE IF( LQUERY ) THEN
 329          RETURN
 330       END IF
 331 *
 332 *     Quick return if possible
 333 *
 334       IF( N.EQ.0 ) THEN
 335          SDIM = 0
 336          RETURN
 337       END IF
 338 *
 339 *     Get machine constants
 340 *
 341       EPS = DLAMCH( 'P' )
 342       SMLNUM = DLAMCH( 'S' )
 343       BIGNUM = ONE / SMLNUM
 344       CALL DLABAD( SMLNUM, BIGNUM )
 345       SMLNUM = SQRT( SMLNUM ) / EPS
 346       BIGNUM = ONE / SMLNUM
 347 *
 348 *     Scale A if max element outside range [SMLNUM,BIGNUM]
 349 *
 350       ANRM = DLANGE( 'M', N, N, A, LDA, DUM )
 351       SCALEA = .FALSE.
 352       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
 353          SCALEA = .TRUE.
 354          CSCALE = SMLNUM
 355       ELSE IF( ANRM.GT.BIGNUM ) THEN
 356          SCALEA = .TRUE.
 357          CSCALE = BIGNUM
 358       END IF
 359       IF( SCALEA )
 360      $   CALL DLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR )
 361 *
 362 *     Permute the matrix to make it more nearly triangular
 363 *     (Workspace: need N)
 364 *
 365       IBAL = 1
 366       CALL DGEBAL( 'P', N, A, LDA, ILO, IHI, WORK( IBAL ), IERR )
 367 *
 368 *     Reduce to upper Hessenberg form
 369 *     (Workspace: need 3*N, prefer 2*N+N*NB)
 370 *
 371       ITAU = N + IBAL
 372       IWRK = N + ITAU
 373       CALL DGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ),
 374      $             LWORK-IWRK+1, IERR )
 375 *
 376       IF( WANTVS ) THEN
 377 *
 378 *        Copy Householder vectors to VS
 379 *
 380          CALL DLACPY( 'L', N, N, A, LDA, VS, LDVS )
 381 *
 382 *        Generate orthogonal matrix in VS
 383 *        (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB)
 384 *
 385          CALL DORGHR( N, ILO, IHI, VS, LDVS, WORK( ITAU ), WORK( IWRK ),
 386      $                LWORK-IWRK+1, IERR )
 387       END IF
 388 *
 389       SDIM = 0
 390 *
 391 *     Perform QR iteration, accumulating Schur vectors in VS if desired
 392 *     (Workspace: need N+1, prefer N+HSWORK (see comments) )
 393 *
 394       IWRK = ITAU
 395       CALL DHSEQR( 'S', JOBVS, N, ILO, IHI, A, LDA, WR, WI, VS, LDVS,
 396      $             WORK( IWRK ), LWORK-IWRK+1, IEVAL )
 397       IF( IEVAL.GT.0 )
 398      $   INFO = IEVAL
 399 *
 400 *     Sort eigenvalues if desired
 401 *
 402       IF( WANTST .AND. INFO.EQ.0 ) THEN
 403          IF( SCALEA ) THEN
 404             CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N, 1, WR, N, IERR )
 405             CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N, 1, WI, N, IERR )
 406          END IF
 407          DO 10 I = 1, N
 408             BWORK( I ) = SELECT( WR( I ), WI( I ) )
 409    10    CONTINUE
 410 *
 411 *        Reorder eigenvalues and transform Schur vectors
 412 *        (Workspace: none needed)
 413 *
 414          CALL DTRSEN( 'N', JOBVS, BWORK, N, A, LDA, VS, LDVS, WR, WI,
 415      $                SDIM, S, SEP, WORK( IWRK ), LWORK-IWRK+1, IDUM, 1,
 416      $                ICOND )
 417          IF( ICOND.GT.0 )
 418      $      INFO = N + ICOND
 419       END IF
 420 *
 421       IF( WANTVS ) THEN
 422 *
 423 *        Undo balancing
 424 *        (Workspace: need N)
 425 *
 426          CALL DGEBAK( 'P', 'R', N, ILO, IHI, WORK( IBAL ), N, VS, LDVS,
 427      $                IERR )
 428       END IF
 429 *
 430       IF( SCALEA ) THEN
 431 *
 432 *        Undo scaling for the Schur form of A
 433 *
 434          CALL DLASCL( 'H', 0, 0, CSCALE, ANRM, N, N, A, LDA, IERR )
 435          CALL DCOPY( N, A, LDA+1, WR, 1 )
 436          IF( CSCALE.EQ.SMLNUM ) THEN
 437 *
 438 *           If scaling back towards underflow, adjust WI if an
 439 *           offdiagonal element of a 2-by-2 block in the Schur form
 440 *           underflows.
 441 *
 442             IF( IEVAL.GT.0 ) THEN
 443                I1 = IEVAL + 1
 444                I2 = IHI - 1
 445                CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WI,
 446      $                      MAX( ILO-1, 1 ), IERR )
 447             ELSE IF( WANTST ) THEN
 448                I1 = 1
 449                I2 = N - 1
 450             ELSE
 451                I1 = ILO
 452                I2 = IHI - 1
 453             END IF
 454             INXT = I1 - 1
 455             DO 20 I = I1, I2
 456                IF( I.LT.INXT )
 457      $            GO TO 20
 458                IF( WI( I ).EQ.ZERO ) THEN
 459                   INXT = I + 1
 460                ELSE
 461                   IF( A( I+1, I ).EQ.ZERO ) THEN
 462                      WI( I ) = ZERO
 463                      WI( I+1 ) = ZERO
 464                   ELSE IF( A( I+1, I ).NE.ZERO .AND. A( I, I+1 ).EQ.
 465      $                     ZERO ) THEN
 466                      WI( I ) = ZERO
 467                      WI( I+1 ) = ZERO
 468                      IF( I.GT.1 )
 469      $                  CALL DSWAP( I-1, A( 1, I ), 1, A( 1, I+1 ), 1 )
 470                      IF( N.GT.I+1 )
 471      $                  CALL DSWAP( N-I-1, A( I, I+2 ), LDA,
 472      $                              A( I+1, I+2 ), LDA )
 473                      IF( WANTVS ) THEN
 474                         CALL DSWAP( N, VS( 1, I ), 1, VS( 1, I+1 ), 1 )
 475                      END IF
 476                      A( I, I+1 ) = A( I+1, I )
 477                      A( I+1, I ) = ZERO
 478                   END IF
 479                   INXT = I + 2
 480                END IF
 481    20       CONTINUE
 482          END IF
 483 *
 484 *        Undo scaling for the imaginary part of the eigenvalues
 485 *
 486          CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N-IEVAL, 1,
 487      $                WI( IEVAL+1 ), MAX( N-IEVAL, 1 ), IERR )
 488       END IF
 489 *
 490       IF( WANTST .AND. INFO.EQ.0 ) THEN
 491 *
 492 *        Check if reordering successful
 493 *
 494          LASTSL = .TRUE.
 495          LST2SL = .TRUE.
 496          SDIM = 0
 497          IP = 0
 498          DO 30 I = 1, N
 499             CURSL = SELECT( WR( I ), WI( I ) )
 500             IF( WI( I ).EQ.ZERO ) THEN
 501                IF( CURSL )
 502      $            SDIM = SDIM + 1
 503                IP = 0
 504                IF( CURSL .AND. .NOT.LASTSL )
 505      $            INFO = N + 2
 506             ELSE
 507                IF( IP.EQ.1 ) THEN
 508 *
 509 *                 Last eigenvalue of conjugate pair
 510 *
 511                   CURSL = CURSL .OR. LASTSL
 512                   LASTSL = CURSL
 513                   IF( CURSL )
 514      $               SDIM = SDIM + 2
 515                   IP = -1
 516                   IF( CURSL .AND. .NOT.LST2SL )
 517      $               INFO = N + 2
 518                ELSE
 519 *
 520 *                 First eigenvalue of conjugate pair
 521 *
 522                   IP = 1
 523                END IF
 524             END IF
 525             LST2SL = LASTSL
 526             LASTSL = CURSL
 527    30    CONTINUE
 528       END IF
 529 *
 530       WORK( 1 ) = MAXWRK
 531       RETURN
 532 *
 533 *     End of DGEES
 534 *
 535       END