1*> \brief <b> SGEEV 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 SGEEV + dependencies 10*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgeev.f"> 11*> [TGZ]</a> 12*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgeev.f"> 13*> [ZIP]</a> 14*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgeev.f"> 15*> [TXT]</a> 16*> \endhtmlonly 17* 18* Definition: 19* =========== 20* 21* SUBROUTINE SGEEV( 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* REAL A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ), 30* $ WI( * ), WORK( * ), WR( * ) 31* .. 32* 33* 34*> \par Purpose: 35* ============= 36*> 37*> \verbatim 38*> 39*> SGEEV 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 REAL 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 REAL array, dimension (N) 92*> \endverbatim 93*> 94*> \param[out] WI 95*> \verbatim 96*> WI is REAL 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 REAL 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 REAL 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 REAL 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* 185* @generated from dgeev.f, fortran d -> s, Tue Apr 19 01:47:44 2016 186* 187*> \ingroup realGEeigen 188* 189* ===================================================================== 190 SUBROUTINE SGEEV( JOBVL, JOBVR, N, A, LDA, WR, WI, VL, LDVL, VR, 191 $ LDVR, WORK, LWORK, INFO ) 192 implicit none 193* 194* -- LAPACK driver routine -- 195* -- LAPACK is a software package provided by Univ. of Tennessee, -- 196* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 197* 198* .. Scalar Arguments .. 199 CHARACTER JOBVL, JOBVR 200 INTEGER INFO, LDA, LDVL, LDVR, LWORK, N 201* .. 202* .. Array Arguments .. 203 REAL A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ), 204 $ WI( * ), WORK( * ), WR( * ) 205* .. 206* 207* ===================================================================== 208* 209* .. Parameters .. 210 REAL ZERO, ONE 211 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 ) 212* .. 213* .. Local Scalars .. 214 LOGICAL LQUERY, SCALEA, WANTVL, WANTVR 215 CHARACTER SIDE 216 INTEGER HSWORK, I, IBAL, IERR, IHI, ILO, ITAU, IWRK, K, 217 $ LWORK_TREVC, MAXWRK, MINWRK, NOUT 218 REAL ANRM, BIGNUM, CS, CSCALE, EPS, R, SCL, SMLNUM, 219 $ SN 220* .. 221* .. Local Arrays .. 222 LOGICAL SELECT( 1 ) 223 REAL DUM( 1 ) 224* .. 225* .. External Subroutines .. 226 EXTERNAL SGEBAK, SGEBAL, SGEHRD, SHSEQR, SLABAD, SLACPY, 227 $ SLARTG, SLASCL, SORGHR, SROT, SSCAL, STREVC3, 228 $ XERBLA 229* .. 230* .. External Functions .. 231 LOGICAL LSAME 232 INTEGER ISAMAX, ILAENV 233 REAL SLAMCH, SLANGE, SLAPY2, SNRM2 234 EXTERNAL LSAME, ISAMAX, ILAENV, SLAMCH, SLANGE, SLAPY2, 235 $ SNRM2 236* .. 237* .. Intrinsic Functions .. 238 INTRINSIC MAX, SQRT 239* .. 240* .. Executable Statements .. 241* 242* Test the input arguments 243* 244 INFO = 0 245 LQUERY = ( LWORK.EQ.-1 ) 246 WANTVL = LSAME( JOBVL, 'V' ) 247 WANTVR = LSAME( JOBVR, 'V' ) 248 IF( ( .NOT.WANTVL ) .AND. ( .NOT.LSAME( JOBVL, 'N' ) ) ) THEN 249 INFO = -1 250 ELSE IF( ( .NOT.WANTVR ) .AND. ( .NOT.LSAME( JOBVR, 'N' ) ) ) THEN 251 INFO = -2 252 ELSE IF( N.LT.0 ) THEN 253 INFO = -3 254 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 255 INFO = -5 256 ELSE IF( LDVL.LT.1 .OR. ( WANTVL .AND. LDVL.LT.N ) ) THEN 257 INFO = -9 258 ELSE IF( LDVR.LT.1 .OR. ( WANTVR .AND. LDVR.LT.N ) ) THEN 259 INFO = -11 260 END IF 261* 262* Compute workspace 263* (Note: Comments in the code beginning "Workspace:" describe the 264* minimal amount of workspace needed at that point in the code, 265* as well as the preferred amount for good performance. 266* NB refers to the optimal block size for the immediately 267* following subroutine, as returned by ILAENV. 268* HSWORK refers to the workspace preferred by SHSEQR, as 269* calculated below. HSWORK is computed assuming ILO=1 and IHI=N, 270* the worst case.) 271* 272 IF( INFO.EQ.0 ) THEN 273 IF( N.EQ.0 ) THEN 274 MINWRK = 1 275 MAXWRK = 1 276 ELSE 277 MAXWRK = 2*N + N*ILAENV( 1, 'SGEHRD', ' ', N, 1, N, 0 ) 278 IF( WANTVL ) THEN 279 MINWRK = 4*N 280 MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1, 281 $ 'SORGHR', ' ', N, 1, N, -1 ) ) 282 CALL SHSEQR( 'S', 'V', N, 1, N, A, LDA, WR, WI, VL, LDVL, 283 $ WORK, -1, INFO ) 284 HSWORK = INT( WORK(1) ) 285 MAXWRK = MAX( MAXWRK, N + 1, N + HSWORK ) 286 CALL STREVC3( 'L', 'B', SELECT, N, A, LDA, 287 $ VL, LDVL, VR, LDVR, N, NOUT, 288 $ WORK, -1, IERR ) 289 LWORK_TREVC = INT( WORK(1) ) 290 MAXWRK = MAX( MAXWRK, N + LWORK_TREVC ) 291 MAXWRK = MAX( MAXWRK, 4*N ) 292 ELSE IF( WANTVR ) THEN 293 MINWRK = 4*N 294 MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1, 295 $ 'SORGHR', ' ', N, 1, N, -1 ) ) 296 CALL SHSEQR( 'S', 'V', N, 1, N, A, LDA, WR, WI, VR, LDVR, 297 $ WORK, -1, INFO ) 298 HSWORK = INT( WORK(1) ) 299 MAXWRK = MAX( MAXWRK, N + 1, N + HSWORK ) 300 CALL STREVC3( 'R', 'B', SELECT, N, A, LDA, 301 $ VL, LDVL, VR, LDVR, N, NOUT, 302 $ WORK, -1, IERR ) 303 LWORK_TREVC = INT( WORK(1) ) 304 MAXWRK = MAX( MAXWRK, N + LWORK_TREVC ) 305 MAXWRK = MAX( MAXWRK, 4*N ) 306 ELSE 307 MINWRK = 3*N 308 CALL SHSEQR( 'E', 'N', N, 1, N, A, LDA, WR, WI, VR, LDVR, 309 $ WORK, -1, INFO ) 310 HSWORK = INT( WORK(1) ) 311 MAXWRK = MAX( MAXWRK, N + 1, N + HSWORK ) 312 END IF 313 MAXWRK = MAX( MAXWRK, MINWRK ) 314 END IF 315 WORK( 1 ) = MAXWRK 316* 317 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN 318 INFO = -13 319 END IF 320 END IF 321* 322 IF( INFO.NE.0 ) THEN 323 CALL XERBLA( 'SGEEV ', -INFO ) 324 RETURN 325 ELSE IF( LQUERY ) THEN 326 RETURN 327 END IF 328* 329* Quick return if possible 330* 331 IF( N.EQ.0 ) 332 $ RETURN 333* 334* Get machine constants 335* 336 EPS = SLAMCH( 'P' ) 337 SMLNUM = SLAMCH( 'S' ) 338 BIGNUM = ONE / SMLNUM 339 CALL SLABAD( SMLNUM, BIGNUM ) 340 SMLNUM = SQRT( SMLNUM ) / EPS 341 BIGNUM = ONE / SMLNUM 342* 343* Scale A if max element outside range [SMLNUM,BIGNUM] 344* 345 ANRM = SLANGE( 'M', N, N, A, LDA, DUM ) 346 SCALEA = .FALSE. 347 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 348 SCALEA = .TRUE. 349 CSCALE = SMLNUM 350 ELSE IF( ANRM.GT.BIGNUM ) THEN 351 SCALEA = .TRUE. 352 CSCALE = BIGNUM 353 END IF 354 IF( SCALEA ) 355 $ CALL SLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR ) 356* 357* Balance the matrix 358* (Workspace: need N) 359* 360 IBAL = 1 361 CALL SGEBAL( 'B', N, A, LDA, ILO, IHI, WORK( IBAL ), IERR ) 362* 363* Reduce to upper Hessenberg form 364* (Workspace: need 3*N, prefer 2*N+N*NB) 365* 366 ITAU = IBAL + N 367 IWRK = ITAU + N 368 CALL SGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ), 369 $ LWORK-IWRK+1, IERR ) 370* 371 IF( WANTVL ) THEN 372* 373* Want left eigenvectors 374* Copy Householder vectors to VL 375* 376 SIDE = 'L' 377 CALL SLACPY( 'L', N, N, A, LDA, VL, LDVL ) 378* 379* Generate orthogonal matrix in VL 380* (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB) 381* 382 CALL SORGHR( N, ILO, IHI, VL, LDVL, WORK( ITAU ), WORK( IWRK ), 383 $ LWORK-IWRK+1, IERR ) 384* 385* Perform QR iteration, accumulating Schur vectors in VL 386* (Workspace: need N+1, prefer N+HSWORK (see comments) ) 387* 388 IWRK = ITAU 389 CALL SHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, WR, WI, VL, LDVL, 390 $ WORK( IWRK ), LWORK-IWRK+1, INFO ) 391* 392 IF( WANTVR ) THEN 393* 394* Want left and right eigenvectors 395* Copy Schur vectors to VR 396* 397 SIDE = 'B' 398 CALL SLACPY( 'F', N, N, VL, LDVL, VR, LDVR ) 399 END IF 400* 401 ELSE IF( WANTVR ) THEN 402* 403* Want right eigenvectors 404* Copy Householder vectors to VR 405* 406 SIDE = 'R' 407 CALL SLACPY( 'L', N, N, A, LDA, VR, LDVR ) 408* 409* Generate orthogonal matrix in VR 410* (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB) 411* 412 CALL SORGHR( N, ILO, IHI, VR, LDVR, WORK( ITAU ), WORK( IWRK ), 413 $ LWORK-IWRK+1, IERR ) 414* 415* Perform QR iteration, accumulating Schur vectors in VR 416* (Workspace: need N+1, prefer N+HSWORK (see comments) ) 417* 418 IWRK = ITAU 419 CALL SHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, WR, WI, VR, LDVR, 420 $ WORK( IWRK ), LWORK-IWRK+1, INFO ) 421* 422 ELSE 423* 424* Compute eigenvalues only 425* (Workspace: need N+1, prefer N+HSWORK (see comments) ) 426* 427 IWRK = ITAU 428 CALL SHSEQR( 'E', 'N', N, ILO, IHI, A, LDA, WR, WI, VR, LDVR, 429 $ WORK( IWRK ), LWORK-IWRK+1, INFO ) 430 END IF 431* 432* If INFO .NE. 0 from SHSEQR, then quit 433* 434 IF( INFO.NE.0 ) 435 $ GO TO 50 436* 437 IF( WANTVL .OR. WANTVR ) THEN 438* 439* Compute left and/or right eigenvectors 440* (Workspace: need 4*N, prefer N + N + 2*N*NB) 441* 442 CALL STREVC3( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR, 443 $ N, NOUT, WORK( IWRK ), LWORK-IWRK+1, IERR ) 444 END IF 445* 446 IF( WANTVL ) THEN 447* 448* Undo balancing of left eigenvectors 449* (Workspace: need N) 450* 451 CALL SGEBAK( 'B', 'L', N, ILO, IHI, WORK( IBAL ), N, VL, LDVL, 452 $ IERR ) 453* 454* Normalize left eigenvectors and make largest component real 455* 456 DO 20 I = 1, N 457 IF( WI( I ).EQ.ZERO ) THEN 458 SCL = ONE / SNRM2( N, VL( 1, I ), 1 ) 459 CALL SSCAL( N, SCL, VL( 1, I ), 1 ) 460 ELSE IF( WI( I ).GT.ZERO ) THEN 461 SCL = ONE / SLAPY2( SNRM2( N, VL( 1, I ), 1 ), 462 $ SNRM2( N, VL( 1, I+1 ), 1 ) ) 463 CALL SSCAL( N, SCL, VL( 1, I ), 1 ) 464 CALL SSCAL( N, SCL, VL( 1, I+1 ), 1 ) 465 DO 10 K = 1, N 466 WORK( IWRK+K-1 ) = VL( K, I )**2 + VL( K, I+1 )**2 467 10 CONTINUE 468 K = ISAMAX( N, WORK( IWRK ), 1 ) 469 CALL SLARTG( VL( K, I ), VL( K, I+1 ), CS, SN, R ) 470 CALL SROT( N, VL( 1, I ), 1, VL( 1, I+1 ), 1, CS, SN ) 471 VL( K, I+1 ) = ZERO 472 END IF 473 20 CONTINUE 474 END IF 475* 476 IF( WANTVR ) THEN 477* 478* Undo balancing of right eigenvectors 479* (Workspace: need N) 480* 481 CALL SGEBAK( 'B', 'R', N, ILO, IHI, WORK( IBAL ), N, VR, LDVR, 482 $ IERR ) 483* 484* Normalize right eigenvectors and make largest component real 485* 486 DO 40 I = 1, N 487 IF( WI( I ).EQ.ZERO ) THEN 488 SCL = ONE / SNRM2( N, VR( 1, I ), 1 ) 489 CALL SSCAL( N, SCL, VR( 1, I ), 1 ) 490 ELSE IF( WI( I ).GT.ZERO ) THEN 491 SCL = ONE / SLAPY2( SNRM2( N, VR( 1, I ), 1 ), 492 $ SNRM2( N, VR( 1, I+1 ), 1 ) ) 493 CALL SSCAL( N, SCL, VR( 1, I ), 1 ) 494 CALL SSCAL( N, SCL, VR( 1, I+1 ), 1 ) 495 DO 30 K = 1, N 496 WORK( IWRK+K-1 ) = VR( K, I )**2 + VR( K, I+1 )**2 497 30 CONTINUE 498 K = ISAMAX( N, WORK( IWRK ), 1 ) 499 CALL SLARTG( VR( K, I ), VR( K, I+1 ), CS, SN, R ) 500 CALL SROT( N, VR( 1, I ), 1, VR( 1, I+1 ), 1, CS, SN ) 501 VR( K, I+1 ) = ZERO 502 END IF 503 40 CONTINUE 504 END IF 505* 506* Undo scaling if necessary 507* 508 50 CONTINUE 509 IF( SCALEA ) THEN 510 CALL SLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, WR( INFO+1 ), 511 $ MAX( N-INFO, 1 ), IERR ) 512 CALL SLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, WI( INFO+1 ), 513 $ MAX( N-INFO, 1 ), IERR ) 514 IF( INFO.GT.0 ) THEN 515 CALL SLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WR, N, 516 $ IERR ) 517 CALL SLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WI, N, 518 $ IERR ) 519 END IF 520 END IF 521* 522 WORK( 1 ) = MAXWRK 523 RETURN 524* 525* End of SGEEV 526* 527 END 528