LAPACK: SRC/zgeev.f Source File

Go to the documentation of this file.
00001       SUBROUTINE ZGEEV( JOBVL, JOBVR, N, A, LDA, W, VL, LDVL, VR, LDVR,
00002      $                  WORK, LWORK, RWORK, INFO )
00003 *
00004 *  -- LAPACK driver routine (version 3.2) --
00005 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00006 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00007 *     November 2006
00008 *
00009 *     .. Scalar Arguments ..
00010       CHARACTER          JOBVL, JOBVR
00011       INTEGER            INFO, LDA, LDVL, LDVR, LWORK, N
00012 *     ..
00013 *     .. Array Arguments ..
00014       DOUBLE PRECISION   RWORK( * )
00015       COMPLEX*16         A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ),
00016      $                   W( * ), WORK( * )
00017 *     ..
00018 *
00019 *  Purpose
00020 *  =======
00021 *
00022 *  ZGEEV computes for an N-by-N complex nonsymmetric matrix A, the
00023 *  eigenvalues and, optionally, the left and/or right eigenvectors.
00024 *
00025 *  The right eigenvector v(j) of A satisfies
00026 *                   A * v(j) = lambda(j) * v(j)
00027 *  where lambda(j) is its eigenvalue.
00028 *  The left eigenvector u(j) of A satisfies
00029 *                u(j)**H * A = lambda(j) * u(j)**H
00030 *  where u(j)**H denotes the conjugate transpose of u(j).
00031 *
00032 *  The computed eigenvectors are normalized to have Euclidean norm
00033 *  equal to 1 and largest component real.
00034 *
00035 *  Arguments
00036 *  =========
00037 *
00038 *  JOBVL   (input) CHARACTER*1
00039 *          = 'N': left eigenvectors of A are not computed;
00040 *          = 'V': left eigenvectors of are computed.
00041 *
00042 *  JOBVR   (input) CHARACTER*1
00043 *          = 'N': right eigenvectors of A are not computed;
00044 *          = 'V': right eigenvectors of A are computed.
00045 *
00046 *  N       (input) INTEGER
00047 *          The order of the matrix A. N >= 0.
00048 *
00049 *  A       (input/output) COMPLEX*16 array, dimension (LDA,N)
00050 *          On entry, the N-by-N matrix A.
00051 *          On exit, A has been overwritten.
00052 *
00053 *  LDA     (input) INTEGER
00054 *          The leading dimension of the array A.  LDA >= max(1,N).
00055 *
00056 *  W       (output) COMPLEX*16 array, dimension (N)
00057 *          W contains the computed eigenvalues.
00058 *
00059 *  VL      (output) COMPLEX*16 array, dimension (LDVL,N)
00060 *          If JOBVL = 'V', the left eigenvectors u(j) are stored one
00061 *          after another in the columns of VL, in the same order
00062 *          as their eigenvalues.
00063 *          If JOBVL = 'N', VL is not referenced.
00064 *          u(j) = VL(:,j), the j-th column of VL.
00065 *
00066 *  LDVL    (input) INTEGER
00067 *          The leading dimension of the array VL.  LDVL >= 1; if
00068 *          JOBVL = 'V', LDVL >= N.
00069 *
00070 *  VR      (output) COMPLEX*16 array, dimension (LDVR,N)
00071 *          If JOBVR = 'V', the right eigenvectors v(j) are stored one
00072 *          after another in the columns of VR, in the same order
00073 *          as their eigenvalues.
00074 *          If JOBVR = 'N', VR is not referenced.
00075 *          v(j) = VR(:,j), the j-th column of VR.
00076 *
00077 *  LDVR    (input) INTEGER
00078 *          The leading dimension of the array VR.  LDVR >= 1; if
00079 *          JOBVR = 'V', LDVR >= N.
00080 *
00081 *  WORK    (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
00082 *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00083 *
00084 *  LWORK   (input) INTEGER
00085 *          The dimension of the array WORK.  LWORK >= max(1,2*N).
00086 *          For good performance, LWORK must generally be larger.
00087 *
00088 *          If LWORK = -1, then a workspace query is assumed; the routine
00089 *          only calculates the optimal size of the WORK array, returns
00090 *          this value as the first entry of the WORK array, and no error
00091 *          message related to LWORK is issued by XERBLA.
00092 *
00093 *  RWORK   (workspace) DOUBLE PRECISION array, dimension (2*N)
00094 *
00095 *  INFO    (output) INTEGER
00096 *          = 0:  successful exit
00097 *          < 0:  if INFO = -i, the i-th argument had an illegal value.
00098 *          > 0:  if INFO = i, the QR algorithm failed to compute all the
00099 *                eigenvalues, and no eigenvectors have been computed;
00100 *                elements and i+1:N of W contain eigenvalues which have
00101 *                converged.
00102 *
00103 *  =====================================================================
00104 *
00105 *     .. Parameters ..
00106       DOUBLE PRECISION   ZERO, ONE
00107       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
00108 *     ..
00109 *     .. Local Scalars ..
00110       LOGICAL            LQUERY, SCALEA, WANTVL, WANTVR
00111       CHARACTER          SIDE
00112       INTEGER            HSWORK, I, IBAL, IERR, IHI, ILO, IRWORK, ITAU,
00113      $                   IWRK, K, MAXWRK, MINWRK, NOUT
00114       DOUBLE PRECISION   ANRM, BIGNUM, CSCALE, EPS, SCL, SMLNUM
00115       COMPLEX*16         TMP
00116 *     ..
00117 *     .. Local Arrays ..
00118       LOGICAL            SELECT( 1 )
00119       DOUBLE PRECISION   DUM( 1 )
00120 *     ..
00121 *     .. External Subroutines ..
00122       EXTERNAL           DLABAD, XERBLA, ZDSCAL, ZGEBAK, ZGEBAL, ZGEHRD,
00123      $                   ZHSEQR, ZLACPY, ZLASCL, ZSCAL, ZTREVC, ZUNGHR
00124 *     ..
00125 *     .. External Functions ..
00126       LOGICAL            LSAME
00127       INTEGER            IDAMAX, ILAENV
00128       DOUBLE PRECISION   DLAMCH, DZNRM2, ZLANGE
00129       EXTERNAL           LSAME, IDAMAX, ILAENV, DLAMCH, DZNRM2, ZLANGE
00130 *     ..
00131 *     .. Intrinsic Functions ..
00132       INTRINSIC          DBLE, DCMPLX, DCONJG, DIMAG, MAX, SQRT
00133 *     ..
00134 *     .. Executable Statements ..
00135 *
00136 *     Test the input arguments
00137 *
00138       INFO = 0
00139       LQUERY = ( LWORK.EQ.-1 )
00140       WANTVL = LSAME( JOBVL, 'V' )
00141       WANTVR = LSAME( JOBVR, 'V' )
00142       IF( ( .NOT.WANTVL ) .AND. ( .NOT.LSAME( JOBVL, 'N' ) ) ) THEN
00143          INFO = -1
00144       ELSE IF( ( .NOT.WANTVR ) .AND. ( .NOT.LSAME( JOBVR, 'N' ) ) ) THEN
00145          INFO = -2
00146       ELSE IF( N.LT.0 ) THEN
00147          INFO = -3
00148       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00149          INFO = -5
00150       ELSE IF( LDVL.LT.1 .OR. ( WANTVL .AND. LDVL.LT.N ) ) THEN
00151          INFO = -8
00152       ELSE IF( LDVR.LT.1 .OR. ( WANTVR .AND. LDVR.LT.N ) ) THEN
00153          INFO = -10
00154       END IF
00155 *
00156 *     Compute workspace
00157 *      (Note: Comments in the code beginning "Workspace:" describe the
00158 *       minimal amount of workspace needed at that point in the code,
00159 *       as well as the preferred amount for good performance.
00160 *       CWorkspace refers to complex workspace, and RWorkspace to real
00161 *       workspace. NB refers to the optimal block size for the
00162 *       immediately following subroutine, as returned by ILAENV.
00163 *       HSWORK refers to the workspace preferred by ZHSEQR, as
00164 *       calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
00165 *       the worst case.)
00166 *
00167       IF( INFO.EQ.0 ) THEN
00168          IF( N.EQ.0 ) THEN
00169             MINWRK = 1
00170             MAXWRK = 1
00171          ELSE
00172             MAXWRK = N + N*ILAENV( 1, 'ZGEHRD', ' ', N, 1, N, 0 )
00173             MINWRK = 2*N
00174             IF( WANTVL ) THEN
00175                MAXWRK = MAX( MAXWRK, N + ( N - 1 )*ILAENV( 1, 'ZUNGHR',
00176      $                       ' ', N, 1, N, -1 ) )
00177                CALL ZHSEQR( 'S', 'V', N, 1, N, A, LDA, W, VL, LDVL,
00178      $                WORK, -1, INFO )
00179             ELSE IF( WANTVR ) THEN
00180                MAXWRK = MAX( MAXWRK, N + ( N - 1 )*ILAENV( 1, 'ZUNGHR',
00181      $                       ' ', N, 1, N, -1 ) )
00182                CALL ZHSEQR( 'S', 'V', N, 1, N, A, LDA, W, VR, LDVR,
00183      $                WORK, -1, INFO )
00184             ELSE
00185                CALL ZHSEQR( 'E', 'N', N, 1, N, A, LDA, W, VR, LDVR,
00186      $                WORK, -1, INFO )
00187             END IF
00188             HSWORK = WORK( 1 )
00189             MAXWRK = MAX( MAXWRK, HSWORK, MINWRK )
00190          END IF
00191          WORK( 1 ) = MAXWRK
00192 *
00193          IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
00194             INFO = -12
00195          END IF
00196       END IF
00197 *
00198       IF( INFO.NE.0 ) THEN
00199          CALL XERBLA( 'ZGEEV ', -INFO )
00200          RETURN
00201       ELSE IF( LQUERY ) THEN
00202          RETURN
00203       END IF
00204 *
00205 *     Quick return if possible
00206 *
00207       IF( N.EQ.0 )
00208      $   RETURN
00209 *
00210 *     Get machine constants
00211 *
00212       EPS = DLAMCH( 'P' )
00213       SMLNUM = DLAMCH( 'S' )
00214       BIGNUM = ONE / SMLNUM
00215       CALL DLABAD( SMLNUM, BIGNUM )
00216       SMLNUM = SQRT( SMLNUM ) / EPS
00217       BIGNUM = ONE / SMLNUM
00218 *
00219 *     Scale A if max element outside range [SMLNUM,BIGNUM]
00220 *
00221       ANRM = ZLANGE( 'M', N, N, A, LDA, DUM )
00222       SCALEA = .FALSE.
00223       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
00224          SCALEA = .TRUE.
00225          CSCALE = SMLNUM
00226       ELSE IF( ANRM.GT.BIGNUM ) THEN
00227          SCALEA = .TRUE.
00228          CSCALE = BIGNUM
00229       END IF
00230       IF( SCALEA )
00231      $   CALL ZLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR )
00232 *
00233 *     Balance the matrix
00234 *     (CWorkspace: none)
00235 *     (RWorkspace: need N)
00236 *
00237       IBAL = 1
00238       CALL ZGEBAL( 'B', N, A, LDA, ILO, IHI, RWORK( IBAL ), IERR )
00239 *
00240 *     Reduce to upper Hessenberg form
00241 *     (CWorkspace: need 2*N, prefer N+N*NB)
00242 *     (RWorkspace: none)
00243 *
00244       ITAU = 1
00245       IWRK = ITAU + N
00246       CALL ZGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ),
00247      $             LWORK-IWRK+1, IERR )
00248 *
00249       IF( WANTVL ) THEN
00250 *
00251 *        Want left eigenvectors
00252 *        Copy Householder vectors to VL
00253 *
00254          SIDE = 'L'
00255          CALL ZLACPY( 'L', N, N, A, LDA, VL, LDVL )
00256 *
00257 *        Generate unitary matrix in VL
00258 *        (CWorkspace: need 2*N-1, prefer N+(N-1)*NB)
00259 *        (RWorkspace: none)
00260 *
00261          CALL ZUNGHR( N, ILO, IHI, VL, LDVL, WORK( ITAU ), WORK( IWRK ),
00262      $                LWORK-IWRK+1, IERR )
00263 *
00264 *        Perform QR iteration, accumulating Schur vectors in VL
00265 *        (CWorkspace: need 1, prefer HSWORK (see comments) )
00266 *        (RWorkspace: none)
00267 *
00268          IWRK = ITAU
00269          CALL ZHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, W, VL, LDVL,
00270      $                WORK( IWRK ), LWORK-IWRK+1, INFO )
00271 *
00272          IF( WANTVR ) THEN
00273 *
00274 *           Want left and right eigenvectors
00275 *           Copy Schur vectors to VR
00276 *
00277             SIDE = 'B'
00278             CALL ZLACPY( 'F', N, N, VL, LDVL, VR, LDVR )
00279          END IF
00280 *
00281       ELSE IF( WANTVR ) THEN
00282 *
00283 *        Want right eigenvectors
00284 *        Copy Householder vectors to VR
00285 *
00286          SIDE = 'R'
00287          CALL ZLACPY( 'L', N, N, A, LDA, VR, LDVR )
00288 *
00289 *        Generate unitary matrix in VR
00290 *        (CWorkspace: need 2*N-1, prefer N+(N-1)*NB)
00291 *        (RWorkspace: none)
00292 *
00293          CALL ZUNGHR( N, ILO, IHI, VR, LDVR, WORK( ITAU ), WORK( IWRK ),
00294      $                LWORK-IWRK+1, IERR )
00295 *
00296 *        Perform QR iteration, accumulating Schur vectors in VR
00297 *        (CWorkspace: need 1, prefer HSWORK (see comments) )
00298 *        (RWorkspace: none)
00299 *
00300          IWRK = ITAU
00301          CALL ZHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, W, VR, LDVR,
00302      $                WORK( IWRK ), LWORK-IWRK+1, INFO )
00303 *
00304       ELSE
00305 *
00306 *        Compute eigenvalues only
00307 *        (CWorkspace: need 1, prefer HSWORK (see comments) )
00308 *        (RWorkspace: none)
00309 *
00310          IWRK = ITAU
00311          CALL ZHSEQR( 'E', 'N', N, ILO, IHI, A, LDA, W, VR, LDVR,
00312      $                WORK( IWRK ), LWORK-IWRK+1, INFO )
00313       END IF
00314 *
00315 *     If INFO > 0 from ZHSEQR, then quit
00316 *
00317       IF( INFO.GT.0 )
00318      $   GO TO 50
00319 *
00320       IF( WANTVL .OR. WANTVR ) THEN
00321 *
00322 *        Compute left and/or right eigenvectors
00323 *        (CWorkspace: need 2*N)
00324 *        (RWorkspace: need 2*N)
00325 *
00326          IRWORK = IBAL + N
00327          CALL ZTREVC( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
00328      $                N, NOUT, WORK( IWRK ), RWORK( IRWORK ), IERR )
00329       END IF
00330 *
00331       IF( WANTVL ) THEN
00332 *
00333 *        Undo balancing of left eigenvectors
00334 *        (CWorkspace: none)
00335 *        (RWorkspace: need N)
00336 *
00337          CALL ZGEBAK( 'B', 'L', N, ILO, IHI, RWORK( IBAL ), N, VL, LDVL,
00338      $                IERR )
00339 *
00340 *        Normalize left eigenvectors and make largest component real
00341 *
00342          DO 20 I = 1, N
00343             SCL = ONE / DZNRM2( N, VL( 1, I ), 1 )
00344             CALL ZDSCAL( N, SCL, VL( 1, I ), 1 )
00345             DO 10 K = 1, N
00346                RWORK( IRWORK+K-1 ) = DBLE( VL( K, I ) )**2 +
00347      $                               DIMAG( VL( K, I ) )**2
00348    10       CONTINUE
00349             K = IDAMAX( N, RWORK( IRWORK ), 1 )
00350             TMP = DCONJG( VL( K, I ) ) / SQRT( RWORK( IRWORK+K-1 ) )
00351             CALL ZSCAL( N, TMP, VL( 1, I ), 1 )
00352             VL( K, I ) = DCMPLX( DBLE( VL( K, I ) ), ZERO )
00353    20    CONTINUE
00354       END IF
00355 *
00356       IF( WANTVR ) THEN
00357 *
00358 *        Undo balancing of right eigenvectors
00359 *        (CWorkspace: none)
00360 *        (RWorkspace: need N)
00361 *
00362          CALL ZGEBAK( 'B', 'R', N, ILO, IHI, RWORK( IBAL ), N, VR, LDVR,
00363      $                IERR )
00364 *
00365 *        Normalize right eigenvectors and make largest component real
00366 *
00367          DO 40 I = 1, N
00368             SCL = ONE / DZNRM2( N, VR( 1, I ), 1 )
00369             CALL ZDSCAL( N, SCL, VR( 1, I ), 1 )
00370             DO 30 K = 1, N
00371                RWORK( IRWORK+K-1 ) = DBLE( VR( K, I ) )**2 +
00372      $                               DIMAG( VR( K, I ) )**2
00373    30       CONTINUE
00374             K = IDAMAX( N, RWORK( IRWORK ), 1 )
00375             TMP = DCONJG( VR( K, I ) ) / SQRT( RWORK( IRWORK+K-1 ) )
00376             CALL ZSCAL( N, TMP, VR( 1, I ), 1 )
00377             VR( K, I ) = DCMPLX( DBLE( VR( K, I ) ), ZERO )
00378    40    CONTINUE
00379       END IF
00380 *
00381 *     Undo scaling if necessary
00382 *
00383    50 CONTINUE
00384       IF( SCALEA ) THEN
00385          CALL ZLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, W( INFO+1 ),
00386      $                MAX( N-INFO, 1 ), IERR )
00387          IF( INFO.GT.0 ) THEN
00388             CALL ZLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, W, N, IERR )
00389          END IF
00390       END IF
00391 *
00392       WORK( 1 ) = MAXWRK
00393       RETURN
00394 *
00395 *     End of ZGEEV
00396 *
00397       END