1 SUBROUTINE ZGEEV( JOBVL, JOBVR, N, A, LDA, W, VL, LDVL, VR, LDVR, 2 $ WORK, LWORK, RWORK, INFO ) 3* 4* -- LAPACK driver routine (version 3.0) -- 5* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 6* Courant Institute, Argonne National Lab, and Rice University 7* June 30, 1999 8* 9* .. Scalar Arguments .. 10 CHARACTER JOBVL, JOBVR 11 INTEGER INFO, LDA, LDVL, LDVR, LWORK, N 12* .. 13* .. Array Arguments .. 14 DOUBLE PRECISION RWORK( * ) 15 COMPLEX*16 A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ), 16 $ W( * ), WORK( * ) 17* .. 18* 19* Purpose 20* ======= 21* 22* ZGEEV computes for an N-by-N complex nonsymmetric matrix A, the 23* eigenvalues and, optionally, the left and/or right eigenvectors. 24* 25* The right eigenvector v(j) of A satisfies 26* A * v(j) = lambda(j) * v(j) 27* where lambda(j) is its eigenvalue. 28* The left eigenvector u(j) of A satisfies 29* u(j)**H * A = lambda(j) * u(j)**H 30* where u(j)**H denotes the conjugate transpose of u(j). 31* 32* The computed eigenvectors are normalized to have Euclidean norm 33* equal to 1 and largest component real. 34* 35* Arguments 36* ========= 37* 38* JOBVL (input) CHARACTER*1 39* = 'N': left eigenvectors of A are not computed; 40* = 'V': left eigenvectors of are computed. 41* 42* JOBVR (input) CHARACTER*1 43* = 'N': right eigenvectors of A are not computed; 44* = 'V': right eigenvectors of A are computed. 45* 46* N (input) INTEGER 47* The order of the matrix A. N >= 0. 48* 49* A (input/output) COMPLEX*16 array, dimension (LDA,N) 50* On entry, the N-by-N matrix A. 51* On exit, A has been overwritten. 52* 53* LDA (input) INTEGER 54* The leading dimension of the array A. LDA >= max(1,N). 55* 56* W (output) COMPLEX*16 array, dimension (N) 57* W contains the computed eigenvalues. 58* 59* VL (output) COMPLEX*16 array, dimension (LDVL,N) 60* If JOBVL = 'V', the left eigenvectors u(j) are stored one 61* after another in the columns of VL, in the same order 62* as their eigenvalues. 63* If JOBVL = 'N', VL is not referenced. 64* u(j) = VL(:,j), the j-th column of VL. 65* 66* LDVL (input) INTEGER 67* The leading dimension of the array VL. LDVL >= 1; if 68* JOBVL = 'V', LDVL >= N. 69* 70* VR (output) COMPLEX*16 array, dimension (LDVR,N) 71* If JOBVR = 'V', the right eigenvectors v(j) are stored one 72* after another in the columns of VR, in the same order 73* as their eigenvalues. 74* If JOBVR = 'N', VR is not referenced. 75* v(j) = VR(:,j), the j-th column of VR. 76* 77* LDVR (input) INTEGER 78* The leading dimension of the array VR. LDVR >= 1; if 79* JOBVR = 'V', LDVR >= N. 80* 81* WORK (workspace/output) COMPLEX*16 array, dimension (LWORK) 82* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 83* 84* LWORK (input) INTEGER 85* The dimension of the array WORK. LWORK >= max(1,2*N). 86* For good performance, LWORK must generally be larger. 87* 88* If LWORK = -1, then a workspace query is assumed; the routine 89* only calculates the optimal size of the WORK array, returns 90* this value as the first entry of the WORK array, and no error 91* message related to LWORK is issued by XERBLA. 92* 93* RWORK (workspace) DOUBLE PRECISION array, dimension (2*N) 94* 95* INFO (output) INTEGER 96* = 0: successful exit 97* < 0: if INFO = -i, the i-th argument had an illegal value. 98* > 0: if INFO = i, the QR algorithm failed to compute all the 99* eigenvalues, and no eigenvectors have been computed; 100* elements and i+1:N of W contain eigenvalues which have 101* converged. 102* 103* ===================================================================== 104* 105* .. Parameters .. 106 DOUBLE PRECISION ZERO, ONE 107 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 108* .. 109* .. Local Scalars .. 110 LOGICAL LQUERY, SCALEA, WANTVL, WANTVR 111 CHARACTER SIDE 112 INTEGER HSWORK, I, IBAL, IERR, IHI, ILO, IRWORK, ITAU, 113 $ IWRK, K, MAXB, MAXWRK, MINWRK, NOUT 114 DOUBLE PRECISION ANRM, BIGNUM, CSCALE, EPS, SCL, SMLNUM 115 COMPLEX*16 TMP 116* .. 117* .. Local Arrays .. 118 LOGICAL SELECT( 1 ) 119 DOUBLE PRECISION DUM( 1 ) 120* .. 121* .. External Subroutines .. 122 EXTERNAL XERBLA, ZDSCAL, ZGEBAK, ZGEBAL, ZGEHRD, ZHSEQR, 123 $ ZLACPY, ZLASCL, ZSCAL, ZTREVC, ZUNGHR 124* .. 125* .. External Functions .. 126 LOGICAL LSAME 127 INTEGER IDAMAX, ILAENV 128 DOUBLE PRECISION DLAMCH, DZNRM2, ZLANGE 129 EXTERNAL LSAME, IDAMAX, ILAENV, DLAMCH, DZNRM2, ZLANGE 130* .. 131* .. Intrinsic Functions .. 132 INTRINSIC DBLE, DCMPLX, DCONJG, DIMAG, MAX, MIN, SQRT 133* .. 134* .. Executable Statements .. 135* 136* Test the input arguments 137* 138 INFO = 0 139 LQUERY = ( LWORK.EQ.-1 ) 140 WANTVL = LSAME( JOBVL, 'V' ) 141 WANTVR = LSAME( JOBVR, 'V' ) 142 IF( ( .NOT.WANTVL ) .AND. ( .NOT.LSAME( JOBVL, 'N' ) ) ) THEN 143 INFO = -1 144 ELSE IF( ( .NOT.WANTVR ) .AND. ( .NOT.LSAME( JOBVR, 'N' ) ) ) THEN 145 INFO = -2 146 ELSE IF( N.LT.0 ) THEN 147 INFO = -3 148 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 149 INFO = -5 150 ELSE IF( LDVL.LT.1 .OR. ( WANTVL .AND. LDVL.LT.N ) ) THEN 151 INFO = -8 152 ELSE IF( LDVR.LT.1 .OR. ( WANTVR .AND. LDVR.LT.N ) ) THEN 153 INFO = -10 154 END IF 155* 156* Compute workspace 157* (Note: Comments in the code beginning "Workspace:" describe the 158* minimal amount of workspace needed at that point in the code, 159* as well as the preferred amount for good performance. 160* CWorkspace refers to complex workspace, and RWorkspace to real 161* workspace. NB refers to the optimal block size for the 162* immediately following subroutine, as returned by ILAENV. 163* HSWORK refers to the workspace preferred by ZHSEQR, as 164* calculated below. HSWORK is computed assuming ILO=1 and IHI=N, 165* the worst case.) 166* 167 MINWRK = 1 168 IF( INFO.EQ.0 .AND. ( LWORK.GE.1 .OR. LQUERY ) ) THEN 169 MAXWRK = N + N*ILAENV( 1, 'ZGEHRD', ' ', N, 1, N, 0 ) 170 IF( ( .NOT.WANTVL ) .AND. ( .NOT.WANTVR ) ) THEN 171 MINWRK = MAX( 1, 2*N ) 172 MAXB = MAX( ILAENV( 8, 'ZHSEQR', 'EN', N, 1, N, -1 ), 2 ) 173 K = MIN( MAXB, N, MAX( 2, ILAENV( 4, 'ZHSEQR', 'EN', N, 1, 174 $ N, -1 ) ) ) 175 HSWORK = MAX( K*( K+2 ), 2*N ) 176 MAXWRK = MAX( MAXWRK, HSWORK ) 177 ELSE 178 MINWRK = MAX( 1, 2*N ) 179 MAXWRK = MAX( MAXWRK, N+( N-1 )* 180 $ ILAENV( 1, 'ZUNGHR', ' ', N, 1, N, -1 ) ) 181 MAXB = MAX( ILAENV( 8, 'ZHSEQR', 'SV', N, 1, N, -1 ), 2 ) 182 K = MIN( MAXB, N, MAX( 2, ILAENV( 4, 'ZHSEQR', 'SV', N, 1, 183 $ N, -1 ) ) ) 184 HSWORK = MAX( K*( K+2 ), 2*N ) 185 MAXWRK = MAX( MAXWRK, HSWORK, 2*N ) 186 END IF 187 WORK( 1 ) = MAXWRK 188 END IF 189 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN 190 INFO = -12 191 END IF 192 IF( INFO.NE.0 ) THEN 193 CALL XERBLA( 'ZGEEV ', -INFO ) 194 RETURN 195 ELSE IF( LQUERY ) THEN 196 RETURN 197 END IF 198* 199* Quick return if possible 200* 201 IF( N.EQ.0 ) 202 $ RETURN 203* 204* Get machine constants 205* 206 EPS = DLAMCH( 'P' ) 207 SMLNUM = DLAMCH( 'S' ) 208 BIGNUM = ONE / SMLNUM 209 CALL DLABAD( SMLNUM, BIGNUM ) 210 SMLNUM = SQRT( SMLNUM ) / EPS 211 BIGNUM = ONE / SMLNUM 212* 213* Scale A if max element outside range [SMLNUM,BIGNUM] 214* 215 ANRM = ZLANGE( 'M', N, N, A, LDA, DUM ) 216 SCALEA = .FALSE. 217 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 218 SCALEA = .TRUE. 219 CSCALE = SMLNUM 220 ELSE IF( ANRM.GT.BIGNUM ) THEN 221 SCALEA = .TRUE. 222 CSCALE = BIGNUM 223 END IF 224 IF( SCALEA ) 225 $ CALL ZLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR ) 226* 227* Balance the matrix 228* (CWorkspace: none) 229* (RWorkspace: need N) 230* 231 IBAL = 1 232 CALL ZGEBAL( 'B', N, A, LDA, ILO, IHI, RWORK( IBAL ), IERR ) 233* 234* Reduce to upper Hessenberg form 235* (CWorkspace: need 2*N, prefer N+N*NB) 236* (RWorkspace: none) 237* 238 ITAU = 1 239 IWRK = ITAU + N 240 CALL ZGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ), 241 $ LWORK-IWRK+1, IERR ) 242* 243 IF( WANTVL ) THEN 244* 245* Want left eigenvectors 246* Copy Householder vectors to VL 247* 248 SIDE = 'L' 249 CALL ZLACPY( 'L', N, N, A, LDA, VL, LDVL ) 250* 251* Generate unitary matrix in VL 252* (CWorkspace: need 2*N-1, prefer N+(N-1)*NB) 253* (RWorkspace: none) 254* 255 CALL ZUNGHR( N, ILO, IHI, VL, LDVL, WORK( ITAU ), WORK( IWRK ), 256 $ LWORK-IWRK+1, IERR ) 257* 258* Perform QR iteration, accumulating Schur vectors in VL 259* (CWorkspace: need 1, prefer HSWORK (see comments) ) 260* (RWorkspace: none) 261* 262 IWRK = ITAU 263 CALL ZHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, W, VL, LDVL, 264 $ WORK( IWRK ), LWORK-IWRK+1, INFO ) 265* 266 IF( WANTVR ) THEN 267* 268* Want left and right eigenvectors 269* Copy Schur vectors to VR 270* 271 SIDE = 'B' 272 CALL ZLACPY( 'F', N, N, VL, LDVL, VR, LDVR ) 273 END IF 274* 275 ELSE IF( WANTVR ) THEN 276* 277* Want right eigenvectors 278* Copy Householder vectors to VR 279* 280 SIDE = 'R' 281 CALL ZLACPY( 'L', N, N, A, LDA, VR, LDVR ) 282* 283* Generate unitary matrix in VR 284* (CWorkspace: need 2*N-1, prefer N+(N-1)*NB) 285* (RWorkspace: none) 286* 287 CALL ZUNGHR( N, ILO, IHI, VR, LDVR, WORK( ITAU ), WORK( IWRK ), 288 $ LWORK-IWRK+1, IERR ) 289* 290* Perform QR iteration, accumulating Schur vectors in VR 291* (CWorkspace: need 1, prefer HSWORK (see comments) ) 292* (RWorkspace: none) 293* 294 IWRK = ITAU 295 CALL ZHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, W, VR, LDVR, 296 $ WORK( IWRK ), LWORK-IWRK+1, INFO ) 297* 298 ELSE 299* 300* Compute eigenvalues only 301* (CWorkspace: need 1, prefer HSWORK (see comments) ) 302* (RWorkspace: none) 303* 304 IWRK = ITAU 305 CALL ZHSEQR( 'E', 'N', N, ILO, IHI, A, LDA, W, VR, LDVR, 306 $ WORK( IWRK ), LWORK-IWRK+1, INFO ) 307 END IF 308* 309* If INFO > 0 from ZHSEQR, then quit 310* 311 IF( INFO.GT.0 ) 312 $ GO TO 50 313* 314 IF( WANTVL .OR. WANTVR ) THEN 315* 316* Compute left and/or right eigenvectors 317* (CWorkspace: need 2*N) 318* (RWorkspace: need 2*N) 319* 320 IRWORK = IBAL + N 321 CALL ZTREVC( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR, 322 $ N, NOUT, WORK( IWRK ), RWORK( IRWORK ), IERR ) 323 END IF 324* 325 IF( WANTVL ) THEN 326* 327* Undo balancing of left eigenvectors 328* (CWorkspace: none) 329* (RWorkspace: need N) 330* 331 CALL ZGEBAK( 'B', 'L', N, ILO, IHI, RWORK( IBAL ), N, VL, LDVL, 332 $ IERR ) 333* 334* Normalize left eigenvectors and make largest component real 335* 336 DO 20 I = 1, N 337 SCL = ONE / DZNRM2( N, VL( 1, I ), 1 ) 338 CALL ZDSCAL( N, SCL, VL( 1, I ), 1 ) 339 DO 10 K = 1, N 340 RWORK( IRWORK+K-1 ) = DBLE( VL( K, I ) )**2 + 341 $ DIMAG( VL( K, I ) )**2 342 10 CONTINUE 343 K = IDAMAX( N, RWORK( IRWORK ), 1 ) 344 TMP = DCONJG( VL( K, I ) ) / SQRT( RWORK( IRWORK+K-1 ) ) 345 CALL ZSCAL( N, TMP, VL( 1, I ), 1 ) 346 VL( K, I ) = DCMPLX( DBLE( VL( K, I ) ), ZERO ) 347 20 CONTINUE 348 END IF 349* 350 IF( WANTVR ) THEN 351* 352* Undo balancing of right eigenvectors 353* (CWorkspace: none) 354* (RWorkspace: need N) 355* 356 CALL ZGEBAK( 'B', 'R', N, ILO, IHI, RWORK( IBAL ), N, VR, LDVR, 357 $ IERR ) 358* 359* Normalize right eigenvectors and make largest component real 360* 361 DO 40 I = 1, N 362 SCL = ONE / DZNRM2( N, VR( 1, I ), 1 ) 363 CALL ZDSCAL( N, SCL, VR( 1, I ), 1 ) 364 DO 30 K = 1, N 365 RWORK( IRWORK+K-1 ) = DBLE( VR( K, I ) )**2 + 366 $ DIMAG( VR( K, I ) )**2 367 30 CONTINUE 368 K = IDAMAX( N, RWORK( IRWORK ), 1 ) 369 TMP = DCONJG( VR( K, I ) ) / SQRT( RWORK( IRWORK+K-1 ) ) 370 CALL ZSCAL( N, TMP, VR( 1, I ), 1 ) 371 VR( K, I ) = DCMPLX( DBLE( VR( K, I ) ), ZERO ) 372 40 CONTINUE 373 END IF 374* 375* Undo scaling if necessary 376* 377 50 CONTINUE 378 IF( SCALEA ) THEN 379 CALL ZLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, W( INFO+1 ), 380 $ MAX( N-INFO, 1 ), IERR ) 381 IF( INFO.GT.0 ) THEN 382 CALL ZLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, W, N, IERR ) 383 END IF 384 END IF 385* 386 WORK( 1 ) = MAXWRK 387 RETURN 388* 389* End of ZGEEV 390* 391 END 392