SRC/zgeev.f

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