ergo
|
00001 /* Ergo, version 3.2, a program for linear scaling electronic structure 00002 * calculations. 00003 * Copyright (C) 2012 Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek. 00004 * 00005 * This program is free software: you can redistribute it and/or modify 00006 * it under the terms of the GNU General Public License as published by 00007 * the Free Software Foundation, either version 3 of the License, or 00008 * (at your option) any later version. 00009 * 00010 * This program is distributed in the hope that it will be useful, 00011 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00012 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00013 * GNU General Public License for more details. 00014 * 00015 * You should have received a copy of the GNU General Public License 00016 * along with this program. If not, see <http://www.gnu.org/licenses/>. 00017 * 00018 * Primary academic reference: 00019 * KohnâSham Density Functional Theory Electronic Structure Calculations 00020 * with Linearly Scaling Computational Time and Memory Usage, 00021 * Elias Rudberg, Emanuel H. Rubensson, and Pawel Salek, 00022 * J. Chem. Theory Comput. 7, 340 (2011), 00023 * <http://dx.doi.org/10.1021/ct100611z> 00024 * 00025 * For further information about Ergo, see <http://www.ergoscf.org>. 00026 */ 00027 00028 /* This file belongs to the template_lapack part of the Ergo source 00029 * code. The source files in the template_lapack directory are modified 00030 * versions of files originally distributed as CLAPACK, see the 00031 * Copyright/license notice in the file template_lapack/COPYING. 00032 */ 00033 00034 00035 #ifndef TEMPLATE_LAPACK_TGEVC_HEADER 00036 #define TEMPLATE_LAPACK_TGEVC_HEADER 00037 00038 00039 template<class Treal> 00040 int template_lapack_tgevc(const char *side, const char *howmny, const logical *select, 00041 const integer *n, const Treal *a, const integer *lda, const Treal *b, const integer *ldb, 00042 Treal *vl, const integer *ldvl, Treal *vr, const integer *ldvr, const integer 00043 *mm, integer *m, Treal *work, integer *info) 00044 { 00045 /* -- LAPACK routine (version 3.0) -- 00046 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 00047 Courant Institute, Argonne National Lab, and Rice University 00048 June 30, 1999 00049 00050 00051 00052 Purpose 00053 ======= 00054 00055 DTGEVC computes some or all of the right and/or left generalized 00056 eigenvectors of a pair of real upper triangular matrices (A,B). 00057 00058 The right generalized eigenvector x and the left generalized 00059 eigenvector y of (A,B) corresponding to a generalized eigenvalue 00060 w are defined by: 00061 00062 (A - wB) * x = 0 and y**H * (A - wB) = 0 00063 00064 where y**H denotes the conjugate tranpose of y. 00065 00066 If an eigenvalue w is determined by zero diagonal elements of both A 00067 and B, a unit vector is returned as the corresponding eigenvector. 00068 00069 If all eigenvectors are requested, the routine may either return 00070 the matrices X and/or Y of right or left eigenvectors of (A,B), or 00071 the products Z*X and/or Q*Y, where Z and Q are input orthogonal 00072 matrices. If (A,B) was obtained from the generalized real-Schur 00073 factorization of an original pair of matrices 00074 (A0,B0) = (Q*A*Z**H,Q*B*Z**H), 00075 then Z*X and Q*Y are the matrices of right or left eigenvectors of 00076 A. 00077 00078 A must be block upper triangular, with 1-by-1 and 2-by-2 diagonal 00079 blocks. Corresponding to each 2-by-2 diagonal block is a complex 00080 conjugate pair of eigenvalues and eigenvectors; only one 00081 eigenvector of the pair is computed, namely the one corresponding 00082 to the eigenvalue with positive imaginary part. 00083 00084 Arguments 00085 ========= 00086 00087 SIDE (input) CHARACTER*1 00088 = 'R': compute right eigenvectors only; 00089 = 'L': compute left eigenvectors only; 00090 = 'B': compute both right and left eigenvectors. 00091 00092 HOWMNY (input) CHARACTER*1 00093 = 'A': compute all right and/or left eigenvectors; 00094 = 'B': compute all right and/or left eigenvectors, and 00095 backtransform them using the input matrices supplied 00096 in VR and/or VL; 00097 = 'S': compute selected right and/or left eigenvectors, 00098 specified by the logical array SELECT. 00099 00100 SELECT (input) LOGICAL array, dimension (N) 00101 If HOWMNY='S', SELECT specifies the eigenvectors to be 00102 computed. 00103 If HOWMNY='A' or 'B', SELECT is not referenced. 00104 To select the real eigenvector corresponding to the real 00105 eigenvalue w(j), SELECT(j) must be set to .TRUE. To select 00106 the complex eigenvector corresponding to a complex conjugate 00107 pair w(j) and w(j+1), either SELECT(j) or SELECT(j+1) must 00108 be set to .TRUE.. 00109 00110 N (input) INTEGER 00111 The order of the matrices A and B. N >= 0. 00112 00113 A (input) DOUBLE PRECISION array, dimension (LDA,N) 00114 The upper quasi-triangular matrix A. 00115 00116 LDA (input) INTEGER 00117 The leading dimension of array A. LDA >= max(1, N). 00118 00119 B (input) DOUBLE PRECISION array, dimension (LDB,N) 00120 The upper triangular matrix B. If A has a 2-by-2 diagonal 00121 block, then the corresponding 2-by-2 block of B must be 00122 diagonal with positive elements. 00123 00124 LDB (input) INTEGER 00125 The leading dimension of array B. LDB >= max(1,N). 00126 00127 VL (input/output) DOUBLE PRECISION array, dimension (LDVL,MM) 00128 On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must 00129 contain an N-by-N matrix Q (usually the orthogonal matrix Q 00130 of left Schur vectors returned by DHGEQZ). 00131 On exit, if SIDE = 'L' or 'B', VL contains: 00132 if HOWMNY = 'A', the matrix Y of left eigenvectors of (A,B); 00133 if HOWMNY = 'B', the matrix Q*Y; 00134 if HOWMNY = 'S', the left eigenvectors of (A,B) specified by 00135 SELECT, stored consecutively in the columns of 00136 VL, in the same order as their eigenvalues. 00137 If SIDE = 'R', VL is not referenced. 00138 00139 A complex eigenvector corresponding to a complex eigenvalue 00140 is stored in two consecutive columns, the first holding the 00141 real part, and the second the imaginary part. 00142 00143 LDVL (input) INTEGER 00144 The leading dimension of array VL. 00145 LDVL >= max(1,N) if SIDE = 'L' or 'B'; LDVL >= 1 otherwise. 00146 00147 VR (input/output) DOUBLE PRECISION array, dimension (LDVR,MM) 00148 On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must 00149 contain an N-by-N matrix Q (usually the orthogonal matrix Z 00150 of right Schur vectors returned by DHGEQZ). 00151 On exit, if SIDE = 'R' or 'B', VR contains: 00152 if HOWMNY = 'A', the matrix X of right eigenvectors of (A,B); 00153 if HOWMNY = 'B', the matrix Z*X; 00154 if HOWMNY = 'S', the right eigenvectors of (A,B) specified by 00155 SELECT, stored consecutively in the columns of 00156 VR, in the same order as their eigenvalues. 00157 If SIDE = 'L', VR is not referenced. 00158 00159 A complex eigenvector corresponding to a complex eigenvalue 00160 is stored in two consecutive columns, the first holding the 00161 real part and the second the imaginary part. 00162 00163 LDVR (input) INTEGER 00164 The leading dimension of the array VR. 00165 LDVR >= max(1,N) if SIDE = 'R' or 'B'; LDVR >= 1 otherwise. 00166 00167 MM (input) INTEGER 00168 The number of columns in the arrays VL and/or VR. MM >= M. 00169 00170 M (output) INTEGER 00171 The number of columns in the arrays VL and/or VR actually 00172 used to store the eigenvectors. If HOWMNY = 'A' or 'B', M 00173 is set to N. Each selected real eigenvector occupies one 00174 column and each selected complex eigenvector occupies two 00175 columns. 00176 00177 WORK (workspace) DOUBLE PRECISION array, dimension (6*N) 00178 00179 INFO (output) INTEGER 00180 = 0: successful exit. 00181 < 0: if INFO = -i, the i-th argument had an illegal value. 00182 > 0: the 2-by-2 block (INFO:INFO+1) does not have a complex 00183 eigenvalue. 00184 00185 Further Details 00186 =============== 00187 00188 Allocation of workspace: 00189 ---------- -- --------- 00190 00191 WORK( j ) = 1-norm of j-th column of A, above the diagonal 00192 WORK( N+j ) = 1-norm of j-th column of B, above the diagonal 00193 WORK( 2*N+1:3*N ) = real part of eigenvector 00194 WORK( 3*N+1:4*N ) = imaginary part of eigenvector 00195 WORK( 4*N+1:5*N ) = real part of back-transformed eigenvector 00196 WORK( 5*N+1:6*N ) = imaginary part of back-transformed eigenvector 00197 00198 Rowwise vs. columnwise solution methods: 00199 ------- -- ---------- -------- ------- 00200 00201 Finding a generalized eigenvector consists basically of solving the 00202 singular triangular system 00203 00204 (A - w B) x = 0 (for right) or: (A - w B)**H y = 0 (for left) 00205 00206 Consider finding the i-th right eigenvector (assume all eigenvalues 00207 are real). The equation to be solved is: 00208 n i 00209 0 = sum C(j,k) v(k) = sum C(j,k) v(k) for j = i,. . .,1 00210 k=j k=j 00211 00212 where C = (A - w B) (The components v(i+1:n) are 0.) 00213 00214 The "rowwise" method is: 00215 00216 (1) v(i) := 1 00217 for j = i-1,. . .,1: 00218 i 00219 (2) compute s = - sum C(j,k) v(k) and 00220 k=j+1 00221 00222 (3) v(j) := s / C(j,j) 00223 00224 Step 2 is sometimes called the "dot product" step, since it is an 00225 inner product between the j-th row and the portion of the eigenvector 00226 that has been computed so far. 00227 00228 The "columnwise" method consists basically in doing the sums 00229 for all the rows in parallel. As each v(j) is computed, the 00230 contribution of v(j) times the j-th column of C is added to the 00231 partial sums. Since FORTRAN arrays are stored columnwise, this has 00232 the advantage that at each step, the elements of C that are accessed 00233 are adjacent to one another, whereas with the rowwise method, the 00234 elements accessed at a step are spaced LDA (and LDB) words apart. 00235 00236 When finding left eigenvectors, the matrix in question is the 00237 transpose of the one in storage, so the rowwise method then 00238 actually accesses columns of A and B at each step, and so is the 00239 preferred method. 00240 00241 ===================================================================== 00242 00243 00244 Decode and Test the input parameters 00245 00246 Parameter adjustments */ 00247 /* Table of constant values */ 00248 logical c_true = TRUE_; 00249 integer c__2 = 2; 00250 Treal c_b35 = 1.; 00251 integer c__1 = 1; 00252 Treal c_b37 = 0.; 00253 logical c_false = FALSE_; 00254 00255 /* System generated locals */ 00256 integer a_dim1, a_offset, b_dim1, b_offset, vl_dim1, vl_offset, vr_dim1, 00257 vr_offset, i__1, i__2, i__3, i__4, i__5; 00258 Treal d__1, d__2, d__3, d__4, d__5, d__6; 00259 /* Local variables */ 00260 integer ibeg, ieig, iend; 00261 Treal dmin__, temp, suma[4] /* was [2][2] */, sumb[4] 00262 /* was [2][2] */, xmax; 00263 Treal cim2a, cim2b, cre2a, cre2b, temp2, bdiag[2]; 00264 integer i__, j; 00265 Treal acoef, scale; 00266 logical ilall; 00267 integer iside; 00268 Treal sbeta; 00269 logical il2by2; 00270 integer iinfo; 00271 Treal small; 00272 logical compl_AAAA; 00273 Treal anorm, bnorm; 00274 logical compr; 00275 Treal temp2i; 00276 Treal temp2r; 00277 integer ja; 00278 logical ilabad, ilbbad; 00279 integer jc, je, na; 00280 Treal acoefa, bcoefa, cimaga, cimagb; 00281 logical ilback; 00282 integer im; 00283 Treal bcoefi, ascale, bscale, creala; 00284 integer jr; 00285 Treal crealb; 00286 Treal bcoefr; 00287 integer jw, nw; 00288 Treal salfar, safmin; 00289 Treal xscale, bignum; 00290 logical ilcomp, ilcplx; 00291 integer ihwmny; 00292 Treal big; 00293 logical lsa, lsb; 00294 Treal ulp, sum[4] /* was [2][2] */; 00295 #define suma_ref(a_1,a_2) suma[(a_2)*2 + a_1 - 3] 00296 #define sumb_ref(a_1,a_2) sumb[(a_2)*2 + a_1 - 3] 00297 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1] 00298 #define b_ref(a_1,a_2) b[(a_2)*b_dim1 + a_1] 00299 #define vl_ref(a_1,a_2) vl[(a_2)*vl_dim1 + a_1] 00300 #define vr_ref(a_1,a_2) vr[(a_2)*vr_dim1 + a_1] 00301 #define sum_ref(a_1,a_2) sum[(a_2)*2 + a_1 - 3] 00302 00303 00304 --select; 00305 a_dim1 = *lda; 00306 a_offset = 1 + a_dim1 * 1; 00307 a -= a_offset; 00308 b_dim1 = *ldb; 00309 b_offset = 1 + b_dim1 * 1; 00310 b -= b_offset; 00311 vl_dim1 = *ldvl; 00312 vl_offset = 1 + vl_dim1 * 1; 00313 vl -= vl_offset; 00314 vr_dim1 = *ldvr; 00315 vr_offset = 1 + vr_dim1 * 1; 00316 vr -= vr_offset; 00317 --work; 00318 00319 /* Initialization added by Elias to get rid of compiler warnings. */ 00320 ilback = 0; 00321 /* Function Body */ 00322 if (template_blas_lsame(howmny, "A")) { 00323 ihwmny = 1; 00324 ilall = TRUE_; 00325 ilback = FALSE_; 00326 } else if (template_blas_lsame(howmny, "S")) { 00327 ihwmny = 2; 00328 ilall = FALSE_; 00329 ilback = FALSE_; 00330 } else if (template_blas_lsame(howmny, "B") || template_blas_lsame(howmny, 00331 "T")) { 00332 ihwmny = 3; 00333 ilall = TRUE_; 00334 ilback = TRUE_; 00335 } else { 00336 ihwmny = -1; 00337 ilall = TRUE_; 00338 } 00339 00340 if (template_blas_lsame(side, "R")) { 00341 iside = 1; 00342 compl_AAAA = FALSE_; 00343 compr = TRUE_; 00344 } else if (template_blas_lsame(side, "L")) { 00345 iside = 2; 00346 compl_AAAA = TRUE_; 00347 compr = FALSE_; 00348 } else if (template_blas_lsame(side, "B")) { 00349 iside = 3; 00350 compl_AAAA = TRUE_; 00351 compr = TRUE_; 00352 } else { 00353 iside = -1; 00354 } 00355 00356 *info = 0; 00357 if (iside < 0) { 00358 *info = -1; 00359 } else if (ihwmny < 0) { 00360 *info = -2; 00361 } else if (*n < 0) { 00362 *info = -4; 00363 } else if (*lda < maxMACRO(1,*n)) { 00364 *info = -6; 00365 } else if (*ldb < maxMACRO(1,*n)) { 00366 *info = -8; 00367 } 00368 if (*info != 0) { 00369 i__1 = -(*info); 00370 template_blas_erbla("TGEVC ", &i__1); 00371 return 0; 00372 } 00373 00374 /* Count the number of eigenvectors to be computed */ 00375 00376 if (! ilall) { 00377 im = 0; 00378 ilcplx = FALSE_; 00379 i__1 = *n; 00380 for (j = 1; j <= i__1; ++j) { 00381 if (ilcplx) { 00382 ilcplx = FALSE_; 00383 goto L10; 00384 } 00385 if (j < *n) { 00386 if (a_ref(j + 1, j) != 0.) { 00387 ilcplx = TRUE_; 00388 } 00389 } 00390 if (ilcplx) { 00391 if (select[j] || select[j + 1]) { 00392 im += 2; 00393 } 00394 } else { 00395 if (select[j]) { 00396 ++im; 00397 } 00398 } 00399 L10: 00400 ; 00401 } 00402 } else { 00403 im = *n; 00404 } 00405 00406 /* Check 2-by-2 diagonal blocks of A, B */ 00407 00408 ilabad = FALSE_; 00409 ilbbad = FALSE_; 00410 i__1 = *n - 1; 00411 for (j = 1; j <= i__1; ++j) { 00412 if (a_ref(j + 1, j) != 0.) { 00413 if (b_ref(j, j) == 0. || b_ref(j + 1, j + 1) == 0. || b_ref(j, j 00414 + 1) != 0.) { 00415 ilbbad = TRUE_; 00416 } 00417 if (j < *n - 1) { 00418 if (a_ref(j + 2, j + 1) != 0.) { 00419 ilabad = TRUE_; 00420 } 00421 } 00422 } 00423 /* L20: */ 00424 } 00425 00426 if (ilabad) { 00427 *info = -5; 00428 } else if (ilbbad) { 00429 *info = -7; 00430 } else if ( ( compl_AAAA && *ldvl < *n ) || *ldvl < 1) { 00431 *info = -10; 00432 } else if ( ( compr && *ldvr < *n ) || *ldvr < 1) { 00433 *info = -12; 00434 } else if (*mm < im) { 00435 *info = -13; 00436 } 00437 if (*info != 0) { 00438 i__1 = -(*info); 00439 template_blas_erbla("TGEVC ", &i__1); 00440 return 0; 00441 } 00442 00443 /* Quick return if possible */ 00444 00445 *m = im; 00446 if (*n == 0) { 00447 return 0; 00448 } 00449 00450 /* Machine Constants */ 00451 00452 safmin = template_lapack_lamch("Safe minimum", (Treal)0); 00453 big = 1. / safmin; 00454 template_lapack_labad(&safmin, &big); 00455 ulp = template_lapack_lamch("Epsilon", (Treal)0) * template_lapack_lamch("Base", (Treal)0); 00456 small = safmin * *n / ulp; 00457 big = 1. / small; 00458 bignum = 1. / (safmin * *n); 00459 00460 /* Compute the 1-norm of each column of the strictly upper triangular 00461 part (i.e., excluding all elements belonging to the diagonal 00462 blocks) of A and B to check for possible overflow in the 00463 triangular solver. */ 00464 00465 anorm = (d__1 = a_ref(1, 1), absMACRO(d__1)); 00466 if (*n > 1) { 00467 anorm += (d__1 = a_ref(2, 1), absMACRO(d__1)); 00468 } 00469 bnorm = (d__1 = b_ref(1, 1), absMACRO(d__1)); 00470 work[1] = 0.; 00471 work[*n + 1] = 0.; 00472 00473 i__1 = *n; 00474 for (j = 2; j <= i__1; ++j) { 00475 temp = 0.; 00476 temp2 = 0.; 00477 if (a_ref(j, j - 1) == 0.) { 00478 iend = j - 1; 00479 } else { 00480 iend = j - 2; 00481 } 00482 i__2 = iend; 00483 for (i__ = 1; i__ <= i__2; ++i__) { 00484 temp += (d__1 = a_ref(i__, j), absMACRO(d__1)); 00485 temp2 += (d__1 = b_ref(i__, j), absMACRO(d__1)); 00486 /* L30: */ 00487 } 00488 work[j] = temp; 00489 work[*n + j] = temp2; 00490 /* Computing MIN */ 00491 i__3 = j + 1; 00492 i__2 = minMACRO(i__3,*n); 00493 for (i__ = iend + 1; i__ <= i__2; ++i__) { 00494 temp += (d__1 = a_ref(i__, j), absMACRO(d__1)); 00495 temp2 += (d__1 = b_ref(i__, j), absMACRO(d__1)); 00496 /* L40: */ 00497 } 00498 anorm = maxMACRO(anorm,temp); 00499 bnorm = maxMACRO(bnorm,temp2); 00500 /* L50: */ 00501 } 00502 00503 ascale = 1. / maxMACRO(anorm,safmin); 00504 bscale = 1. / maxMACRO(bnorm,safmin); 00505 00506 /* Left eigenvectors */ 00507 00508 if (compl_AAAA) { 00509 ieig = 0; 00510 00511 /* Main loop over eigenvalues */ 00512 00513 ilcplx = FALSE_; 00514 i__1 = *n; 00515 for (je = 1; je <= i__1; ++je) { 00516 00517 /* Skip this iteration if (a) HOWMNY='S' and SELECT=.FALSE., or 00518 (b) this would be the second of a complex pair. 00519 Check for complex eigenvalue, so as to be sure of which 00520 entry(-ies) of SELECT to look at. */ 00521 00522 if (ilcplx) { 00523 ilcplx = FALSE_; 00524 goto L220; 00525 } 00526 nw = 1; 00527 if (je < *n) { 00528 if (a_ref(je + 1, je) != 0.) { 00529 ilcplx = TRUE_; 00530 nw = 2; 00531 } 00532 } 00533 if (ilall) { 00534 ilcomp = TRUE_; 00535 } else if (ilcplx) { 00536 ilcomp = select[je] || select[je + 1]; 00537 } else { 00538 ilcomp = select[je]; 00539 } 00540 if (! ilcomp) { 00541 goto L220; 00542 } 00543 00544 /* Decide if (a) singular pencil, (b) real eigenvalue, or 00545 (c) complex eigenvalue. */ 00546 00547 if (! ilcplx) { 00548 if ((d__1 = a_ref(je, je), absMACRO(d__1)) <= safmin && (d__2 = 00549 b_ref(je, je), absMACRO(d__2)) <= safmin) { 00550 00551 /* Singular matrix pencil -- return unit eigenvector */ 00552 00553 ++ieig; 00554 i__2 = *n; 00555 for (jr = 1; jr <= i__2; ++jr) { 00556 vl_ref(jr, ieig) = 0.; 00557 /* L60: */ 00558 } 00559 vl_ref(ieig, ieig) = 1.; 00560 goto L220; 00561 } 00562 } 00563 00564 /* Clear vector */ 00565 00566 i__2 = nw * *n; 00567 for (jr = 1; jr <= i__2; ++jr) { 00568 work[(*n << 1) + jr] = 0.; 00569 /* L70: */ 00570 } 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 (! ilcplx) { 00577 00578 /* Real eigenvalue 00579 00580 Computing MAX */ 00581 d__3 = (d__1 = a_ref(je, je), absMACRO(d__1)) * ascale, d__4 = ( 00582 d__2 = b_ref(je, je), absMACRO(d__2)) * bscale, d__3 = maxMACRO( 00583 d__3,d__4); 00584 temp = 1. / maxMACRO(d__3,safmin); 00585 salfar = temp * a_ref(je, je) * ascale; 00586 sbeta = temp * b_ref(je, je) * bscale; 00587 acoef = sbeta * ascale; 00588 bcoefr = salfar * bscale; 00589 bcoefi = 0.; 00590 00591 /* Scale to avoid underflow */ 00592 00593 scale = 1.; 00594 lsa = absMACRO(sbeta) >= safmin && absMACRO(acoef) < small; 00595 lsb = absMACRO(salfar) >= safmin && absMACRO(bcoefr) < small; 00596 if (lsa) { 00597 scale = small / absMACRO(sbeta) * minMACRO(anorm,big); 00598 } 00599 if (lsb) { 00600 /* Computing MAX */ 00601 d__1 = scale, d__2 = small / absMACRO(salfar) * minMACRO(bnorm,big); 00602 scale = maxMACRO(d__1,d__2); 00603 } 00604 if (lsa || lsb) { 00605 /* Computing MIN 00606 Computing MAX */ 00607 d__3 = 1., d__4 = absMACRO(acoef), d__3 = maxMACRO(d__3,d__4), d__4 00608 = absMACRO(bcoefr); 00609 d__1 = scale, d__2 = 1. / (safmin * maxMACRO(d__3,d__4)); 00610 scale = minMACRO(d__1,d__2); 00611 if (lsa) { 00612 acoef = ascale * (scale * sbeta); 00613 } else { 00614 acoef = scale * acoef; 00615 } 00616 if (lsb) { 00617 bcoefr = bscale * (scale * salfar); 00618 } else { 00619 bcoefr = scale * bcoefr; 00620 } 00621 } 00622 acoefa = absMACRO(acoef); 00623 bcoefa = absMACRO(bcoefr); 00624 00625 /* First component is 1 */ 00626 00627 work[(*n << 1) + je] = 1.; 00628 xmax = 1.; 00629 } else { 00630 00631 /* Complex eigenvalue */ 00632 00633 d__1 = safmin * 100.; 00634 template_lapack_lag2(&a_ref(je, je), lda, &b_ref(je, je), ldb, &d__1, & 00635 acoef, &temp, &bcoefr, &temp2, &bcoefi); 00636 bcoefi = -bcoefi; 00637 if (bcoefi == 0.) { 00638 *info = je; 00639 return 0; 00640 } 00641 00642 /* Scale to avoid over/underflow */ 00643 00644 acoefa = absMACRO(acoef); 00645 bcoefa = absMACRO(bcoefr) + absMACRO(bcoefi); 00646 scale = 1.; 00647 if (acoefa * ulp < safmin && acoefa >= safmin) { 00648 scale = safmin / ulp / acoefa; 00649 } 00650 if (bcoefa * ulp < safmin && bcoefa >= safmin) { 00651 /* Computing MAX */ 00652 d__1 = scale, d__2 = safmin / ulp / bcoefa; 00653 scale = maxMACRO(d__1,d__2); 00654 } 00655 if (safmin * acoefa > ascale) { 00656 scale = ascale / (safmin * acoefa); 00657 } 00658 if (safmin * bcoefa > bscale) { 00659 /* Computing MIN */ 00660 d__1 = scale, d__2 = bscale / (safmin * bcoefa); 00661 scale = minMACRO(d__1,d__2); 00662 } 00663 if (scale != 1.) { 00664 acoef = scale * acoef; 00665 acoefa = absMACRO(acoef); 00666 bcoefr = scale * bcoefr; 00667 bcoefi = scale * bcoefi; 00668 bcoefa = absMACRO(bcoefr) + absMACRO(bcoefi); 00669 } 00670 00671 /* Compute first two components of eigenvector */ 00672 00673 temp = acoef * a_ref(je + 1, je); 00674 temp2r = acoef * a_ref(je, je) - bcoefr * b_ref(je, je); 00675 temp2i = -bcoefi * b_ref(je, je); 00676 if (absMACRO(temp) > absMACRO(temp2r) + absMACRO(temp2i)) { 00677 work[(*n << 1) + je] = 1.; 00678 work[*n * 3 + je] = 0.; 00679 work[(*n << 1) + je + 1] = -temp2r / temp; 00680 work[*n * 3 + je + 1] = -temp2i / temp; 00681 } else { 00682 work[(*n << 1) + je + 1] = 1.; 00683 work[*n * 3 + je + 1] = 0.; 00684 temp = acoef * a_ref(je, je + 1); 00685 work[(*n << 1) + je] = (bcoefr * b_ref(je + 1, je + 1) - 00686 acoef * a_ref(je + 1, je + 1)) / temp; 00687 work[*n * 3 + je] = bcoefi * b_ref(je + 1, je + 1) / temp; 00688 } 00689 /* Computing MAX */ 00690 d__5 = (d__1 = work[(*n << 1) + je], absMACRO(d__1)) + (d__2 = 00691 work[*n * 3 + je], absMACRO(d__2)), d__6 = (d__3 = work[(* 00692 n << 1) + je + 1], absMACRO(d__3)) + (d__4 = work[*n * 3 + 00693 je + 1], absMACRO(d__4)); 00694 xmax = maxMACRO(d__5,d__6); 00695 } 00696 00697 /* Computing MAX */ 00698 d__1 = ulp * acoefa * anorm, d__2 = ulp * bcoefa * bnorm, d__1 = 00699 maxMACRO(d__1,d__2); 00700 dmin__ = maxMACRO(d__1,safmin); 00701 00702 /* T 00703 Triangular solve of (a A - b B) y = 0 00704 00705 T 00706 (rowwise in (a A - b B) , or columnwise in (a A - b B) ) */ 00707 00708 il2by2 = FALSE_; 00709 00710 i__2 = *n; 00711 for (j = je + nw; j <= i__2; ++j) { 00712 if (il2by2) { 00713 il2by2 = FALSE_; 00714 goto L160; 00715 } 00716 00717 na = 1; 00718 bdiag[0] = b_ref(j, j); 00719 if (j < *n) { 00720 if (a_ref(j + 1, j) != 0.) { 00721 il2by2 = TRUE_; 00722 bdiag[1] = b_ref(j + 1, j + 1); 00723 na = 2; 00724 } 00725 } 00726 00727 /* Check whether scaling is necessary for dot products */ 00728 00729 xscale = 1. / maxMACRO(1.,xmax); 00730 /* Computing MAX */ 00731 d__1 = work[j], d__2 = work[*n + j], d__1 = maxMACRO(d__1,d__2), 00732 d__2 = acoefa * work[j] + bcoefa * work[*n + j]; 00733 temp = maxMACRO(d__1,d__2); 00734 if (il2by2) { 00735 /* Computing MAX */ 00736 d__1 = temp, d__2 = work[j + 1], d__1 = maxMACRO(d__1,d__2), 00737 d__2 = work[*n + j + 1], d__1 = maxMACRO(d__1,d__2), 00738 d__2 = acoefa * work[j + 1] + bcoefa * work[*n + 00739 j + 1]; 00740 temp = maxMACRO(d__1,d__2); 00741 } 00742 if (temp > bignum * xscale) { 00743 i__3 = nw - 1; 00744 for (jw = 0; jw <= i__3; ++jw) { 00745 i__4 = j - 1; 00746 for (jr = je; jr <= i__4; ++jr) { 00747 work[(jw + 2) * *n + jr] = xscale * work[(jw + 2) 00748 * *n + jr]; 00749 /* L80: */ 00750 } 00751 /* L90: */ 00752 } 00753 xmax *= xscale; 00754 } 00755 00756 /* Compute dot products 00757 00758 j-1 00759 SUM = sum conjg( a*A(k,j) - b*B(k,j) )*x(k) 00760 k=je 00761 00762 To reduce the op count, this is done as 00763 00764 _ j-1 _ j-1 00765 a*conjg( sum A(k,j)*x(k) ) - b*conjg( sum B(k,j)*x(k) ) 00766 k=je k=je 00767 00768 which may cause underflow problems if A or B are close 00769 to underflow. (E.g., less than SMALL.) 00770 00771 00772 A series of compiler directives to defeat vectorization 00773 for the next loop 00774 00775 $PL$ CMCHAR=' ' 00776 DIR$ NEXTSCALAR 00777 $DIR SCALAR 00778 DIR$ NEXT SCALAR 00779 VD$L NOVECTOR 00780 DEC$ NOVECTOR 00781 VD$ NOVECTOR 00782 VDIR NOVECTOR 00783 VOCL LOOP,SCALAR 00784 IBM PREFER SCALAR 00785 $PL$ CMCHAR='*' */ 00786 00787 i__3 = nw; 00788 for (jw = 1; jw <= i__3; ++jw) { 00789 00790 /* $PL$ CMCHAR=' ' 00791 DIR$ NEXTSCALAR 00792 $DIR SCALAR 00793 DIR$ NEXT SCALAR 00794 VD$L NOVECTOR 00795 DEC$ NOVECTOR 00796 VD$ NOVECTOR 00797 VDIR NOVECTOR 00798 VOCL LOOP,SCALAR 00799 IBM PREFER SCALAR 00800 $PL$ CMCHAR='*' */ 00801 00802 i__4 = na; 00803 for (ja = 1; ja <= i__4; ++ja) { 00804 suma_ref(ja, jw) = 0.; 00805 sumb_ref(ja, jw) = 0.; 00806 00807 i__5 = j - 1; 00808 for (jr = je; jr <= i__5; ++jr) { 00809 suma_ref(ja, jw) = suma_ref(ja, jw) + a_ref(jr, j 00810 + ja - 1) * work[(jw + 1) * *n + jr]; 00811 sumb_ref(ja, jw) = sumb_ref(ja, jw) + b_ref(jr, j 00812 + ja - 1) * work[(jw + 1) * *n + jr]; 00813 /* L100: */ 00814 } 00815 /* L110: */ 00816 } 00817 /* L120: */ 00818 } 00819 00820 /* $PL$ CMCHAR=' ' 00821 DIR$ NEXTSCALAR 00822 $DIR SCALAR 00823 DIR$ NEXT SCALAR 00824 VD$L NOVECTOR 00825 DEC$ NOVECTOR 00826 VD$ NOVECTOR 00827 VDIR NOVECTOR 00828 VOCL LOOP,SCALAR 00829 IBM PREFER SCALAR 00830 $PL$ CMCHAR='*' */ 00831 00832 i__3 = na; 00833 for (ja = 1; ja <= i__3; ++ja) { 00834 if (ilcplx) { 00835 sum_ref(ja, 1) = -acoef * suma_ref(ja, 1) + bcoefr * 00836 sumb_ref(ja, 1) - bcoefi * sumb_ref(ja, 2); 00837 sum_ref(ja, 2) = -acoef * suma_ref(ja, 2) + bcoefr * 00838 sumb_ref(ja, 2) + bcoefi * sumb_ref(ja, 1); 00839 } else { 00840 sum_ref(ja, 1) = -acoef * suma_ref(ja, 1) + bcoefr * 00841 sumb_ref(ja, 1); 00842 } 00843 /* L130: */ 00844 } 00845 00846 /* T 00847 Solve ( a A - b B ) y = SUM(,) 00848 with scaling and perturbation of the denominator */ 00849 00850 template_lapack_laln2(&c_true, &na, &nw, &dmin__, &acoef, &a_ref(j, j), lda, 00851 bdiag, &bdiag[1], sum, &c__2, &bcoefr, &bcoefi, & 00852 work[(*n << 1) + j], n, &scale, &temp, &iinfo); 00853 if (scale < 1.) { 00854 i__3 = nw - 1; 00855 for (jw = 0; jw <= i__3; ++jw) { 00856 i__4 = j - 1; 00857 for (jr = je; jr <= i__4; ++jr) { 00858 work[(jw + 2) * *n + jr] = scale * work[(jw + 2) * 00859 *n + jr]; 00860 /* L140: */ 00861 } 00862 /* L150: */ 00863 } 00864 xmax = scale * xmax; 00865 } 00866 xmax = maxMACRO(xmax,temp); 00867 L160: 00868 ; 00869 } 00870 00871 /* Copy eigenvector to VL, back transforming if 00872 HOWMNY='B'. */ 00873 00874 ++ieig; 00875 if (ilback) { 00876 i__2 = nw - 1; 00877 for (jw = 0; jw <= i__2; ++jw) { 00878 i__3 = *n + 1 - je; 00879 template_blas_gemv("N", n, &i__3, &c_b35, &vl_ref(1, je), ldvl, &work[ 00880 (jw + 2) * *n + je], &c__1, &c_b37, &work[(jw + 4) 00881 * *n + 1], &c__1); 00882 /* L170: */ 00883 } 00884 template_lapack_lacpy(" ", n, &nw, &work[(*n << 2) + 1], n, &vl_ref(1, je), 00885 ldvl); 00886 ibeg = 1; 00887 } else { 00888 template_lapack_lacpy(" ", n, &nw, &work[(*n << 1) + 1], n, &vl_ref(1, ieig) 00889 , ldvl); 00890 ibeg = je; 00891 } 00892 00893 /* Scale eigenvector */ 00894 00895 xmax = 0.; 00896 if (ilcplx) { 00897 i__2 = *n; 00898 for (j = ibeg; j <= i__2; ++j) { 00899 /* Computing MAX */ 00900 d__3 = xmax, d__4 = (d__1 = vl_ref(j, ieig), absMACRO(d__1)) + 00901 (d__2 = vl_ref(j, ieig + 1), absMACRO(d__2)); 00902 xmax = maxMACRO(d__3,d__4); 00903 /* L180: */ 00904 } 00905 } else { 00906 i__2 = *n; 00907 for (j = ibeg; j <= i__2; ++j) { 00908 /* Computing MAX */ 00909 d__2 = xmax, d__3 = (d__1 = vl_ref(j, ieig), absMACRO(d__1)); 00910 xmax = maxMACRO(d__2,d__3); 00911 /* L190: */ 00912 } 00913 } 00914 00915 if (xmax > safmin) { 00916 xscale = 1. / xmax; 00917 00918 i__2 = nw - 1; 00919 for (jw = 0; jw <= i__2; ++jw) { 00920 i__3 = *n; 00921 for (jr = ibeg; jr <= i__3; ++jr) { 00922 vl_ref(jr, ieig + jw) = xscale * vl_ref(jr, ieig + jw) 00923 ; 00924 /* L200: */ 00925 } 00926 /* L210: */ 00927 } 00928 } 00929 ieig = ieig + nw - 1; 00930 00931 L220: 00932 ; 00933 } 00934 } 00935 00936 /* Right eigenvectors */ 00937 00938 if (compr) { 00939 ieig = im + 1; 00940 00941 /* Main loop over eigenvalues */ 00942 00943 ilcplx = FALSE_; 00944 for (je = *n; je >= 1; --je) { 00945 00946 /* Skip this iteration if (a) HOWMNY='S' and SELECT=.FALSE., or 00947 (b) this would be the second of a complex pair. 00948 Check for complex eigenvalue, so as to be sure of which 00949 entry(-ies) of SELECT to look at -- if complex, SELECT(JE) 00950 or SELECT(JE-1). 00951 If this is a complex pair, the 2-by-2 diagonal block 00952 corresponding to the eigenvalue is in rows/columns JE-1:JE */ 00953 00954 if (ilcplx) { 00955 ilcplx = FALSE_; 00956 goto L500; 00957 } 00958 nw = 1; 00959 if (je > 1) { 00960 if (a_ref(je, je - 1) != 0.) { 00961 ilcplx = TRUE_; 00962 nw = 2; 00963 } 00964 } 00965 if (ilall) { 00966 ilcomp = TRUE_; 00967 } else if (ilcplx) { 00968 ilcomp = select[je] || select[je - 1]; 00969 } else { 00970 ilcomp = select[je]; 00971 } 00972 if (! ilcomp) { 00973 goto L500; 00974 } 00975 00976 /* Decide if (a) singular pencil, (b) real eigenvalue, or 00977 (c) complex eigenvalue. */ 00978 00979 if (! ilcplx) { 00980 if ((d__1 = a_ref(je, je), absMACRO(d__1)) <= safmin && (d__2 = 00981 b_ref(je, je), absMACRO(d__2)) <= safmin) { 00982 00983 /* Singular matrix pencil -- unit eigenvector */ 00984 00985 --ieig; 00986 i__1 = *n; 00987 for (jr = 1; jr <= i__1; ++jr) { 00988 vr_ref(jr, ieig) = 0.; 00989 /* L230: */ 00990 } 00991 vr_ref(ieig, ieig) = 1.; 00992 goto L500; 00993 } 00994 } 00995 00996 /* Clear vector */ 00997 00998 i__1 = nw - 1; 00999 for (jw = 0; jw <= i__1; ++jw) { 01000 i__2 = *n; 01001 for (jr = 1; jr <= i__2; ++jr) { 01002 work[(jw + 2) * *n + jr] = 0.; 01003 /* L240: */ 01004 } 01005 /* L250: */ 01006 } 01007 01008 /* Compute coefficients in ( a A - b B ) x = 0 01009 a is ACOEF 01010 b is BCOEFR + i*BCOEFI */ 01011 01012 if (! ilcplx) { 01013 01014 /* Real eigenvalue 01015 01016 Computing MAX */ 01017 d__3 = (d__1 = a_ref(je, je), absMACRO(d__1)) * ascale, d__4 = ( 01018 d__2 = b_ref(je, je), absMACRO(d__2)) * bscale, d__3 = maxMACRO( 01019 d__3,d__4); 01020 temp = 1. / maxMACRO(d__3,safmin); 01021 salfar = temp * a_ref(je, je) * ascale; 01022 sbeta = temp * b_ref(je, je) * bscale; 01023 acoef = sbeta * ascale; 01024 bcoefr = salfar * bscale; 01025 bcoefi = 0.; 01026 01027 /* Scale to avoid underflow */ 01028 01029 scale = 1.; 01030 lsa = absMACRO(sbeta) >= safmin && absMACRO(acoef) < small; 01031 lsb = absMACRO(salfar) >= safmin && absMACRO(bcoefr) < small; 01032 if (lsa) { 01033 scale = small / absMACRO(sbeta) * minMACRO(anorm,big); 01034 } 01035 if (lsb) { 01036 /* Computing MAX */ 01037 d__1 = scale, d__2 = small / absMACRO(salfar) * minMACRO(bnorm,big); 01038 scale = maxMACRO(d__1,d__2); 01039 } 01040 if (lsa || lsb) { 01041 /* Computing MIN 01042 Computing MAX */ 01043 d__3 = 1., d__4 = absMACRO(acoef), d__3 = maxMACRO(d__3,d__4), d__4 01044 = absMACRO(bcoefr); 01045 d__1 = scale, d__2 = 1. / (safmin * maxMACRO(d__3,d__4)); 01046 scale = minMACRO(d__1,d__2); 01047 if (lsa) { 01048 acoef = ascale * (scale * sbeta); 01049 } else { 01050 acoef = scale * acoef; 01051 } 01052 if (lsb) { 01053 bcoefr = bscale * (scale * salfar); 01054 } else { 01055 bcoefr = scale * bcoefr; 01056 } 01057 } 01058 acoefa = absMACRO(acoef); 01059 bcoefa = absMACRO(bcoefr); 01060 01061 /* First component is 1 */ 01062 01063 work[(*n << 1) + je] = 1.; 01064 xmax = 1.; 01065 01066 /* Compute contribution from column JE of A and B to sum 01067 (See "Further Details", above.) */ 01068 01069 i__1 = je - 1; 01070 for (jr = 1; jr <= i__1; ++jr) { 01071 work[(*n << 1) + jr] = bcoefr * b_ref(jr, je) - acoef * 01072 a_ref(jr, je); 01073 /* L260: */ 01074 } 01075 } else { 01076 01077 /* Complex eigenvalue */ 01078 01079 d__1 = safmin * 100.; 01080 template_lapack_lag2(&a_ref(je - 1, je - 1), lda, &b_ref(je - 1, je - 1), 01081 ldb, &d__1, &acoef, &temp, &bcoefr, &temp2, &bcoefi); 01082 if (bcoefi == 0.) { 01083 *info = je - 1; 01084 return 0; 01085 } 01086 01087 /* Scale to avoid over/underflow */ 01088 01089 acoefa = absMACRO(acoef); 01090 bcoefa = absMACRO(bcoefr) + absMACRO(bcoefi); 01091 scale = 1.; 01092 if (acoefa * ulp < safmin && acoefa >= safmin) { 01093 scale = safmin / ulp / acoefa; 01094 } 01095 if (bcoefa * ulp < safmin && bcoefa >= safmin) { 01096 /* Computing MAX */ 01097 d__1 = scale, d__2 = safmin / ulp / bcoefa; 01098 scale = maxMACRO(d__1,d__2); 01099 } 01100 if (safmin * acoefa > ascale) { 01101 scale = ascale / (safmin * acoefa); 01102 } 01103 if (safmin * bcoefa > bscale) { 01104 /* Computing MIN */ 01105 d__1 = scale, d__2 = bscale / (safmin * bcoefa); 01106 scale = minMACRO(d__1,d__2); 01107 } 01108 if (scale != 1.) { 01109 acoef = scale * acoef; 01110 acoefa = absMACRO(acoef); 01111 bcoefr = scale * bcoefr; 01112 bcoefi = scale * bcoefi; 01113 bcoefa = absMACRO(bcoefr) + absMACRO(bcoefi); 01114 } 01115 01116 /* Compute first two components of eigenvector 01117 and contribution to sums */ 01118 01119 temp = acoef * a_ref(je, je - 1); 01120 temp2r = acoef * a_ref(je, je) - bcoefr * b_ref(je, je); 01121 temp2i = -bcoefi * b_ref(je, je); 01122 if (absMACRO(temp) >= absMACRO(temp2r) + absMACRO(temp2i)) { 01123 work[(*n << 1) + je] = 1.; 01124 work[*n * 3 + je] = 0.; 01125 work[(*n << 1) + je - 1] = -temp2r / temp; 01126 work[*n * 3 + je - 1] = -temp2i / temp; 01127 } else { 01128 work[(*n << 1) + je - 1] = 1.; 01129 work[*n * 3 + je - 1] = 0.; 01130 temp = acoef * a_ref(je - 1, je); 01131 work[(*n << 1) + je] = (bcoefr * b_ref(je - 1, je - 1) - 01132 acoef * a_ref(je - 1, je - 1)) / temp; 01133 work[*n * 3 + je] = bcoefi * b_ref(je - 1, je - 1) / temp; 01134 } 01135 01136 /* Computing MAX */ 01137 d__5 = (d__1 = work[(*n << 1) + je], absMACRO(d__1)) + (d__2 = 01138 work[*n * 3 + je], absMACRO(d__2)), d__6 = (d__3 = work[(* 01139 n << 1) + je - 1], absMACRO(d__3)) + (d__4 = work[*n * 3 + 01140 je - 1], absMACRO(d__4)); 01141 xmax = maxMACRO(d__5,d__6); 01142 01143 /* Compute contribution from columns JE and JE-1 01144 of A and B to the sums. */ 01145 01146 creala = acoef * work[(*n << 1) + je - 1]; 01147 cimaga = acoef * work[*n * 3 + je - 1]; 01148 crealb = bcoefr * work[(*n << 1) + je - 1] - bcoefi * work[*n 01149 * 3 + je - 1]; 01150 cimagb = bcoefi * work[(*n << 1) + je - 1] + bcoefr * work[*n 01151 * 3 + je - 1]; 01152 cre2a = acoef * work[(*n << 1) + je]; 01153 cim2a = acoef * work[*n * 3 + je]; 01154 cre2b = bcoefr * work[(*n << 1) + je] - bcoefi * work[*n * 3 01155 + je]; 01156 cim2b = bcoefi * work[(*n << 1) + je] + bcoefr * work[*n * 3 01157 + je]; 01158 i__1 = je - 2; 01159 for (jr = 1; jr <= i__1; ++jr) { 01160 work[(*n << 1) + jr] = -creala * a_ref(jr, je - 1) + 01161 crealb * b_ref(jr, je - 1) - cre2a * a_ref(jr, je) 01162 + cre2b * b_ref(jr, je); 01163 work[*n * 3 + jr] = -cimaga * a_ref(jr, je - 1) + cimagb * 01164 b_ref(jr, je - 1) - cim2a * a_ref(jr, je) + 01165 cim2b * b_ref(jr, je); 01166 /* L270: */ 01167 } 01168 } 01169 01170 /* Computing MAX */ 01171 d__1 = ulp * acoefa * anorm, d__2 = ulp * bcoefa * bnorm, d__1 = 01172 maxMACRO(d__1,d__2); 01173 dmin__ = maxMACRO(d__1,safmin); 01174 01175 /* Columnwise triangular solve of (a A - b B) x = 0 */ 01176 01177 il2by2 = FALSE_; 01178 for (j = je - nw; j >= 1; --j) { 01179 01180 /* If a 2-by-2 block, is in position j-1:j, wait until 01181 next iteration to process it (when it will be j:j+1) */ 01182 01183 if (! il2by2 && j > 1) { 01184 if (a_ref(j, j - 1) != 0.) { 01185 il2by2 = TRUE_; 01186 goto L370; 01187 } 01188 } 01189 bdiag[0] = b_ref(j, j); 01190 if (il2by2) { 01191 na = 2; 01192 bdiag[1] = b_ref(j + 1, j + 1); 01193 } else { 01194 na = 1; 01195 } 01196 01197 /* Compute x(j) (and x(j+1), if 2-by-2 block) */ 01198 01199 template_lapack_laln2(&c_false, &na, &nw, &dmin__, &acoef, &a_ref(j, j), 01200 lda, bdiag, &bdiag[1], &work[(*n << 1) + j], n, & 01201 bcoefr, &bcoefi, sum, &c__2, &scale, &temp, &iinfo); 01202 if (scale < 1.) { 01203 01204 i__1 = nw - 1; 01205 for (jw = 0; jw <= i__1; ++jw) { 01206 i__2 = je; 01207 for (jr = 1; jr <= i__2; ++jr) { 01208 work[(jw + 2) * *n + jr] = scale * work[(jw + 2) * 01209 *n + jr]; 01210 /* L280: */ 01211 } 01212 /* L290: */ 01213 } 01214 } 01215 /* Computing MAX */ 01216 d__1 = scale * xmax; 01217 xmax = maxMACRO(d__1,temp); 01218 01219 i__1 = nw; 01220 for (jw = 1; jw <= i__1; ++jw) { 01221 i__2 = na; 01222 for (ja = 1; ja <= i__2; ++ja) { 01223 work[(jw + 1) * *n + j + ja - 1] = sum_ref(ja, jw); 01224 /* L300: */ 01225 } 01226 /* L310: */ 01227 } 01228 01229 /* w = w + x(j)*(a A(*,j) - b B(*,j) ) with scaling */ 01230 01231 if (j > 1) { 01232 01233 /* Check whether scaling is necessary for sum. */ 01234 01235 xscale = 1. / maxMACRO(1.,xmax); 01236 temp = acoefa * work[j] + bcoefa * work[*n + j]; 01237 if (il2by2) { 01238 /* Computing MAX */ 01239 d__1 = temp, d__2 = acoefa * work[j + 1] + bcoefa * 01240 work[*n + j + 1]; 01241 temp = maxMACRO(d__1,d__2); 01242 } 01243 /* Computing MAX */ 01244 d__1 = maxMACRO(temp,acoefa); 01245 temp = maxMACRO(d__1,bcoefa); 01246 if (temp > bignum * xscale) { 01247 01248 i__1 = nw - 1; 01249 for (jw = 0; jw <= i__1; ++jw) { 01250 i__2 = je; 01251 for (jr = 1; jr <= i__2; ++jr) { 01252 work[(jw + 2) * *n + jr] = xscale * work[(jw 01253 + 2) * *n + jr]; 01254 /* L320: */ 01255 } 01256 /* L330: */ 01257 } 01258 xmax *= xscale; 01259 } 01260 01261 /* Compute the contributions of the off-diagonals of 01262 column j (and j+1, if 2-by-2 block) of A and B to the 01263 sums. */ 01264 01265 01266 i__1 = na; 01267 for (ja = 1; ja <= i__1; ++ja) { 01268 if (ilcplx) { 01269 creala = acoef * work[(*n << 1) + j + ja - 1]; 01270 cimaga = acoef * work[*n * 3 + j + ja - 1]; 01271 crealb = bcoefr * work[(*n << 1) + j + ja - 1] - 01272 bcoefi * work[*n * 3 + j + ja - 1]; 01273 cimagb = bcoefi * work[(*n << 1) + j + ja - 1] + 01274 bcoefr * work[*n * 3 + j + ja - 1]; 01275 i__2 = j - 1; 01276 for (jr = 1; jr <= i__2; ++jr) { 01277 work[(*n << 1) + jr] = work[(*n << 1) + jr] - 01278 creala * a_ref(jr, j + ja - 1) + 01279 crealb * b_ref(jr, j + ja - 1); 01280 work[*n * 3 + jr] = work[*n * 3 + jr] - 01281 cimaga * a_ref(jr, j + ja - 1) + 01282 cimagb * b_ref(jr, j + ja - 1); 01283 /* L340: */ 01284 } 01285 } else { 01286 creala = acoef * work[(*n << 1) + j + ja - 1]; 01287 crealb = bcoefr * work[(*n << 1) + j + ja - 1]; 01288 i__2 = j - 1; 01289 for (jr = 1; jr <= i__2; ++jr) { 01290 work[(*n << 1) + jr] = work[(*n << 1) + jr] - 01291 creala * a_ref(jr, j + ja - 1) + 01292 crealb * b_ref(jr, j + ja - 1); 01293 /* L350: */ 01294 } 01295 } 01296 /* L360: */ 01297 } 01298 } 01299 01300 il2by2 = FALSE_; 01301 L370: 01302 ; 01303 } 01304 01305 /* Copy eigenvector to VR, back transforming if 01306 HOWMNY='B'. */ 01307 01308 ieig -= nw; 01309 if (ilback) { 01310 01311 i__1 = nw - 1; 01312 for (jw = 0; jw <= i__1; ++jw) { 01313 i__2 = *n; 01314 for (jr = 1; jr <= i__2; ++jr) { 01315 work[(jw + 4) * *n + jr] = work[(jw + 2) * *n + 1] * 01316 vr_ref(jr, 1); 01317 /* L380: */ 01318 } 01319 01320 /* A series of compiler directives to defeat 01321 vectorization for the next loop */ 01322 01323 01324 i__2 = je; 01325 for (jc = 2; jc <= i__2; ++jc) { 01326 i__3 = *n; 01327 for (jr = 1; jr <= i__3; ++jr) { 01328 work[(jw + 4) * *n + jr] += work[(jw + 2) * *n + 01329 jc] * vr_ref(jr, jc); 01330 /* L390: */ 01331 } 01332 /* L400: */ 01333 } 01334 /* L410: */ 01335 } 01336 01337 i__1 = nw - 1; 01338 for (jw = 0; jw <= i__1; ++jw) { 01339 i__2 = *n; 01340 for (jr = 1; jr <= i__2; ++jr) { 01341 vr_ref(jr, ieig + jw) = work[(jw + 4) * *n + jr]; 01342 /* L420: */ 01343 } 01344 /* L430: */ 01345 } 01346 01347 iend = *n; 01348 } else { 01349 i__1 = nw - 1; 01350 for (jw = 0; jw <= i__1; ++jw) { 01351 i__2 = *n; 01352 for (jr = 1; jr <= i__2; ++jr) { 01353 vr_ref(jr, ieig + jw) = work[(jw + 2) * *n + jr]; 01354 /* L440: */ 01355 } 01356 /* L450: */ 01357 } 01358 01359 iend = je; 01360 } 01361 01362 /* Scale eigenvector */ 01363 01364 xmax = 0.; 01365 if (ilcplx) { 01366 i__1 = iend; 01367 for (j = 1; j <= i__1; ++j) { 01368 /* Computing MAX */ 01369 d__3 = xmax, d__4 = (d__1 = vr_ref(j, ieig), absMACRO(d__1)) + 01370 (d__2 = vr_ref(j, ieig + 1), absMACRO(d__2)); 01371 xmax = maxMACRO(d__3,d__4); 01372 /* L460: */ 01373 } 01374 } else { 01375 i__1 = iend; 01376 for (j = 1; j <= i__1; ++j) { 01377 /* Computing MAX */ 01378 d__2 = xmax, d__3 = (d__1 = vr_ref(j, ieig), absMACRO(d__1)); 01379 xmax = maxMACRO(d__2,d__3); 01380 /* L470: */ 01381 } 01382 } 01383 01384 if (xmax > safmin) { 01385 xscale = 1. / xmax; 01386 i__1 = nw - 1; 01387 for (jw = 0; jw <= i__1; ++jw) { 01388 i__2 = iend; 01389 for (jr = 1; jr <= i__2; ++jr) { 01390 vr_ref(jr, ieig + jw) = xscale * vr_ref(jr, ieig + jw) 01391 ; 01392 /* L480: */ 01393 } 01394 /* L490: */ 01395 } 01396 } 01397 L500: 01398 ; 01399 } 01400 } 01401 01402 return 0; 01403 01404 /* End of DTGEVC */ 01405 01406 } /* dtgevc_ */ 01407 01408 #undef sum_ref 01409 #undef vr_ref 01410 #undef vl_ref 01411 #undef b_ref 01412 #undef a_ref 01413 #undef sumb_ref 01414 #undef suma_ref 01415 01416 01417 #endif