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