SRC/zgeesx.f

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