![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief \b DTGEVC 00002 * 00003 * =========== DOCUMENTATION =========== 00004 * 00005 * Online html documentation available at 00006 * http://www.netlib.org/lapack/explore-html/ 00007 * 00008 *> \htmlonly 00009 *> Download DTGEVC + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtgevc.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtgevc.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtgevc.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE DTGEVC( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL, 00022 * LDVL, VR, LDVR, MM, M, WORK, INFO ) 00023 * 00024 * .. Scalar Arguments .. 00025 * CHARACTER HOWMNY, SIDE 00026 * INTEGER INFO, LDP, LDS, LDVL, LDVR, M, MM, N 00027 * .. 00028 * .. Array Arguments .. 00029 * LOGICAL SELECT( * ) 00030 * DOUBLE PRECISION P( LDP, * ), S( LDS, * ), VL( LDVL, * ), 00031 * $ VR( LDVR, * ), WORK( * ) 00032 * .. 00033 * 00034 * 00035 * 00036 *> \par Purpose: 00037 * ============= 00038 *> 00039 *> \verbatim 00040 *> 00041 *> DTGEVC computes some or all of the right and/or left eigenvectors of 00042 *> a pair of real matrices (S,P), where S is a quasi-triangular matrix 00043 *> and P is upper triangular. Matrix pairs of this type are produced by 00044 *> the generalized Schur factorization of a matrix pair (A,B): 00045 *> 00046 *> A = Q*S*Z**T, B = Q*P*Z**T 00047 *> 00048 *> as computed by DGGHRD + DHGEQZ. 00049 *> 00050 *> The right eigenvector x and the left eigenvector y of (S,P) 00051 *> corresponding to an eigenvalue w are defined by: 00052 *> 00053 *> S*x = w*P*x, (y**H)*S = w*(y**H)*P, 00054 *> 00055 *> where y**H denotes the conjugate tranpose of y. 00056 *> The eigenvalues are not input to this routine, but are computed 00057 *> directly from the diagonal blocks of S and P. 00058 *> 00059 *> This routine returns the matrices X and/or Y of right and left 00060 *> eigenvectors of (S,P), or the products Z*X and/or Q*Y, 00061 *> where Z and Q are input matrices. 00062 *> If Q and Z are the orthogonal factors from the generalized Schur 00063 *> factorization of a matrix pair (A,B), then Z*X and Q*Y 00064 *> are the matrices of right and left eigenvectors of (A,B). 00065 *> 00066 *> \endverbatim 00067 * 00068 * Arguments: 00069 * ========== 00070 * 00071 *> \param[in] SIDE 00072 *> \verbatim 00073 *> SIDE is CHARACTER*1 00074 *> = 'R': compute right eigenvectors only; 00075 *> = 'L': compute left eigenvectors only; 00076 *> = 'B': compute both right and left eigenvectors. 00077 *> \endverbatim 00078 *> 00079 *> \param[in] HOWMNY 00080 *> \verbatim 00081 *> HOWMNY is CHARACTER*1 00082 *> = 'A': compute all right and/or left eigenvectors; 00083 *> = 'B': compute all right and/or left eigenvectors, 00084 *> backtransformed by the matrices in VR and/or VL; 00085 *> = 'S': compute selected right and/or left eigenvectors, 00086 *> specified by the logical array SELECT. 00087 *> \endverbatim 00088 *> 00089 *> \param[in] SELECT 00090 *> \verbatim 00091 *> SELECT is LOGICAL array, dimension (N) 00092 *> If HOWMNY='S', SELECT specifies the eigenvectors to be 00093 *> computed. If w(j) is a real eigenvalue, the corresponding 00094 *> real eigenvector is computed if SELECT(j) is .TRUE.. 00095 *> If w(j) and w(j+1) are the real and imaginary parts of a 00096 *> complex eigenvalue, the corresponding complex eigenvector 00097 *> is computed if either SELECT(j) or SELECT(j+1) is .TRUE., 00098 *> and on exit SELECT(j) is set to .TRUE. and SELECT(j+1) is 00099 *> set to .FALSE.. 00100 *> Not referenced if HOWMNY = 'A' or 'B'. 00101 *> \endverbatim 00102 *> 00103 *> \param[in] N 00104 *> \verbatim 00105 *> N is INTEGER 00106 *> The order of the matrices S and P. N >= 0. 00107 *> \endverbatim 00108 *> 00109 *> \param[in] S 00110 *> \verbatim 00111 *> S is DOUBLE PRECISION array, dimension (LDS,N) 00112 *> The upper quasi-triangular matrix S from a generalized Schur 00113 *> factorization, as computed by DHGEQZ. 00114 *> \endverbatim 00115 *> 00116 *> \param[in] LDS 00117 *> \verbatim 00118 *> LDS is INTEGER 00119 *> The leading dimension of array S. LDS >= max(1,N). 00120 *> \endverbatim 00121 *> 00122 *> \param[in] P 00123 *> \verbatim 00124 *> P is DOUBLE PRECISION array, dimension (LDP,N) 00125 *> The upper triangular matrix P from a generalized Schur 00126 *> factorization, as computed by DHGEQZ. 00127 *> 2-by-2 diagonal blocks of P corresponding to 2-by-2 blocks 00128 *> of S must be in positive diagonal form. 00129 *> \endverbatim 00130 *> 00131 *> \param[in] LDP 00132 *> \verbatim 00133 *> LDP is INTEGER 00134 *> The leading dimension of array P. LDP >= max(1,N). 00135 *> \endverbatim 00136 *> 00137 *> \param[in,out] VL 00138 *> \verbatim 00139 *> VL is DOUBLE PRECISION array, dimension (LDVL,MM) 00140 *> On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must 00141 *> contain an N-by-N matrix Q (usually the orthogonal matrix Q 00142 *> of left Schur vectors returned by DHGEQZ). 00143 *> On exit, if SIDE = 'L' or 'B', VL contains: 00144 *> if HOWMNY = 'A', the matrix Y of left eigenvectors of (S,P); 00145 *> if HOWMNY = 'B', the matrix Q*Y; 00146 *> if HOWMNY = 'S', the left eigenvectors of (S,P) specified by 00147 *> SELECT, stored consecutively in the columns of 00148 *> VL, in the same order as their eigenvalues. 00149 *> 00150 *> A complex eigenvector corresponding to a complex eigenvalue 00151 *> is stored in two consecutive columns, the first holding the 00152 *> real part, and the second the imaginary part. 00153 *> 00154 *> Not referenced if SIDE = 'R'. 00155 *> \endverbatim 00156 *> 00157 *> \param[in] LDVL 00158 *> \verbatim 00159 *> LDVL is INTEGER 00160 *> The leading dimension of array VL. LDVL >= 1, and if 00161 *> SIDE = 'L' or 'B', LDVL >= N. 00162 *> \endverbatim 00163 *> 00164 *> \param[in,out] VR 00165 *> \verbatim 00166 *> VR is DOUBLE PRECISION array, dimension (LDVR,MM) 00167 *> On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must 00168 *> contain an N-by-N matrix Z (usually the orthogonal matrix Z 00169 *> of right Schur vectors returned by DHGEQZ). 00170 *> 00171 *> On exit, if SIDE = 'R' or 'B', VR contains: 00172 *> if HOWMNY = 'A', the matrix X of right eigenvectors of (S,P); 00173 *> if HOWMNY = 'B' or 'b', the matrix Z*X; 00174 *> if HOWMNY = 'S' or 's', the right eigenvectors of (S,P) 00175 *> specified by SELECT, stored consecutively in the 00176 *> columns of VR, in the same order as their 00177 *> eigenvalues. 00178 *> 00179 *> A complex eigenvector corresponding to a complex eigenvalue 00180 *> is stored in two consecutive columns, the first holding the 00181 *> real part and the second the imaginary part. 00182 *> 00183 *> Not referenced if SIDE = 'L'. 00184 *> \endverbatim 00185 *> 00186 *> \param[in] LDVR 00187 *> \verbatim 00188 *> LDVR is INTEGER 00189 *> The leading dimension of the array VR. LDVR >= 1, and if 00190 *> SIDE = 'R' or 'B', LDVR >= N. 00191 *> \endverbatim 00192 *> 00193 *> \param[in] MM 00194 *> \verbatim 00195 *> MM is INTEGER 00196 *> The number of columns in the arrays VL and/or VR. MM >= M. 00197 *> \endverbatim 00198 *> 00199 *> \param[out] M 00200 *> \verbatim 00201 *> M is INTEGER 00202 *> The number of columns in the arrays VL and/or VR actually 00203 *> used to store the eigenvectors. If HOWMNY = 'A' or 'B', M 00204 *> is set to N. Each selected real eigenvector occupies one 00205 *> column and each selected complex eigenvector occupies two 00206 *> columns. 00207 *> \endverbatim 00208 *> 00209 *> \param[out] WORK 00210 *> \verbatim 00211 *> WORK is DOUBLE PRECISION array, dimension (6*N) 00212 *> \endverbatim 00213 *> 00214 *> \param[out] INFO 00215 *> \verbatim 00216 *> INFO is INTEGER 00217 *> = 0: successful exit. 00218 *> < 0: if INFO = -i, the i-th argument had an illegal value. 00219 *> > 0: the 2-by-2 block (INFO:INFO+1) does not have a complex 00220 *> eigenvalue. 00221 *> \endverbatim 00222 * 00223 * Authors: 00224 * ======== 00225 * 00226 *> \author Univ. of Tennessee 00227 *> \author Univ. of California Berkeley 00228 *> \author Univ. of Colorado Denver 00229 *> \author NAG Ltd. 00230 * 00231 *> \date November 2011 00232 * 00233 *> \ingroup doubleGEcomputational 00234 * 00235 *> \par Further Details: 00236 * ===================== 00237 *> 00238 *> \verbatim 00239 *> 00240 *> Allocation of workspace: 00241 *> ---------- -- --------- 00242 *> 00243 *> WORK( j ) = 1-norm of j-th column of A, above the diagonal 00244 *> WORK( N+j ) = 1-norm of j-th column of B, above the diagonal 00245 *> WORK( 2*N+1:3*N ) = real part of eigenvector 00246 *> WORK( 3*N+1:4*N ) = imaginary part of eigenvector 00247 *> WORK( 4*N+1:5*N ) = real part of back-transformed eigenvector 00248 *> WORK( 5*N+1:6*N ) = imaginary part of back-transformed eigenvector 00249 *> 00250 *> Rowwise vs. columnwise solution methods: 00251 *> ------- -- ---------- -------- ------- 00252 *> 00253 *> Finding a generalized eigenvector consists basically of solving the 00254 *> singular triangular system 00255 *> 00256 *> (A - w B) x = 0 (for right) or: (A - w B)**H y = 0 (for left) 00257 *> 00258 *> Consider finding the i-th right eigenvector (assume all eigenvalues 00259 *> are real). The equation to be solved is: 00260 *> n i 00261 *> 0 = sum C(j,k) v(k) = sum C(j,k) v(k) for j = i,. . .,1 00262 *> k=j k=j 00263 *> 00264 *> where C = (A - w B) (The components v(i+1:n) are 0.) 00265 *> 00266 *> The "rowwise" method is: 00267 *> 00268 *> (1) v(i) := 1 00269 *> for j = i-1,. . .,1: 00270 *> i 00271 *> (2) compute s = - sum C(j,k) v(k) and 00272 *> k=j+1 00273 *> 00274 *> (3) v(j) := s / C(j,j) 00275 *> 00276 *> Step 2 is sometimes called the "dot product" step, since it is an 00277 *> inner product between the j-th row and the portion of the eigenvector 00278 *> that has been computed so far. 00279 *> 00280 *> The "columnwise" method consists basically in doing the sums 00281 *> for all the rows in parallel. As each v(j) is computed, the 00282 *> contribution of v(j) times the j-th column of C is added to the 00283 *> partial sums. Since FORTRAN arrays are stored columnwise, this has 00284 *> the advantage that at each step, the elements of C that are accessed 00285 *> are adjacent to one another, whereas with the rowwise method, the 00286 *> elements accessed at a step are spaced LDS (and LDP) words apart. 00287 *> 00288 *> When finding left eigenvectors, the matrix in question is the 00289 *> transpose of the one in storage, so the rowwise method then 00290 *> actually accesses columns of A and B at each step, and so is the 00291 *> preferred method. 00292 *> \endverbatim 00293 *> 00294 * ===================================================================== 00295 SUBROUTINE DTGEVC( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL, 00296 $ LDVL, VR, LDVR, MM, M, WORK, INFO ) 00297 * 00298 * -- LAPACK computational routine (version 3.4.0) -- 00299 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00300 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00301 * November 2011 00302 * 00303 * .. Scalar Arguments .. 00304 CHARACTER HOWMNY, SIDE 00305 INTEGER INFO, LDP, LDS, LDVL, LDVR, M, MM, N 00306 * .. 00307 * .. Array Arguments .. 00308 LOGICAL SELECT( * ) 00309 DOUBLE PRECISION P( LDP, * ), S( LDS, * ), VL( LDVL, * ), 00310 $ VR( LDVR, * ), WORK( * ) 00311 * .. 00312 * 00313 * 00314 * ===================================================================== 00315 * 00316 * .. Parameters .. 00317 DOUBLE PRECISION ZERO, ONE, SAFETY 00318 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0, 00319 $ SAFETY = 1.0D+2 ) 00320 * .. 00321 * .. Local Scalars .. 00322 LOGICAL COMPL, COMPR, IL2BY2, ILABAD, ILALL, ILBACK, 00323 $ ILBBAD, ILCOMP, ILCPLX, LSA, LSB 00324 INTEGER I, IBEG, IEIG, IEND, IHWMNY, IINFO, IM, ISIDE, 00325 $ J, JA, JC, JE, JR, JW, NA, NW 00326 DOUBLE PRECISION ACOEF, ACOEFA, ANORM, ASCALE, BCOEFA, BCOEFI, 00327 $ BCOEFR, BIG, BIGNUM, BNORM, BSCALE, CIM2A, 00328 $ CIM2B, CIMAGA, CIMAGB, CRE2A, CRE2B, CREALA, 00329 $ CREALB, DMIN, SAFMIN, SALFAR, SBETA, SCALE, 00330 $ SMALL, TEMP, TEMP2, TEMP2I, TEMP2R, ULP, XMAX, 00331 $ XSCALE 00332 * .. 00333 * .. Local Arrays .. 00334 DOUBLE PRECISION BDIAG( 2 ), SUM( 2, 2 ), SUMS( 2, 2 ), 00335 $ SUMP( 2, 2 ) 00336 * .. 00337 * .. External Functions .. 00338 LOGICAL LSAME 00339 DOUBLE PRECISION DLAMCH 00340 EXTERNAL LSAME, DLAMCH 00341 * .. 00342 * .. External Subroutines .. 00343 EXTERNAL DGEMV, DLABAD, DLACPY, DLAG2, DLALN2, XERBLA 00344 * .. 00345 * .. Intrinsic Functions .. 00346 INTRINSIC ABS, MAX, MIN 00347 * .. 00348 * .. Executable Statements .. 00349 * 00350 * Decode and Test the input parameters 00351 * 00352 IF( LSAME( HOWMNY, 'A' ) ) THEN 00353 IHWMNY = 1 00354 ILALL = .TRUE. 00355 ILBACK = .FALSE. 00356 ELSE IF( LSAME( HOWMNY, 'S' ) ) THEN 00357 IHWMNY = 2 00358 ILALL = .FALSE. 00359 ILBACK = .FALSE. 00360 ELSE IF( LSAME( HOWMNY, 'B' ) ) THEN 00361 IHWMNY = 3 00362 ILALL = .TRUE. 00363 ILBACK = .TRUE. 00364 ELSE 00365 IHWMNY = -1 00366 ILALL = .TRUE. 00367 END IF 00368 * 00369 IF( LSAME( SIDE, 'R' ) ) THEN 00370 ISIDE = 1 00371 COMPL = .FALSE. 00372 COMPR = .TRUE. 00373 ELSE IF( LSAME( SIDE, 'L' ) ) THEN 00374 ISIDE = 2 00375 COMPL = .TRUE. 00376 COMPR = .FALSE. 00377 ELSE IF( LSAME( SIDE, 'B' ) ) THEN 00378 ISIDE = 3 00379 COMPL = .TRUE. 00380 COMPR = .TRUE. 00381 ELSE 00382 ISIDE = -1 00383 END IF 00384 * 00385 INFO = 0 00386 IF( ISIDE.LT.0 ) THEN 00387 INFO = -1 00388 ELSE IF( IHWMNY.LT.0 ) THEN 00389 INFO = -2 00390 ELSE IF( N.LT.0 ) THEN 00391 INFO = -4 00392 ELSE IF( LDS.LT.MAX( 1, N ) ) THEN 00393 INFO = -6 00394 ELSE IF( LDP.LT.MAX( 1, N ) ) THEN 00395 INFO = -8 00396 END IF 00397 IF( INFO.NE.0 ) THEN 00398 CALL XERBLA( 'DTGEVC', -INFO ) 00399 RETURN 00400 END IF 00401 * 00402 * Count the number of eigenvectors to be computed 00403 * 00404 IF( .NOT.ILALL ) THEN 00405 IM = 0 00406 ILCPLX = .FALSE. 00407 DO 10 J = 1, N 00408 IF( ILCPLX ) THEN 00409 ILCPLX = .FALSE. 00410 GO TO 10 00411 END IF 00412 IF( J.LT.N ) THEN 00413 IF( S( J+1, J ).NE.ZERO ) 00414 $ ILCPLX = .TRUE. 00415 END IF 00416 IF( ILCPLX ) THEN 00417 IF( SELECT( J ) .OR. SELECT( J+1 ) ) 00418 $ IM = IM + 2 00419 ELSE 00420 IF( SELECT( J ) ) 00421 $ IM = IM + 1 00422 END IF 00423 10 CONTINUE 00424 ELSE 00425 IM = N 00426 END IF 00427 * 00428 * Check 2-by-2 diagonal blocks of A, B 00429 * 00430 ILABAD = .FALSE. 00431 ILBBAD = .FALSE. 00432 DO 20 J = 1, N - 1 00433 IF( S( J+1, J ).NE.ZERO ) THEN 00434 IF( P( J, J ).EQ.ZERO .OR. P( J+1, J+1 ).EQ.ZERO .OR. 00435 $ P( J, J+1 ).NE.ZERO )ILBBAD = .TRUE. 00436 IF( J.LT.N-1 ) THEN 00437 IF( S( J+2, J+1 ).NE.ZERO ) 00438 $ ILABAD = .TRUE. 00439 END IF 00440 END IF 00441 20 CONTINUE 00442 * 00443 IF( ILABAD ) THEN 00444 INFO = -5 00445 ELSE IF( ILBBAD ) THEN 00446 INFO = -7 00447 ELSE IF( COMPL .AND. LDVL.LT.N .OR. LDVL.LT.1 ) THEN 00448 INFO = -10 00449 ELSE IF( COMPR .AND. LDVR.LT.N .OR. LDVR.LT.1 ) THEN 00450 INFO = -12 00451 ELSE IF( MM.LT.IM ) THEN 00452 INFO = -13 00453 END IF 00454 IF( INFO.NE.0 ) THEN 00455 CALL XERBLA( 'DTGEVC', -INFO ) 00456 RETURN 00457 END IF 00458 * 00459 * Quick return if possible 00460 * 00461 M = IM 00462 IF( N.EQ.0 ) 00463 $ RETURN 00464 * 00465 * Machine Constants 00466 * 00467 SAFMIN = DLAMCH( 'Safe minimum' ) 00468 BIG = ONE / SAFMIN 00469 CALL DLABAD( SAFMIN, BIG ) 00470 ULP = DLAMCH( 'Epsilon' )*DLAMCH( 'Base' ) 00471 SMALL = SAFMIN*N / ULP 00472 BIG = ONE / SMALL 00473 BIGNUM = ONE / ( SAFMIN*N ) 00474 * 00475 * Compute the 1-norm of each column of the strictly upper triangular 00476 * part (i.e., excluding all elements belonging to the diagonal 00477 * blocks) of A and B to check for possible overflow in the 00478 * triangular solver. 00479 * 00480 ANORM = ABS( S( 1, 1 ) ) 00481 IF( N.GT.1 ) 00482 $ ANORM = ANORM + ABS( S( 2, 1 ) ) 00483 BNORM = ABS( P( 1, 1 ) ) 00484 WORK( 1 ) = ZERO 00485 WORK( N+1 ) = ZERO 00486 * 00487 DO 50 J = 2, N 00488 TEMP = ZERO 00489 TEMP2 = ZERO 00490 IF( S( J, J-1 ).EQ.ZERO ) THEN 00491 IEND = J - 1 00492 ELSE 00493 IEND = J - 2 00494 END IF 00495 DO 30 I = 1, IEND 00496 TEMP = TEMP + ABS( S( I, J ) ) 00497 TEMP2 = TEMP2 + ABS( P( I, J ) ) 00498 30 CONTINUE 00499 WORK( J ) = TEMP 00500 WORK( N+J ) = TEMP2 00501 DO 40 I = IEND + 1, MIN( J+1, N ) 00502 TEMP = TEMP + ABS( S( I, J ) ) 00503 TEMP2 = TEMP2 + ABS( P( I, J ) ) 00504 40 CONTINUE 00505 ANORM = MAX( ANORM, TEMP ) 00506 BNORM = MAX( BNORM, TEMP2 ) 00507 50 CONTINUE 00508 * 00509 ASCALE = ONE / MAX( ANORM, SAFMIN ) 00510 BSCALE = ONE / MAX( BNORM, SAFMIN ) 00511 * 00512 * Left eigenvectors 00513 * 00514 IF( COMPL ) THEN 00515 IEIG = 0 00516 * 00517 * Main loop over eigenvalues 00518 * 00519 ILCPLX = .FALSE. 00520 DO 220 JE = 1, N 00521 * 00522 * Skip this iteration if (a) HOWMNY='S' and SELECT=.FALSE., or 00523 * (b) this would be the second of a complex pair. 00524 * Check for complex eigenvalue, so as to be sure of which 00525 * entry(-ies) of SELECT to look at. 00526 * 00527 IF( ILCPLX ) THEN 00528 ILCPLX = .FALSE. 00529 GO TO 220 00530 END IF 00531 NW = 1 00532 IF( JE.LT.N ) THEN 00533 IF( S( JE+1, JE ).NE.ZERO ) THEN 00534 ILCPLX = .TRUE. 00535 NW = 2 00536 END IF 00537 END IF 00538 IF( ILALL ) THEN 00539 ILCOMP = .TRUE. 00540 ELSE IF( ILCPLX ) THEN 00541 ILCOMP = SELECT( JE ) .OR. SELECT( JE+1 ) 00542 ELSE 00543 ILCOMP = SELECT( JE ) 00544 END IF 00545 IF( .NOT.ILCOMP ) 00546 $ GO TO 220 00547 * 00548 * Decide if (a) singular pencil, (b) real eigenvalue, or 00549 * (c) complex eigenvalue. 00550 * 00551 IF( .NOT.ILCPLX ) THEN 00552 IF( ABS( S( JE, JE ) ).LE.SAFMIN .AND. 00553 $ ABS( P( JE, JE ) ).LE.SAFMIN ) THEN 00554 * 00555 * Singular matrix pencil -- return unit eigenvector 00556 * 00557 IEIG = IEIG + 1 00558 DO 60 JR = 1, N 00559 VL( JR, IEIG ) = ZERO 00560 60 CONTINUE 00561 VL( IEIG, IEIG ) = ONE 00562 GO TO 220 00563 END IF 00564 END IF 00565 * 00566 * Clear vector 00567 * 00568 DO 70 JR = 1, NW*N 00569 WORK( 2*N+JR ) = ZERO 00570 70 CONTINUE 00571 * T 00572 * Compute coefficients in ( a A - b B ) y = 0 00573 * a is ACOEF 00574 * b is BCOEFR + i*BCOEFI 00575 * 00576 IF( .NOT.ILCPLX ) THEN 00577 * 00578 * Real eigenvalue 00579 * 00580 TEMP = ONE / MAX( ABS( S( JE, JE ) )*ASCALE, 00581 $ ABS( P( JE, JE ) )*BSCALE, SAFMIN ) 00582 SALFAR = ( TEMP*S( JE, JE ) )*ASCALE 00583 SBETA = ( TEMP*P( JE, JE ) )*BSCALE 00584 ACOEF = SBETA*ASCALE 00585 BCOEFR = SALFAR*BSCALE 00586 BCOEFI = ZERO 00587 * 00588 * Scale to avoid underflow 00589 * 00590 SCALE = ONE 00591 LSA = ABS( SBETA ).GE.SAFMIN .AND. ABS( ACOEF ).LT.SMALL 00592 LSB = ABS( SALFAR ).GE.SAFMIN .AND. ABS( BCOEFR ).LT. 00593 $ SMALL 00594 IF( LSA ) 00595 $ SCALE = ( SMALL / ABS( SBETA ) )*MIN( ANORM, BIG ) 00596 IF( LSB ) 00597 $ SCALE = MAX( SCALE, ( SMALL / ABS( SALFAR ) )* 00598 $ MIN( BNORM, BIG ) ) 00599 IF( LSA .OR. LSB ) THEN 00600 SCALE = MIN( SCALE, ONE / 00601 $ ( SAFMIN*MAX( ONE, ABS( ACOEF ), 00602 $ ABS( BCOEFR ) ) ) ) 00603 IF( LSA ) THEN 00604 ACOEF = ASCALE*( SCALE*SBETA ) 00605 ELSE 00606 ACOEF = SCALE*ACOEF 00607 END IF 00608 IF( LSB ) THEN 00609 BCOEFR = BSCALE*( SCALE*SALFAR ) 00610 ELSE 00611 BCOEFR = SCALE*BCOEFR 00612 END IF 00613 END IF 00614 ACOEFA = ABS( ACOEF ) 00615 BCOEFA = ABS( BCOEFR ) 00616 * 00617 * First component is 1 00618 * 00619 WORK( 2*N+JE ) = ONE 00620 XMAX = ONE 00621 ELSE 00622 * 00623 * Complex eigenvalue 00624 * 00625 CALL DLAG2( S( JE, JE ), LDS, P( JE, JE ), LDP, 00626 $ SAFMIN*SAFETY, ACOEF, TEMP, BCOEFR, TEMP2, 00627 $ BCOEFI ) 00628 BCOEFI = -BCOEFI 00629 IF( BCOEFI.EQ.ZERO ) THEN 00630 INFO = JE 00631 RETURN 00632 END IF 00633 * 00634 * Scale to avoid over/underflow 00635 * 00636 ACOEFA = ABS( ACOEF ) 00637 BCOEFA = ABS( BCOEFR ) + ABS( BCOEFI ) 00638 SCALE = ONE 00639 IF( ACOEFA*ULP.LT.SAFMIN .AND. ACOEFA.GE.SAFMIN ) 00640 $ SCALE = ( SAFMIN / ULP ) / ACOEFA 00641 IF( BCOEFA*ULP.LT.SAFMIN .AND. BCOEFA.GE.SAFMIN ) 00642 $ SCALE = MAX( SCALE, ( SAFMIN / ULP ) / BCOEFA ) 00643 IF( SAFMIN*ACOEFA.GT.ASCALE ) 00644 $ SCALE = ASCALE / ( SAFMIN*ACOEFA ) 00645 IF( SAFMIN*BCOEFA.GT.BSCALE ) 00646 $ SCALE = MIN( SCALE, BSCALE / ( SAFMIN*BCOEFA ) ) 00647 IF( SCALE.NE.ONE ) THEN 00648 ACOEF = SCALE*ACOEF 00649 ACOEFA = ABS( ACOEF ) 00650 BCOEFR = SCALE*BCOEFR 00651 BCOEFI = SCALE*BCOEFI 00652 BCOEFA = ABS( BCOEFR ) + ABS( BCOEFI ) 00653 END IF 00654 * 00655 * Compute first two components of eigenvector 00656 * 00657 TEMP = ACOEF*S( JE+1, JE ) 00658 TEMP2R = ACOEF*S( JE, JE ) - BCOEFR*P( JE, JE ) 00659 TEMP2I = -BCOEFI*P( JE, JE ) 00660 IF( ABS( TEMP ).GT.ABS( TEMP2R )+ABS( TEMP2I ) ) THEN 00661 WORK( 2*N+JE ) = ONE 00662 WORK( 3*N+JE ) = ZERO 00663 WORK( 2*N+JE+1 ) = -TEMP2R / TEMP 00664 WORK( 3*N+JE+1 ) = -TEMP2I / TEMP 00665 ELSE 00666 WORK( 2*N+JE+1 ) = ONE 00667 WORK( 3*N+JE+1 ) = ZERO 00668 TEMP = ACOEF*S( JE, JE+1 ) 00669 WORK( 2*N+JE ) = ( BCOEFR*P( JE+1, JE+1 )-ACOEF* 00670 $ S( JE+1, JE+1 ) ) / TEMP 00671 WORK( 3*N+JE ) = BCOEFI*P( JE+1, JE+1 ) / TEMP 00672 END IF 00673 XMAX = MAX( ABS( WORK( 2*N+JE ) )+ABS( WORK( 3*N+JE ) ), 00674 $ ABS( WORK( 2*N+JE+1 ) )+ABS( WORK( 3*N+JE+1 ) ) ) 00675 END IF 00676 * 00677 DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN ) 00678 * 00679 * T 00680 * Triangular solve of (a A - b B) y = 0 00681 * 00682 * T 00683 * (rowwise in (a A - b B) , or columnwise in (a A - b B) ) 00684 * 00685 IL2BY2 = .FALSE. 00686 * 00687 DO 160 J = JE + NW, N 00688 IF( IL2BY2 ) THEN 00689 IL2BY2 = .FALSE. 00690 GO TO 160 00691 END IF 00692 * 00693 NA = 1 00694 BDIAG( 1 ) = P( J, J ) 00695 IF( J.LT.N ) THEN 00696 IF( S( J+1, J ).NE.ZERO ) THEN 00697 IL2BY2 = .TRUE. 00698 BDIAG( 2 ) = P( J+1, J+1 ) 00699 NA = 2 00700 END IF 00701 END IF 00702 * 00703 * Check whether scaling is necessary for dot products 00704 * 00705 XSCALE = ONE / MAX( ONE, XMAX ) 00706 TEMP = MAX( WORK( J ), WORK( N+J ), 00707 $ ACOEFA*WORK( J )+BCOEFA*WORK( N+J ) ) 00708 IF( IL2BY2 ) 00709 $ TEMP = MAX( TEMP, WORK( J+1 ), WORK( N+J+1 ), 00710 $ ACOEFA*WORK( J+1 )+BCOEFA*WORK( N+J+1 ) ) 00711 IF( TEMP.GT.BIGNUM*XSCALE ) THEN 00712 DO 90 JW = 0, NW - 1 00713 DO 80 JR = JE, J - 1 00714 WORK( ( JW+2 )*N+JR ) = XSCALE* 00715 $ WORK( ( JW+2 )*N+JR ) 00716 80 CONTINUE 00717 90 CONTINUE 00718 XMAX = XMAX*XSCALE 00719 END IF 00720 * 00721 * Compute dot products 00722 * 00723 * j-1 00724 * SUM = sum conjg( a*S(k,j) - b*P(k,j) )*x(k) 00725 * k=je 00726 * 00727 * To reduce the op count, this is done as 00728 * 00729 * _ j-1 _ j-1 00730 * a*conjg( sum S(k,j)*x(k) ) - b*conjg( sum P(k,j)*x(k) ) 00731 * k=je k=je 00732 * 00733 * which may cause underflow problems if A or B are close 00734 * to underflow. (E.g., less than SMALL.) 00735 * 00736 * 00737 DO 120 JW = 1, NW 00738 DO 110 JA = 1, NA 00739 SUMS( JA, JW ) = ZERO 00740 SUMP( JA, JW ) = ZERO 00741 * 00742 DO 100 JR = JE, J - 1 00743 SUMS( JA, JW ) = SUMS( JA, JW ) + 00744 $ S( JR, J+JA-1 )* 00745 $ WORK( ( JW+1 )*N+JR ) 00746 SUMP( JA, JW ) = SUMP( JA, JW ) + 00747 $ P( JR, J+JA-1 )* 00748 $ WORK( ( JW+1 )*N+JR ) 00749 100 CONTINUE 00750 110 CONTINUE 00751 120 CONTINUE 00752 * 00753 DO 130 JA = 1, NA 00754 IF( ILCPLX ) THEN 00755 SUM( JA, 1 ) = -ACOEF*SUMS( JA, 1 ) + 00756 $ BCOEFR*SUMP( JA, 1 ) - 00757 $ BCOEFI*SUMP( JA, 2 ) 00758 SUM( JA, 2 ) = -ACOEF*SUMS( JA, 2 ) + 00759 $ BCOEFR*SUMP( JA, 2 ) + 00760 $ BCOEFI*SUMP( JA, 1 ) 00761 ELSE 00762 SUM( JA, 1 ) = -ACOEF*SUMS( JA, 1 ) + 00763 $ BCOEFR*SUMP( JA, 1 ) 00764 END IF 00765 130 CONTINUE 00766 * 00767 * T 00768 * Solve ( a A - b B ) y = SUM(,) 00769 * with scaling and perturbation of the denominator 00770 * 00771 CALL DLALN2( .TRUE., NA, NW, DMIN, ACOEF, S( J, J ), LDS, 00772 $ BDIAG( 1 ), BDIAG( 2 ), SUM, 2, BCOEFR, 00773 $ BCOEFI, WORK( 2*N+J ), N, SCALE, TEMP, 00774 $ IINFO ) 00775 IF( SCALE.LT.ONE ) THEN 00776 DO 150 JW = 0, NW - 1 00777 DO 140 JR = JE, J - 1 00778 WORK( ( JW+2 )*N+JR ) = SCALE* 00779 $ WORK( ( JW+2 )*N+JR ) 00780 140 CONTINUE 00781 150 CONTINUE 00782 XMAX = SCALE*XMAX 00783 END IF 00784 XMAX = MAX( XMAX, TEMP ) 00785 160 CONTINUE 00786 * 00787 * Copy eigenvector to VL, back transforming if 00788 * HOWMNY='B'. 00789 * 00790 IEIG = IEIG + 1 00791 IF( ILBACK ) THEN 00792 DO 170 JW = 0, NW - 1 00793 CALL DGEMV( 'N', N, N+1-JE, ONE, VL( 1, JE ), LDVL, 00794 $ WORK( ( JW+2 )*N+JE ), 1, ZERO, 00795 $ WORK( ( JW+4 )*N+1 ), 1 ) 00796 170 CONTINUE 00797 CALL DLACPY( ' ', N, NW, WORK( 4*N+1 ), N, VL( 1, JE ), 00798 $ LDVL ) 00799 IBEG = 1 00800 ELSE 00801 CALL DLACPY( ' ', N, NW, WORK( 2*N+1 ), N, VL( 1, IEIG ), 00802 $ LDVL ) 00803 IBEG = JE 00804 END IF 00805 * 00806 * Scale eigenvector 00807 * 00808 XMAX = ZERO 00809 IF( ILCPLX ) THEN 00810 DO 180 J = IBEG, N 00811 XMAX = MAX( XMAX, ABS( VL( J, IEIG ) )+ 00812 $ ABS( VL( J, IEIG+1 ) ) ) 00813 180 CONTINUE 00814 ELSE 00815 DO 190 J = IBEG, N 00816 XMAX = MAX( XMAX, ABS( VL( J, IEIG ) ) ) 00817 190 CONTINUE 00818 END IF 00819 * 00820 IF( XMAX.GT.SAFMIN ) THEN 00821 XSCALE = ONE / XMAX 00822 * 00823 DO 210 JW = 0, NW - 1 00824 DO 200 JR = IBEG, N 00825 VL( JR, IEIG+JW ) = XSCALE*VL( JR, IEIG+JW ) 00826 200 CONTINUE 00827 210 CONTINUE 00828 END IF 00829 IEIG = IEIG + NW - 1 00830 * 00831 220 CONTINUE 00832 END IF 00833 * 00834 * Right eigenvectors 00835 * 00836 IF( COMPR ) THEN 00837 IEIG = IM + 1 00838 * 00839 * Main loop over eigenvalues 00840 * 00841 ILCPLX = .FALSE. 00842 DO 500 JE = N, 1, -1 00843 * 00844 * Skip this iteration if (a) HOWMNY='S' and SELECT=.FALSE., or 00845 * (b) this would be the second of a complex pair. 00846 * Check for complex eigenvalue, so as to be sure of which 00847 * entry(-ies) of SELECT to look at -- if complex, SELECT(JE) 00848 * or SELECT(JE-1). 00849 * If this is a complex pair, the 2-by-2 diagonal block 00850 * corresponding to the eigenvalue is in rows/columns JE-1:JE 00851 * 00852 IF( ILCPLX ) THEN 00853 ILCPLX = .FALSE. 00854 GO TO 500 00855 END IF 00856 NW = 1 00857 IF( JE.GT.1 ) THEN 00858 IF( S( JE, JE-1 ).NE.ZERO ) THEN 00859 ILCPLX = .TRUE. 00860 NW = 2 00861 END IF 00862 END IF 00863 IF( ILALL ) THEN 00864 ILCOMP = .TRUE. 00865 ELSE IF( ILCPLX ) THEN 00866 ILCOMP = SELECT( JE ) .OR. SELECT( JE-1 ) 00867 ELSE 00868 ILCOMP = SELECT( JE ) 00869 END IF 00870 IF( .NOT.ILCOMP ) 00871 $ GO TO 500 00872 * 00873 * Decide if (a) singular pencil, (b) real eigenvalue, or 00874 * (c) complex eigenvalue. 00875 * 00876 IF( .NOT.ILCPLX ) THEN 00877 IF( ABS( S( JE, JE ) ).LE.SAFMIN .AND. 00878 $ ABS( P( JE, JE ) ).LE.SAFMIN ) THEN 00879 * 00880 * Singular matrix pencil -- unit eigenvector 00881 * 00882 IEIG = IEIG - 1 00883 DO 230 JR = 1, N 00884 VR( JR, IEIG ) = ZERO 00885 230 CONTINUE 00886 VR( IEIG, IEIG ) = ONE 00887 GO TO 500 00888 END IF 00889 END IF 00890 * 00891 * Clear vector 00892 * 00893 DO 250 JW = 0, NW - 1 00894 DO 240 JR = 1, N 00895 WORK( ( JW+2 )*N+JR ) = ZERO 00896 240 CONTINUE 00897 250 CONTINUE 00898 * 00899 * Compute coefficients in ( a A - b B ) x = 0 00900 * a is ACOEF 00901 * b is BCOEFR + i*BCOEFI 00902 * 00903 IF( .NOT.ILCPLX ) THEN 00904 * 00905 * Real eigenvalue 00906 * 00907 TEMP = ONE / MAX( ABS( S( JE, JE ) )*ASCALE, 00908 $ ABS( P( JE, JE ) )*BSCALE, SAFMIN ) 00909 SALFAR = ( TEMP*S( JE, JE ) )*ASCALE 00910 SBETA = ( TEMP*P( JE, JE ) )*BSCALE 00911 ACOEF = SBETA*ASCALE 00912 BCOEFR = SALFAR*BSCALE 00913 BCOEFI = ZERO 00914 * 00915 * Scale to avoid underflow 00916 * 00917 SCALE = ONE 00918 LSA = ABS( SBETA ).GE.SAFMIN .AND. ABS( ACOEF ).LT.SMALL 00919 LSB = ABS( SALFAR ).GE.SAFMIN .AND. ABS( BCOEFR ).LT. 00920 $ SMALL 00921 IF( LSA ) 00922 $ SCALE = ( SMALL / ABS( SBETA ) )*MIN( ANORM, BIG ) 00923 IF( LSB ) 00924 $ SCALE = MAX( SCALE, ( SMALL / ABS( SALFAR ) )* 00925 $ MIN( BNORM, BIG ) ) 00926 IF( LSA .OR. LSB ) THEN 00927 SCALE = MIN( SCALE, ONE / 00928 $ ( SAFMIN*MAX( ONE, ABS( ACOEF ), 00929 $ ABS( BCOEFR ) ) ) ) 00930 IF( LSA ) THEN 00931 ACOEF = ASCALE*( SCALE*SBETA ) 00932 ELSE 00933 ACOEF = SCALE*ACOEF 00934 END IF 00935 IF( LSB ) THEN 00936 BCOEFR = BSCALE*( SCALE*SALFAR ) 00937 ELSE 00938 BCOEFR = SCALE*BCOEFR 00939 END IF 00940 END IF 00941 ACOEFA = ABS( ACOEF ) 00942 BCOEFA = ABS( BCOEFR ) 00943 * 00944 * First component is 1 00945 * 00946 WORK( 2*N+JE ) = ONE 00947 XMAX = ONE 00948 * 00949 * Compute contribution from column JE of A and B to sum 00950 * (See "Further Details", above.) 00951 * 00952 DO 260 JR = 1, JE - 1 00953 WORK( 2*N+JR ) = BCOEFR*P( JR, JE ) - 00954 $ ACOEF*S( JR, JE ) 00955 260 CONTINUE 00956 ELSE 00957 * 00958 * Complex eigenvalue 00959 * 00960 CALL DLAG2( S( JE-1, JE-1 ), LDS, P( JE-1, JE-1 ), LDP, 00961 $ SAFMIN*SAFETY, ACOEF, TEMP, BCOEFR, TEMP2, 00962 $ BCOEFI ) 00963 IF( BCOEFI.EQ.ZERO ) THEN 00964 INFO = JE - 1 00965 RETURN 00966 END IF 00967 * 00968 * Scale to avoid over/underflow 00969 * 00970 ACOEFA = ABS( ACOEF ) 00971 BCOEFA = ABS( BCOEFR ) + ABS( BCOEFI ) 00972 SCALE = ONE 00973 IF( ACOEFA*ULP.LT.SAFMIN .AND. ACOEFA.GE.SAFMIN ) 00974 $ SCALE = ( SAFMIN / ULP ) / ACOEFA 00975 IF( BCOEFA*ULP.LT.SAFMIN .AND. BCOEFA.GE.SAFMIN ) 00976 $ SCALE = MAX( SCALE, ( SAFMIN / ULP ) / BCOEFA ) 00977 IF( SAFMIN*ACOEFA.GT.ASCALE ) 00978 $ SCALE = ASCALE / ( SAFMIN*ACOEFA ) 00979 IF( SAFMIN*BCOEFA.GT.BSCALE ) 00980 $ SCALE = MIN( SCALE, BSCALE / ( SAFMIN*BCOEFA ) ) 00981 IF( SCALE.NE.ONE ) THEN 00982 ACOEF = SCALE*ACOEF 00983 ACOEFA = ABS( ACOEF ) 00984 BCOEFR = SCALE*BCOEFR 00985 BCOEFI = SCALE*BCOEFI 00986 BCOEFA = ABS( BCOEFR ) + ABS( BCOEFI ) 00987 END IF 00988 * 00989 * Compute first two components of eigenvector 00990 * and contribution to sums 00991 * 00992 TEMP = ACOEF*S( JE, JE-1 ) 00993 TEMP2R = ACOEF*S( JE, JE ) - BCOEFR*P( JE, JE ) 00994 TEMP2I = -BCOEFI*P( JE, JE ) 00995 IF( ABS( TEMP ).GE.ABS( TEMP2R )+ABS( TEMP2I ) ) THEN 00996 WORK( 2*N+JE ) = ONE 00997 WORK( 3*N+JE ) = ZERO 00998 WORK( 2*N+JE-1 ) = -TEMP2R / TEMP 00999 WORK( 3*N+JE-1 ) = -TEMP2I / TEMP 01000 ELSE 01001 WORK( 2*N+JE-1 ) = ONE 01002 WORK( 3*N+JE-1 ) = ZERO 01003 TEMP = ACOEF*S( JE-1, JE ) 01004 WORK( 2*N+JE ) = ( BCOEFR*P( JE-1, JE-1 )-ACOEF* 01005 $ S( JE-1, JE-1 ) ) / TEMP 01006 WORK( 3*N+JE ) = BCOEFI*P( JE-1, JE-1 ) / TEMP 01007 END IF 01008 * 01009 XMAX = MAX( ABS( WORK( 2*N+JE ) )+ABS( WORK( 3*N+JE ) ), 01010 $ ABS( WORK( 2*N+JE-1 ) )+ABS( WORK( 3*N+JE-1 ) ) ) 01011 * 01012 * Compute contribution from columns JE and JE-1 01013 * of A and B to the sums. 01014 * 01015 CREALA = ACOEF*WORK( 2*N+JE-1 ) 01016 CIMAGA = ACOEF*WORK( 3*N+JE-1 ) 01017 CREALB = BCOEFR*WORK( 2*N+JE-1 ) - 01018 $ BCOEFI*WORK( 3*N+JE-1 ) 01019 CIMAGB = BCOEFI*WORK( 2*N+JE-1 ) + 01020 $ BCOEFR*WORK( 3*N+JE-1 ) 01021 CRE2A = ACOEF*WORK( 2*N+JE ) 01022 CIM2A = ACOEF*WORK( 3*N+JE ) 01023 CRE2B = BCOEFR*WORK( 2*N+JE ) - BCOEFI*WORK( 3*N+JE ) 01024 CIM2B = BCOEFI*WORK( 2*N+JE ) + BCOEFR*WORK( 3*N+JE ) 01025 DO 270 JR = 1, JE - 2 01026 WORK( 2*N+JR ) = -CREALA*S( JR, JE-1 ) + 01027 $ CREALB*P( JR, JE-1 ) - 01028 $ CRE2A*S( JR, JE ) + CRE2B*P( JR, JE ) 01029 WORK( 3*N+JR ) = -CIMAGA*S( JR, JE-1 ) + 01030 $ CIMAGB*P( JR, JE-1 ) - 01031 $ CIM2A*S( JR, JE ) + CIM2B*P( JR, JE ) 01032 270 CONTINUE 01033 END IF 01034 * 01035 DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN ) 01036 * 01037 * Columnwise triangular solve of (a A - b B) x = 0 01038 * 01039 IL2BY2 = .FALSE. 01040 DO 370 J = JE - NW, 1, -1 01041 * 01042 * If a 2-by-2 block, is in position j-1:j, wait until 01043 * next iteration to process it (when it will be j:j+1) 01044 * 01045 IF( .NOT.IL2BY2 .AND. J.GT.1 ) THEN 01046 IF( S( J, J-1 ).NE.ZERO ) THEN 01047 IL2BY2 = .TRUE. 01048 GO TO 370 01049 END IF 01050 END IF 01051 BDIAG( 1 ) = P( J, J ) 01052 IF( IL2BY2 ) THEN 01053 NA = 2 01054 BDIAG( 2 ) = P( J+1, J+1 ) 01055 ELSE 01056 NA = 1 01057 END IF 01058 * 01059 * Compute x(j) (and x(j+1), if 2-by-2 block) 01060 * 01061 CALL DLALN2( .FALSE., NA, NW, DMIN, ACOEF, S( J, J ), 01062 $ LDS, BDIAG( 1 ), BDIAG( 2 ), WORK( 2*N+J ), 01063 $ N, BCOEFR, BCOEFI, SUM, 2, SCALE, TEMP, 01064 $ IINFO ) 01065 IF( SCALE.LT.ONE ) THEN 01066 * 01067 DO 290 JW = 0, NW - 1 01068 DO 280 JR = 1, JE 01069 WORK( ( JW+2 )*N+JR ) = SCALE* 01070 $ WORK( ( JW+2 )*N+JR ) 01071 280 CONTINUE 01072 290 CONTINUE 01073 END IF 01074 XMAX = MAX( SCALE*XMAX, TEMP ) 01075 * 01076 DO 310 JW = 1, NW 01077 DO 300 JA = 1, NA 01078 WORK( ( JW+1 )*N+J+JA-1 ) = SUM( JA, JW ) 01079 300 CONTINUE 01080 310 CONTINUE 01081 * 01082 * w = w + x(j)*(a S(*,j) - b P(*,j) ) with scaling 01083 * 01084 IF( J.GT.1 ) THEN 01085 * 01086 * Check whether scaling is necessary for sum. 01087 * 01088 XSCALE = ONE / MAX( ONE, XMAX ) 01089 TEMP = ACOEFA*WORK( J ) + BCOEFA*WORK( N+J ) 01090 IF( IL2BY2 ) 01091 $ TEMP = MAX( TEMP, ACOEFA*WORK( J+1 )+BCOEFA* 01092 $ WORK( N+J+1 ) ) 01093 TEMP = MAX( TEMP, ACOEFA, BCOEFA ) 01094 IF( TEMP.GT.BIGNUM*XSCALE ) THEN 01095 * 01096 DO 330 JW = 0, NW - 1 01097 DO 320 JR = 1, JE 01098 WORK( ( JW+2 )*N+JR ) = XSCALE* 01099 $ WORK( ( JW+2 )*N+JR ) 01100 320 CONTINUE 01101 330 CONTINUE 01102 XMAX = XMAX*XSCALE 01103 END IF 01104 * 01105 * Compute the contributions of the off-diagonals of 01106 * column j (and j+1, if 2-by-2 block) of A and B to the 01107 * sums. 01108 * 01109 * 01110 DO 360 JA = 1, NA 01111 IF( ILCPLX ) THEN 01112 CREALA = ACOEF*WORK( 2*N+J+JA-1 ) 01113 CIMAGA = ACOEF*WORK( 3*N+J+JA-1 ) 01114 CREALB = BCOEFR*WORK( 2*N+J+JA-1 ) - 01115 $ BCOEFI*WORK( 3*N+J+JA-1 ) 01116 CIMAGB = BCOEFI*WORK( 2*N+J+JA-1 ) + 01117 $ BCOEFR*WORK( 3*N+J+JA-1 ) 01118 DO 340 JR = 1, J - 1 01119 WORK( 2*N+JR ) = WORK( 2*N+JR ) - 01120 $ CREALA*S( JR, J+JA-1 ) + 01121 $ CREALB*P( JR, J+JA-1 ) 01122 WORK( 3*N+JR ) = WORK( 3*N+JR ) - 01123 $ CIMAGA*S( JR, J+JA-1 ) + 01124 $ CIMAGB*P( JR, J+JA-1 ) 01125 340 CONTINUE 01126 ELSE 01127 CREALA = ACOEF*WORK( 2*N+J+JA-1 ) 01128 CREALB = BCOEFR*WORK( 2*N+J+JA-1 ) 01129 DO 350 JR = 1, J - 1 01130 WORK( 2*N+JR ) = WORK( 2*N+JR ) - 01131 $ CREALA*S( JR, J+JA-1 ) + 01132 $ CREALB*P( JR, J+JA-1 ) 01133 350 CONTINUE 01134 END IF 01135 360 CONTINUE 01136 END IF 01137 * 01138 IL2BY2 = .FALSE. 01139 370 CONTINUE 01140 * 01141 * Copy eigenvector to VR, back transforming if 01142 * HOWMNY='B'. 01143 * 01144 IEIG = IEIG - NW 01145 IF( ILBACK ) THEN 01146 * 01147 DO 410 JW = 0, NW - 1 01148 DO 380 JR = 1, N 01149 WORK( ( JW+4 )*N+JR ) = WORK( ( JW+2 )*N+1 )* 01150 $ VR( JR, 1 ) 01151 380 CONTINUE 01152 * 01153 * A series of compiler directives to defeat 01154 * vectorization for the next loop 01155 * 01156 * 01157 DO 400 JC = 2, JE 01158 DO 390 JR = 1, N 01159 WORK( ( JW+4 )*N+JR ) = WORK( ( JW+4 )*N+JR ) + 01160 $ WORK( ( JW+2 )*N+JC )*VR( JR, JC ) 01161 390 CONTINUE 01162 400 CONTINUE 01163 410 CONTINUE 01164 * 01165 DO 430 JW = 0, NW - 1 01166 DO 420 JR = 1, N 01167 VR( JR, IEIG+JW ) = WORK( ( JW+4 )*N+JR ) 01168 420 CONTINUE 01169 430 CONTINUE 01170 * 01171 IEND = N 01172 ELSE 01173 DO 450 JW = 0, NW - 1 01174 DO 440 JR = 1, N 01175 VR( JR, IEIG+JW ) = WORK( ( JW+2 )*N+JR ) 01176 440 CONTINUE 01177 450 CONTINUE 01178 * 01179 IEND = JE 01180 END IF 01181 * 01182 * Scale eigenvector 01183 * 01184 XMAX = ZERO 01185 IF( ILCPLX ) THEN 01186 DO 460 J = 1, IEND 01187 XMAX = MAX( XMAX, ABS( VR( J, IEIG ) )+ 01188 $ ABS( VR( J, IEIG+1 ) ) ) 01189 460 CONTINUE 01190 ELSE 01191 DO 470 J = 1, IEND 01192 XMAX = MAX( XMAX, ABS( VR( J, IEIG ) ) ) 01193 470 CONTINUE 01194 END IF 01195 * 01196 IF( XMAX.GT.SAFMIN ) THEN 01197 XSCALE = ONE / XMAX 01198 DO 490 JW = 0, NW - 1 01199 DO 480 JR = 1, IEND 01200 VR( JR, IEIG+JW ) = XSCALE*VR( JR, IEIG+JW ) 01201 480 CONTINUE 01202 490 CONTINUE 01203 END IF 01204 500 CONTINUE 01205 END IF 01206 * 01207 RETURN 01208 * 01209 * End of DTGEVC 01210 * 01211 END