1#if 1 2 SUBROUTINE UTIL_DGEEV 3 $ ( JOBVL, JOBVR, N, A, LDA, WR, WI, VL, LDVL, VR, 4 $ LDVR, WORK, LWORK, INFO ) 5c $Id$ 6 implicit none 7 CHARACTER JOBVL, JOBVR 8 INTEGER INFO, LDA, LDVL, LDVR, LWORK, N 9 DOUBLE PRECISION A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ), 10 $ WI( * ), WORK( * ), WR( * ) 11c 12#include "errquit.fh" 13#include "mafdecls.fh" 14 character*1 balanc,sense 15 integer iwork 16 integer ilo,ihi 17 double precision abnrm, rconde, rcondv 18 integer l_scale, k_scale 19c 20 balanc='n' 21 sense='n' 22 if (.not. ma_push_get(mt_dbl, n, 'util_dgeev: scale', 23 $ l_scale, k_scale)) call errquit('util_dgeev: maget failed', 24 $ n, MA_ERR) 25 CALL dgeevx( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, WR, WI, 26 $ VL, LDVL, VR, LDVR, ILO, IHI, 27 D dbl_mb(k_scale), ABNRM, 28 $ RCONDE, RCONDV, WORK, LWORK, IWORK, INFO ) 29 if (.not. ma_pop_stack(l_scale)) call errquit 30 $ ('util_dgeev: ma error popping sc', 0, MA_ERR) 31#else 32 SUBROUTINE UTIL_DGEEV 33 $ ( JOBVL, JOBVR, N, A, LDA, WR, WI, VL, LDVL, VR, 34 $ LDVR, WORK, LWORK, INFO ) 35c $Id$ 36c The subdoutine is identical to the original dgeev but had to be renamed, 37c since otherwise, for IBM SP2, dynamically linked dgeev having different 38c interface syntax is going to be used rather than this. 39* 40* -- LAPACK driver routine (version 2.0) -- 41* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 42* Courant Institute, Argonne National Lab, and Rice University 43* September 30, 1994 44* 45* .. Scalar Arguments .. 46 CHARACTER JOBVL, JOBVR 47 INTEGER INFO, LDA, LDVL, LDVR, LWORK, N 48* .. 49* .. Array Arguments .. 50 DOUBLE PRECISION A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ), 51 $ WI( * ), WORK( * ), WR( * ) 52* .. 53* 54* Purpose 55* ======= 56* 57* DGEEV computes for an N-by-N real nonsymmetric matrix A, the 58* eigenvalues and, optionally, the left and/or right eigenvectors. 59* 60* The right eigenvector v(j) of A satisfies 61* A * v(j) = lambda(j) * v(j) 62* where lambda(j) is its eigenvalue. 63* The left eigenvector u(j) of A satisfies 64* u(j)**H * A = lambda(j) * u(j)**H 65* where u(j)**H denotes the conjugate transpose of u(j). 66* 67* The computed eigenvectors are normalized to have Euclidean norm 68* equal to 1 and largest component real. 69* 70* Arguments 71* ========= 72* 73* JOBVL (input) CHARACTER*1 74* = 'N': left eigenvectors of A are not computed; 75* = 'V': left eigenvectors of A are computed. 76* 77* JOBVR (input) CHARACTER*1 78* = 'N': right eigenvectors of A are not computed; 79* = 'V': right eigenvectors of A are computed. 80* 81* N (input) INTEGER 82* The order of the matrix A. N >= 0. 83* 84* A (input/output) DOUBLE PRECISION array, dimension (LDA,N) 85* On entry, the N-by-N matrix A. 86* On exit, A has been overwritten. 87* 88* LDA (input) INTEGER 89* The leading dimension of the array A. LDA >= max(1,N). 90* 91* WR (output) DOUBLE PRECISION array, dimension (N) 92* WI (output) DOUBLE PRECISION array, dimension (N) 93* WR and WI contain the real and imaginary parts, 94* respectively, of the computed eigenvalues. Complex 95* conjugate pairs of eigenvalues appear consecutively 96* with the eigenvalue having the positive imaginary part 97* first. 98* 99* VL (output) DOUBLE PRECISION array, dimension (LDVL,N) 100* If JOBVL = 'V', the left eigenvectors u(j) are stored one 101* after another in the columns of VL, in the same order 102* as their eigenvalues. 103* If JOBVL = 'N', VL is not referenced. 104* If the j-th eigenvalue is real, then u(j) = VL(:,j), 105* the j-th column of VL. 106* If the j-th and (j+1)-st eigenvalues form a complex 107* conjugate pair, then u(j) = VL(:,j) + i*VL(:,j+1) and 108* u(j+1) = VL(:,j) - i*VL(:,j+1). 109* 110* LDVL (input) INTEGER 111* The leading dimension of the array VL. LDVL >= 1; if 112* JOBVL = 'V', LDVL >= N. 113* 114* VR (output) DOUBLE PRECISION array, dimension (LDVR,N) 115* If JOBVR = 'V', the right eigenvectors v(j) are stored one 116* after another in the columns of VR, in the same order 117* as their eigenvalues. 118* If JOBVR = 'N', VR is not referenced. 119* If the j-th eigenvalue is real, then v(j) = VR(:,j), 120* the j-th column of VR. 121* If the j-th and (j+1)-st eigenvalues form a complex 122* conjugate pair, then v(j) = VR(:,j) + i*VR(:,j+1) and 123* v(j+1) = VR(:,j) - i*VR(:,j+1). 124* 125* LDVR (input) INTEGER 126* The leading dimension of the array VR. LDVR >= 1; if 127* JOBVR = 'V', LDVR >= N. 128* 129* WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK) 130* On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 131* 132* LWORK (input) INTEGER 133* The dimension of the array WORK. LWORK >= max(1,3*N), and 134* if JOBVL = 'V' or JOBVR = 'V', LWORK >= 4*N. For good 135* performance, LWORK must generally be larger. 136* 137* INFO (output) INTEGER 138* = 0: successful exit 139* < 0: if INFO = -i, the i-th argument had an illegal value. 140* > 0: if INFO = i, the QR algorithm failed to compute all the 141* eigenvalues, and no eigenvectors have been computed; 142* elements i+1:N of WR and WI contain eigenvalues which 143* have converged. 144* 145* ===================================================================== 146* 147* .. Parameters .. 148 DOUBLE PRECISION ZERO, ONE 149 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 150* .. 151* .. Local Scalars .. 152 LOGICAL SCALEA, WANTVL, WANTVR 153 CHARACTER SIDE 154 INTEGER HSWORK, I, IBAL, IERR, IHI, ILO, ITAU, IWRK, K, 155 $ MAXB, MAXWRK, MINWRK, NOUT 156 DOUBLE PRECISION ANRM, BIGNUM, CS, CSCALE, EPS, R, SCL, SMLNUM, 157 $ SN 158* .. 159* .. Local Arrays .. 160 LOGICAL SELECT( 1 ) 161 DOUBLE PRECISION DUM( 1 ) 162* .. 163* .. External Subroutines .. 164 EXTERNAL DGEBAK, DGEBAL, DGEHRD, DHSEQR, DLABAD, DLACPY, 165 $ DLARTG, DLASCL, DORGHR, DROT, DSCAL, DTREVC, 166 $ XERBLA 167* .. 168* .. External Functions .. 169 LOGICAL LSAME 170 INTEGER IDAMAX, ILAENV 171 DOUBLE PRECISION DLAMCH, DLANGE, DLAPY2, DNRM2 172 EXTERNAL LSAME, IDAMAX, ILAENV, DLAMCH, DLANGE, DLAPY2, 173 $ DNRM2 174* .. 175* .. Intrinsic Functions .. 176 INTRINSIC MAX, MIN, SQRT 177* .. 178* .. Executable Statements .. 179* 180* Test the input arguments 181* 182 INFO = 0 183 WANTVL = LSAME( JOBVL, 'V' ) 184 WANTVR = LSAME( JOBVR, 'V' ) 185 IF( ( .NOT.WANTVL ) .AND. ( .NOT.LSAME( JOBVL, 'N' ) ) ) THEN 186 INFO = -1 187 ELSE IF( ( .NOT.WANTVR ) .AND. ( .NOT.LSAME( JOBVR, 'N' ) ) ) THEN 188 INFO = -2 189 ELSE IF( N.LT.0 ) THEN 190 INFO = -3 191 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 192 INFO = -5 193 ELSE IF( LDVL.LT.1 .OR. ( WANTVL .AND. LDVL.LT.N ) ) THEN 194 INFO = -9 195 ELSE IF( LDVR.LT.1 .OR. ( WANTVR .AND. LDVR.LT.N ) ) THEN 196 INFO = -11 197 END IF 198* 199* Compute workspace 200* (Note: Comments in the code beginning "Workspace:" describe the 201* minimal amount of workspace needed at that point in the code, 202* as well as the preferred amount for good performance. 203* NB refers to the optimal block size for the immediately 204* following subroutine, as returned by ILAENV. 205* HSWORK refers to the workspace preferred by DHSEQR, as 206* calculated below. HSWORK is computed assuming ILO=1 and IHI=N, 207* the worst case.) 208* 209 MINWRK = 1 210 IF( INFO.EQ.0 .AND. LWORK.GE.1 ) THEN 211 MAXWRK = 2*N + N*ILAENV( 1, 'DGEHRD', ' ', N, 1, N, 0 ) 212 IF( ( .NOT.WANTVL ) .AND. ( .NOT.WANTVR ) ) THEN 213 MINWRK = MAX( 1, 3*N ) 214 MAXB = MAX( ILAENV( 8, 'DHSEQR', 'EN', N, 1, N, -1 ), 2 ) 215 K = MIN( MAXB, N, MAX( 2, ILAENV( 4, 'DHSEQR', 'EN', N, 1, 216 $ N, -1 ) ) ) 217 HSWORK = MAX( K*( K+2 ), 2*N ) 218 MAXWRK = MAX( MAXWRK, N+1, N+HSWORK ) 219 ELSE 220 MINWRK = MAX( 1, 4*N ) 221 MAXWRK = MAX( MAXWRK, 2*N+( N-1 )* 222 $ ILAENV( 1, 'DORGHR', ' ', N, 1, N, -1 ) ) 223 MAXB = MAX( ILAENV( 8, 'DHSEQR', 'SV', N, 1, N, -1 ), 2 ) 224 K = MIN( MAXB, N, MAX( 2, ILAENV( 4, 'DHSEQR', 'SV', N, 1, 225 $ N, -1 ) ) ) 226 HSWORK = MAX( K*( K+2 ), 2*N ) 227 MAXWRK = MAX( MAXWRK, N+1, N+HSWORK ) 228 MAXWRK = MAX( MAXWRK, 4*N ) 229 END IF 230 WORK( 1 ) = MAXWRK 231 END IF 232 IF( LWORK.LT.MINWRK ) THEN 233 INFO = -13 234 END IF 235 IF( INFO.NE.0 ) THEN 236 CALL XERBLA( 'DGEEV ', -INFO ) 237 RETURN 238 END IF 239* 240* Quick return if possible 241* 242 IF( N.EQ.0 ) 243 $ RETURN 244* 245* Get machine constants 246* 247 EPS = DLAMCH( 'P' ) 248 SMLNUM = DLAMCH( 'S' ) 249 BIGNUM = ONE / SMLNUM 250 CALL DLABAD( SMLNUM, BIGNUM ) 251 SMLNUM = SQRT( SMLNUM ) / EPS 252 BIGNUM = ONE / SMLNUM 253* 254* Scale A if max element outside range [SMLNUM,BIGNUM] 255* 256 ANRM = DLANGE( 'M', N, N, A, LDA, DUM ) 257 SCALEA = .FALSE. 258 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 259 SCALEA = .TRUE. 260 CSCALE = SMLNUM 261 ELSE IF( ANRM.GT.BIGNUM ) THEN 262 SCALEA = .TRUE. 263 CSCALE = BIGNUM 264 END IF 265 IF( SCALEA ) 266 $ CALL DLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR ) 267* 268* Balance the matrix 269* (Workspace: need N) 270* 271 IBAL = 1 272 CALL DGEBAL( 'B', N, A, LDA, ILO, IHI, WORK( IBAL ), IERR ) 273* 274* Reduce to upper Hessenberg form 275* (Workspace: need 3*N, prefer 2*N+N*NB) 276* 277 ITAU = IBAL + N 278 IWRK = ITAU + N 279 CALL DGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ), 280 $ LWORK-IWRK+1, IERR ) 281* 282 IF( WANTVL ) THEN 283* 284* Want left eigenvectors 285* Copy Householder vectors to VL 286* 287 SIDE = 'L' 288 CALL DLACPY( 'L', N, N, A, LDA, VL, LDVL ) 289* 290* Generate orthogonal matrix in VL 291* (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB) 292* 293 CALL DORGHR( N, ILO, IHI, VL, LDVL, WORK( ITAU ), WORK( IWRK ), 294 $ LWORK-IWRK+1, IERR ) 295* 296* Perform QR iteration, accumulating Schur vectors in VL 297* (Workspace: need N+1, prefer N+HSWORK (see comments) ) 298* 299 IWRK = ITAU 300 CALL DHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, WR, WI, VL, LDVL, 301 $ WORK( IWRK ), LWORK-IWRK+1, INFO ) 302* 303 IF( WANTVR ) THEN 304* 305* Want left and right eigenvectors 306* Copy Schur vectors to VR 307* 308 SIDE = 'B' 309 CALL DLACPY( 'F', N, N, VL, LDVL, VR, LDVR ) 310 END IF 311* 312 ELSE IF( WANTVR ) THEN 313* 314* Want right eigenvectors 315* Copy Householder vectors to VR 316* 317 SIDE = 'R' 318 CALL DLACPY( 'L', N, N, A, LDA, VR, LDVR ) 319* 320* Generate orthogonal matrix in VR 321* (Workspace: need 3*N-1, prefer 2*N+(N-1)*NB) 322* 323 CALL DORGHR( N, ILO, IHI, VR, LDVR, WORK( ITAU ), WORK( IWRK ), 324 $ LWORK-IWRK+1, IERR ) 325* 326* Perform QR iteration, accumulating Schur vectors in VR 327* (Workspace: need N+1, prefer N+HSWORK (see comments) ) 328* 329 IWRK = ITAU 330 CALL DHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, WR, WI, VR, LDVR, 331 $ WORK( IWRK ), LWORK-IWRK+1, INFO ) 332* 333 ELSE 334* 335* Compute eigenvalues only 336* (Workspace: need N+1, prefer N+HSWORK (see comments) ) 337* 338 IWRK = ITAU 339 CALL DHSEQR( 'E', 'N', N, ILO, IHI, A, LDA, WR, WI, VR, LDVR, 340 $ WORK( IWRK ), LWORK-IWRK+1, INFO ) 341 END IF 342* 343* If INFO > 0 from DHSEQR, then quit 344* 345 IF( INFO.GT.0 ) 346 $ GO TO 50 347* 348 IF( WANTVL .OR. WANTVR ) THEN 349* 350* Compute left and/or right eigenvectors 351* (Workspace: need 4*N) 352* 353 CALL DTREVC( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR, 354 $ N, NOUT, WORK( IWRK ), IERR ) 355 END IF 356* 357 IF( WANTVL ) THEN 358* 359* Undo balancing of left eigenvectors 360* (Workspace: need N) 361* 362 CALL DGEBAK( 'B', 'L', N, ILO, IHI, WORK( IBAL ), N, VL, LDVL, 363 $ IERR ) 364* 365* Normalize left eigenvectors and make largest component real 366* 367 DO 20 I = 1, N 368 IF( WI( I ).EQ.ZERO ) THEN 369 SCL = ONE / DNRM2( N, VL( 1, I ), 1 ) 370 CALL DSCAL( N, SCL, VL( 1, I ), 1 ) 371 ELSE IF( WI( I ).GT.ZERO ) THEN 372 SCL = ONE / DLAPY2( DNRM2( N, VL( 1, I ), 1 ), 373 $ DNRM2( N, VL( 1, I+1 ), 1 ) ) 374 CALL DSCAL( N, SCL, VL( 1, I ), 1 ) 375 CALL DSCAL( N, SCL, VL( 1, I+1 ), 1 ) 376 DO 10 K = 1, N 377 WORK( IWRK+K-1 ) = VL( K, I )**2 + VL( K, I+1 )**2 378 10 CONTINUE 379 K = IDAMAX( N, WORK( IWRK ), 1 ) 380 CALL DLARTG( VL( K, I ), VL( K, I+1 ), CS, SN, R ) 381 CALL DROT( N, VL( 1, I ), 1, VL( 1, I+1 ), 1, CS, SN ) 382 VL( K, I+1 ) = ZERO 383 END IF 384 20 CONTINUE 385 END IF 386* 387 IF( WANTVR ) THEN 388* 389* Undo balancing of right eigenvectors 390* (Workspace: need N) 391* 392 CALL DGEBAK( 'B', 'R', N, ILO, IHI, WORK( IBAL ), N, VR, LDVR, 393 $ IERR ) 394* 395* Normalize right eigenvectors and make largest component real 396* 397 DO 40 I = 1, N 398 IF( WI( I ).EQ.ZERO ) THEN 399 SCL = ONE / DNRM2( N, VR( 1, I ), 1 ) 400 CALL DSCAL( N, SCL, VR( 1, I ), 1 ) 401 ELSE IF( WI( I ).GT.ZERO ) THEN 402 SCL = ONE / DLAPY2( DNRM2( N, VR( 1, I ), 1 ), 403 $ DNRM2( N, VR( 1, I+1 ), 1 ) ) 404 CALL DSCAL( N, SCL, VR( 1, I ), 1 ) 405 CALL DSCAL( N, SCL, VR( 1, I+1 ), 1 ) 406 DO 30 K = 1, N 407 WORK( IWRK+K-1 ) = VR( K, I )**2 + VR( K, I+1 )**2 408 30 CONTINUE 409 K = IDAMAX( N, WORK( IWRK ), 1 ) 410 CALL DLARTG( VR( K, I ), VR( K, I+1 ), CS, SN, R ) 411 CALL DROT( N, VR( 1, I ), 1, VR( 1, I+1 ), 1, CS, SN ) 412 VR( K, I+1 ) = ZERO 413 END IF 414 40 CONTINUE 415 END IF 416* 417* Undo scaling if necessary 418* 419 50 CONTINUE 420 IF( SCALEA ) THEN 421 CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, WR( INFO+1 ), 422 $ MAX( N-INFO, 1 ), IERR ) 423 CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, WI( INFO+1 ), 424 $ MAX( N-INFO, 1 ), IERR ) 425 IF( INFO.GT.0 ) THEN 426 CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WR, N, 427 $ IERR ) 428 CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WI, N, 429 $ IERR ) 430 END IF 431 END IF 432* 433 WORK( 1 ) = MAXWRK 434 RETURN 435* 436* End of DGEEV 437* 438#endif 439 END 440