SRC/zgees.f

   1 *> \brief <b> ZGEES 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 ZGEES + dependencies
  10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgees.f">
  11 *> [TGZ]</a>
  12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgees.f">
  13 *> [ZIP]</a>
  14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgees.f">
  15 *> [TXT]</a>
  16 *> \endhtmlonly
  17 *
  18 *  Definition:
  19 *  ===========
  20 *
  21 *       SUBROUTINE ZGEES( 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 *       DOUBLE PRECISION   RWORK( * )
  31 *       COMPLEX*16         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 *> ZGEES 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*16 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*16 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*16 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*16 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*16 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 DOUBLE PRECISION 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 complex16GEeigen
 195 *
 196 *  =====================================================================
 197       SUBROUTINE ZGEES( 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       DOUBLE PRECISION   RWORK( * )
 212       COMPLEX*16         A( LDA, * ), VS( LDVS, * ), W( * ), WORK( * )
 213 *     ..
 214 *     .. Function Arguments ..
 215       LOGICAL            SELECT
 216       EXTERNAL           SELECT
 217 *     ..
 218 *
 219 *  =====================================================================
 220 *
 221 *     .. Parameters ..
 222       DOUBLE PRECISION   ZERO, ONE
 223       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
 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       DOUBLE PRECISION   ANRM, BIGNUM, CSCALE, EPS, S, SEP, SMLNUM
 230 *     ..
 231 *     .. Local Arrays ..
 232       DOUBLE PRECISION   DUM( 1 )
 233 *     ..
 234 *     .. External Subroutines ..
 235       EXTERNAL           DLABAD, XERBLA, ZCOPY, ZGEBAK, ZGEBAL, ZGEHRD,
 236      $                   ZHSEQR, ZLACPY, ZLASCL, ZTRSEN, ZUNGHR
 237 *     ..
 238 *     .. External Functions ..
 239       LOGICAL            LSAME
 240       INTEGER            ILAENV
 241       DOUBLE PRECISION   DLAMCH, ZLANGE
 242       EXTERNAL           LSAME, ILAENV, DLAMCH, ZLANGE
 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 ZHSEQR, 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, 'ZGEHRD', ' ', N, 1, N, 0 )
 284             MINWRK = 2*N
 285 *
 286             CALL ZHSEQR( '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, 'ZUNGHR',
 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( 'ZGEES ', -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 = DLAMCH( 'P' )
 322       SMLNUM = DLAMCH( 'S' )
 323       BIGNUM = ONE / SMLNUM
 324       CALL DLABAD( 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 = ZLANGE( '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 ZLASCL( '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 ZGEBAL( '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 ZGEHRD( 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 ZLACPY( '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 ZUNGHR( 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 ZHSEQR( '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 ZLASCL( '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 ZTRSEN( '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 ZGEBAK( '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 ZLASCL( 'U', 0, 0, CSCALE, ANRM, N, N, A, LDA, IERR )
 416          CALL ZCOPY( N, A, LDA+1, W, 1 )
 417       END IF
 418 *
 419       WORK( 1 ) = MAXWRK
 420       RETURN
 421 *
 422 *     End of ZGEES
 423 *
 424       END