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_HGEQZ_HEADER 00036 #define TEMPLATE_LAPACK_HGEQZ_HEADER 00037 00038 00039 template<class Treal> 00040 int template_lapack_hgeqz(const char *job, const char *compq, const char *compz, const integer *n, 00041 const integer *ilo, const integer *ihi, Treal *a, const integer *lda, Treal * 00042 b, const integer *ldb, Treal *alphar, Treal *alphai, Treal * 00043 beta, Treal *q, const integer *ldq, Treal *z__, const integer *ldz, 00044 Treal *work, const integer *lwork, integer *info) 00045 { 00046 /* -- LAPACK routine (version 3.0) -- 00047 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., 00048 Courant Institute, Argonne National Lab, and Rice University 00049 June 30, 1999 00050 00051 00052 Purpose 00053 ======= 00054 00055 DHGEQZ implements a single-/double-shift version of the QZ method for 00056 finding the generalized eigenvalues 00057 00058 w(j)=(ALPHAR(j) + i*ALPHAI(j))/BETAR(j) of the equation 00059 00060 det( A - w(i) B ) = 0 00061 00062 In addition, the pair A,B may be reduced to generalized Schur form: 00063 B is upper triangular, and A is block upper triangular, where the 00064 diagonal blocks are either 1-by-1 or 2-by-2, the 2-by-2 blocks having 00065 complex generalized eigenvalues (see the description of the argument 00066 JOB.) 00067 00068 If JOB='S', then the pair (A,B) is simultaneously reduced to Schur 00069 form by applying one orthogonal tranformation (usually called Q) on 00070 the left and another (usually called Z) on the right. The 2-by-2 00071 upper-triangular diagonal blocks of B corresponding to 2-by-2 blocks 00072 of A will be reduced to positive diagonal matrices. (I.e., 00073 if A(j+1,j) is non-zero, then B(j+1,j)=B(j,j+1)=0 and B(j,j) and 00074 B(j+1,j+1) will be positive.) 00075 00076 If JOB='E', then at each iteration, the same transformations 00077 are computed, but they are only applied to those parts of A and B 00078 which are needed to compute ALPHAR, ALPHAI, and BETAR. 00079 00080 If JOB='S' and COMPQ and COMPZ are 'V' or 'I', then the orthogonal 00081 transformations used to reduce (A,B) are accumulated into the arrays 00082 Q and Z s.t.: 00083 00084 Q(in) A(in) Z(in)* = Q(out) A(out) Z(out)* 00085 Q(in) B(in) Z(in)* = Q(out) B(out) Z(out)* 00086 00087 Ref: C.B. Moler & G.W. Stewart, "An Algorithm for Generalized Matrix 00088 Eigenvalue Problems", SIAM J. Numer. Anal., 10(1973), 00089 pp. 241--256. 00090 00091 Arguments 00092 ========= 00093 00094 JOB (input) CHARACTER*1 00095 = 'E': compute only ALPHAR, ALPHAI, and BETA. A and B will 00096 not necessarily be put into generalized Schur form. 00097 = 'S': put A and B into generalized Schur form, as well 00098 as computing ALPHAR, ALPHAI, and BETA. 00099 00100 COMPQ (input) CHARACTER*1 00101 = 'N': do not modify Q. 00102 = 'V': multiply the array Q on the right by the transpose of 00103 the orthogonal tranformation that is applied to the 00104 left side of A and B to reduce them to Schur form. 00105 = 'I': like COMPQ='V', except that Q will be initialized to 00106 the identity first. 00107 00108 COMPZ (input) CHARACTER*1 00109 = 'N': do not modify Z. 00110 = 'V': multiply the array Z on the right by the orthogonal 00111 tranformation that is applied to the right side of 00112 A and B to reduce them to Schur form. 00113 = 'I': like COMPZ='V', except that Z will be initialized to 00114 the identity first. 00115 00116 N (input) INTEGER 00117 The order of the matrices A, B, Q, and Z. N >= 0. 00118 00119 ILO (input) INTEGER 00120 IHI (input) INTEGER 00121 It is assumed that A is already upper triangular in rows and 00122 columns 1:ILO-1 and IHI+1:N. 00123 1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0. 00124 00125 A (input/output) DOUBLE PRECISION array, dimension (LDA, N) 00126 On entry, the N-by-N upper Hessenberg matrix A. Elements 00127 below the subdiagonal must be zero. 00128 If JOB='S', then on exit A and B will have been 00129 simultaneously reduced to generalized Schur form. 00130 If JOB='E', then on exit A will have been destroyed. 00131 The diagonal blocks will be correct, but the off-diagonal 00132 portion will be meaningless. 00133 00134 LDA (input) INTEGER 00135 The leading dimension of the array A. LDA >= max( 1, N ). 00136 00137 B (input/output) DOUBLE PRECISION array, dimension (LDB, N) 00138 On entry, the N-by-N upper triangular matrix B. Elements 00139 below the diagonal must be zero. 2-by-2 blocks in B 00140 corresponding to 2-by-2 blocks in A will be reduced to 00141 positive diagonal form. (I.e., if A(j+1,j) is non-zero, 00142 then B(j+1,j)=B(j,j+1)=0 and B(j,j) and B(j+1,j+1) will be 00143 positive.) 00144 If JOB='S', then on exit A and B will have been 00145 simultaneously reduced to Schur form. 00146 If JOB='E', then on exit B will have been destroyed. 00147 Elements corresponding to diagonal blocks of A will be 00148 correct, but the off-diagonal portion will be meaningless. 00149 00150 LDB (input) INTEGER 00151 The leading dimension of the array B. LDB >= max( 1, N ). 00152 00153 ALPHAR (output) DOUBLE PRECISION array, dimension (N) 00154 ALPHAR(1:N) will be set to real parts of the diagonal 00155 elements of A that would result from reducing A and B to 00156 Schur form and then further reducing them both to triangular 00157 form using unitary transformations s.t. the diagonal of B 00158 was non-negative real. Thus, if A(j,j) is in a 1-by-1 block 00159 (i.e., A(j+1,j)=A(j,j+1)=0), then ALPHAR(j)=A(j,j). 00160 Note that the (real or complex) values 00161 (ALPHAR(j) + i*ALPHAI(j))/BETA(j), j=1,...,N, are the 00162 generalized eigenvalues of the matrix pencil A - wB. 00163 00164 ALPHAI (output) DOUBLE PRECISION array, dimension (N) 00165 ALPHAI(1:N) will be set to imaginary parts of the diagonal 00166 elements of A that would result from reducing A and B to 00167 Schur form and then further reducing them both to triangular 00168 form using unitary transformations s.t. the diagonal of B 00169 was non-negative real. Thus, if A(j,j) is in a 1-by-1 block 00170 (i.e., A(j+1,j)=A(j,j+1)=0), then ALPHAR(j)=0. 00171 Note that the (real or complex) values 00172 (ALPHAR(j) + i*ALPHAI(j))/BETA(j), j=1,...,N, are the 00173 generalized eigenvalues of the matrix pencil A - wB. 00174 00175 BETA (output) DOUBLE PRECISION array, dimension (N) 00176 BETA(1:N) will be set to the (real) diagonal elements of B 00177 that would result from reducing A and B to Schur form and 00178 then further reducing them both to triangular form using 00179 unitary transformations s.t. the diagonal of B was 00180 non-negative real. Thus, if A(j,j) is in a 1-by-1 block 00181 (i.e., A(j+1,j)=A(j,j+1)=0), then BETA(j)=B(j,j). 00182 Note that the (real or complex) values 00183 (ALPHAR(j) + i*ALPHAI(j))/BETA(j), j=1,...,N, are the 00184 generalized eigenvalues of the matrix pencil A - wB. 00185 (Note that BETA(1:N) will always be non-negative, and no 00186 BETAI is necessary.) 00187 00188 Q (input/output) DOUBLE PRECISION array, dimension (LDQ, N) 00189 If COMPQ='N', then Q will not be referenced. 00190 If COMPQ='V' or 'I', then the transpose of the orthogonal 00191 transformations which are applied to A and B on the left 00192 will be applied to the array Q on the right. 00193 00194 LDQ (input) INTEGER 00195 The leading dimension of the array Q. LDQ >= 1. 00196 If COMPQ='V' or 'I', then LDQ >= N. 00197 00198 Z (input/output) DOUBLE PRECISION array, dimension (LDZ, N) 00199 If COMPZ='N', then Z will not be referenced. 00200 If COMPZ='V' or 'I', then the orthogonal transformations 00201 which are applied to A and B on the right will be applied 00202 to the array Z on the right. 00203 00204 LDZ (input) INTEGER 00205 The leading dimension of the array Z. LDZ >= 1. 00206 If COMPZ='V' or 'I', then LDZ >= N. 00207 00208 WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK) 00209 On exit, if INFO >= 0, WORK(1) returns the optimal LWORK. 00210 00211 LWORK (input) INTEGER 00212 The dimension of the array WORK. LWORK >= max(1,N). 00213 00214 If LWORK = -1, then a workspace query is assumed; the routine 00215 only calculates the optimal size of the WORK array, returns 00216 this value as the first entry of the WORK array, and no error 00217 message related to LWORK is issued by XERBLA. 00218 00219 INFO (output) INTEGER 00220 = 0: successful exit 00221 < 0: if INFO = -i, the i-th argument had an illegal value 00222 = 1,...,N: the QZ iteration did not converge. (A,B) is not 00223 in Schur form, but ALPHAR(i), ALPHAI(i), and 00224 BETA(i), i=INFO+1,...,N should be correct. 00225 = N+1,...,2*N: the shift calculation failed. (A,B) is not 00226 in Schur form, but ALPHAR(i), ALPHAI(i), and 00227 BETA(i), i=INFO-N+1,...,N should be correct. 00228 > 2*N: various "impossible" errors. 00229 00230 Further Details 00231 =============== 00232 00233 Iteration counters: 00234 00235 JITER -- counts iterations. 00236 IITER -- counts iterations run since ILAST was last 00237 changed. This is therefore reset only when a 1-by-1 or 00238 2-by-2 block deflates off the bottom. 00239 00240 ===================================================================== 00241 00242 $ SAFETY = 1.0E+0 ) 00243 00244 Decode JOB, COMPQ, COMPZ 00245 00246 Parameter adjustments */ 00247 /* Table of constant values */ 00248 Treal c_b12 = 0.; 00249 Treal c_b13 = 1.; 00250 integer c__1 = 1; 00251 integer c__3 = 3; 00252 00253 /* System generated locals */ 00254 integer a_dim1, a_offset, b_dim1, b_offset, q_dim1, q_offset, z_dim1, 00255 z_offset, i__1, i__2, i__3, i__4; 00256 Treal d__1, d__2, d__3, d__4; 00257 /* Local variables */ 00258 Treal ad11l, ad12l, ad21l, ad22l, ad32l, wabs, atol, btol, 00259 temp; 00260 Treal temp2, s1inv, c__; 00261 integer j; 00262 Treal s, t, v[3], scale; 00263 integer iiter, ilast, jiter; 00264 Treal anorm, bnorm; 00265 integer maxit; 00266 Treal tempi, tempr, s1, s2, u1, u2; 00267 logical ilazr2; 00268 Treal a11, a12, a21, a22, b11, b22, c12, c21; 00269 integer jc; 00270 Treal an, bn, cl, cq, cr; 00271 integer in; 00272 Treal ascale, bscale, u12, w11; 00273 integer jr; 00274 Treal cz, sl, w12, w21, w22, wi; 00275 Treal sr; 00276 Treal vs, wr; 00277 Treal safmin; 00278 Treal safmax; 00279 Treal eshift; 00280 logical ilschr; 00281 Treal b1a, b2a; 00282 integer icompq, ilastm; 00283 Treal a1i; 00284 integer ischur; 00285 Treal a2i, b1i; 00286 logical ilazro; 00287 integer icompz, ifirst; 00288 Treal b2i; 00289 integer ifrstm; 00290 Treal a1r; 00291 integer istart; 00292 logical ilpivt; 00293 Treal a2r, b1r, b2r; 00294 logical lquery; 00295 Treal wr2, ad11, ad12, ad21, ad22, c11i, c22i; 00296 integer jch; 00297 Treal c11r, c22r, u12l; 00298 logical ilq; 00299 Treal tau, sqi; 00300 logical ilz; 00301 Treal ulp, sqr, szi, szr; 00302 #define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1] 00303 #define b_ref(a_1,a_2) b[(a_2)*b_dim1 + a_1] 00304 #define q_ref(a_1,a_2) q[(a_2)*q_dim1 + a_1] 00305 #define z___ref(a_1,a_2) z__[(a_2)*z_dim1 + a_1] 00306 00307 00308 a_dim1 = *lda; 00309 a_offset = 1 + a_dim1 * 1; 00310 a -= a_offset; 00311 b_dim1 = *ldb; 00312 b_offset = 1 + b_dim1 * 1; 00313 b -= b_offset; 00314 --alphar; 00315 --alphai; 00316 --beta; 00317 q_dim1 = *ldq; 00318 q_offset = 1 + q_dim1 * 1; 00319 q -= q_offset; 00320 z_dim1 = *ldz; 00321 z_offset = 1 + z_dim1 * 1; 00322 z__ -= z_offset; 00323 --work; 00324 00325 /* Initialization added by Elias to get rid of compiler warnings. */ 00326 ilschr = ilq = ilz = 0; 00327 /* Function Body */ 00328 if (template_blas_lsame(job, "E")) { 00329 ilschr = FALSE_; 00330 ischur = 1; 00331 } else if (template_blas_lsame(job, "S")) { 00332 ilschr = TRUE_; 00333 ischur = 2; 00334 } else { 00335 ischur = 0; 00336 } 00337 00338 if (template_blas_lsame(compq, "N")) { 00339 ilq = FALSE_; 00340 icompq = 1; 00341 } else if (template_blas_lsame(compq, "V")) { 00342 ilq = TRUE_; 00343 icompq = 2; 00344 } else if (template_blas_lsame(compq, "I")) { 00345 ilq = TRUE_; 00346 icompq = 3; 00347 } else { 00348 icompq = 0; 00349 } 00350 00351 if (template_blas_lsame(compz, "N")) { 00352 ilz = FALSE_; 00353 icompz = 1; 00354 } else if (template_blas_lsame(compz, "V")) { 00355 ilz = TRUE_; 00356 icompz = 2; 00357 } else if (template_blas_lsame(compz, "I")) { 00358 ilz = TRUE_; 00359 icompz = 3; 00360 } else { 00361 icompz = 0; 00362 } 00363 00364 /* Check Argument Values */ 00365 00366 *info = 0; 00367 work[1] = (Treal) maxMACRO(1,*n); 00368 lquery = *lwork == -1; 00369 if (ischur == 0) { 00370 *info = -1; 00371 } else if (icompq == 0) { 00372 *info = -2; 00373 } else if (icompz == 0) { 00374 *info = -3; 00375 } else if (*n < 0) { 00376 *info = -4; 00377 } else if (*ilo < 1) { 00378 *info = -5; 00379 } else if (*ihi > *n || *ihi < *ilo - 1) { 00380 *info = -6; 00381 } else if (*lda < *n) { 00382 *info = -8; 00383 } else if (*ldb < *n) { 00384 *info = -10; 00385 } else if (*ldq < 1 || ( ilq && *ldq < *n ) ) { 00386 *info = -15; 00387 } else if (*ldz < 1 || ( ilz && *ldz < *n ) ) { 00388 *info = -17; 00389 } else if (*lwork < maxMACRO(1,*n) && ! lquery) { 00390 *info = -19; 00391 } 00392 if (*info != 0) { 00393 i__1 = -(*info); 00394 template_blas_erbla("HGEQZ ", &i__1); 00395 return 0; 00396 } else if (lquery) { 00397 return 0; 00398 } 00399 00400 /* Quick return if possible */ 00401 00402 if (*n <= 0) { 00403 work[1] = 1.; 00404 return 0; 00405 } 00406 00407 /* Initialize Q and Z */ 00408 00409 if (icompq == 3) { 00410 template_lapack_laset("Full", n, n, &c_b12, &c_b13, &q[q_offset], ldq); 00411 } 00412 if (icompz == 3) { 00413 template_lapack_laset("Full", n, n, &c_b12, &c_b13, &z__[z_offset], ldz); 00414 } 00415 00416 /* Machine Constants */ 00417 00418 in = *ihi + 1 - *ilo; 00419 safmin = template_lapack_lamch("S", (Treal)0); 00420 safmax = 1. / safmin; 00421 ulp = template_lapack_lamch("E", (Treal)0) * template_lapack_lamch("B", (Treal)0); 00422 anorm = dlanhs_("F", &in, &a_ref(*ilo, *ilo), lda, &work[1]); 00423 bnorm = dlanhs_("F", &in, &b_ref(*ilo, *ilo), ldb, &work[1]); 00424 /* Computing MAX */ 00425 d__1 = safmin, d__2 = ulp * anorm; 00426 atol = maxMACRO(d__1,d__2); 00427 /* Computing MAX */ 00428 d__1 = safmin, d__2 = ulp * bnorm; 00429 btol = maxMACRO(d__1,d__2); 00430 ascale = 1. / maxMACRO(safmin,anorm); 00431 bscale = 1. / maxMACRO(safmin,bnorm); 00432 00433 /* Set Eigenvalues IHI+1:N */ 00434 00435 i__1 = *n; 00436 for (j = *ihi + 1; j <= i__1; ++j) { 00437 if (b_ref(j, j) < 0.) { 00438 if (ilschr) { 00439 i__2 = j; 00440 for (jr = 1; jr <= i__2; ++jr) { 00441 a_ref(jr, j) = -a_ref(jr, j); 00442 b_ref(jr, j) = -b_ref(jr, j); 00443 /* L10: */ 00444 } 00445 } else { 00446 a_ref(j, j) = -a_ref(j, j); 00447 b_ref(j, j) = -b_ref(j, j); 00448 } 00449 if (ilz) { 00450 i__2 = *n; 00451 for (jr = 1; jr <= i__2; ++jr) { 00452 z___ref(jr, j) = -z___ref(jr, j); 00453 /* L20: */ 00454 } 00455 } 00456 } 00457 alphar[j] = a_ref(j, j); 00458 alphai[j] = 0.; 00459 beta[j] = b_ref(j, j); 00460 /* L30: */ 00461 } 00462 00463 /* If IHI < ILO, skip QZ steps */ 00464 00465 if (*ihi < *ilo) { 00466 goto L380; 00467 } 00468 00469 /* MAIN QZ ITERATION LOOP 00470 00471 Initialize dynamic indices 00472 00473 Eigenvalues ILAST+1:N have been found. 00474 Column operations modify rows IFRSTM:whatever. 00475 Row operations modify columns whatever:ILASTM. 00476 00477 If only eigenvalues are being computed, then 00478 IFRSTM is the row of the last splitting row above row ILAST; 00479 this is always at least ILO. 00480 IITER counts iterations since the last eigenvalue was found, 00481 to tell when to use an extraordinary shift. 00482 MAXIT is the maximum number of QZ sweeps allowed. */ 00483 00484 ilast = *ihi; 00485 if (ilschr) { 00486 ifrstm = 1; 00487 ilastm = *n; 00488 } else { 00489 ifrstm = *ilo; 00490 ilastm = *ihi; 00491 } 00492 iiter = 0; 00493 eshift = 0.; 00494 maxit = (*ihi - *ilo + 1) * 30; 00495 00496 i__1 = maxit; 00497 for (jiter = 1; jiter <= i__1; ++jiter) { 00498 00499 /* Split the matrix if possible. 00500 00501 Two tests: 00502 1: A(j,j-1)=0 or j=ILO 00503 2: B(j,j)=0 */ 00504 00505 if (ilast == *ilo) { 00506 00507 /* Special case: j=ILAST */ 00508 00509 goto L80; 00510 } else { 00511 if ((d__1 = a_ref(ilast, ilast - 1), absMACRO(d__1)) <= atol) { 00512 a_ref(ilast, ilast - 1) = 0.; 00513 goto L80; 00514 } 00515 } 00516 00517 if ((d__1 = b_ref(ilast, ilast), absMACRO(d__1)) <= btol) { 00518 b_ref(ilast, ilast) = 0.; 00519 goto L70; 00520 } 00521 00522 /* General case: j<ILAST */ 00523 00524 i__2 = *ilo; 00525 for (j = ilast - 1; j >= i__2; --j) { 00526 00527 /* Test 1: for A(j,j-1)=0 or j=ILO */ 00528 00529 if (j == *ilo) { 00530 ilazro = TRUE_; 00531 } else { 00532 if ((d__1 = a_ref(j, j - 1), absMACRO(d__1)) <= atol) { 00533 a_ref(j, j - 1) = 0.; 00534 ilazro = TRUE_; 00535 } else { 00536 ilazro = FALSE_; 00537 } 00538 } 00539 00540 /* Test 2: for B(j,j)=0 */ 00541 00542 if ((d__1 = b_ref(j, j), absMACRO(d__1)) < btol) { 00543 b_ref(j, j) = 0.; 00544 00545 /* Test 1a: Check for 2 consecutive small subdiagonals in A */ 00546 00547 ilazr2 = FALSE_; 00548 if (! ilazro) { 00549 temp = (d__1 = a_ref(j, j - 1), absMACRO(d__1)); 00550 temp2 = (d__1 = a_ref(j, j), absMACRO(d__1)); 00551 tempr = maxMACRO(temp,temp2); 00552 if (tempr < 1. && tempr != 0.) { 00553 temp /= tempr; 00554 temp2 /= tempr; 00555 } 00556 if (temp * (ascale * (d__1 = a_ref(j + 1, j), absMACRO(d__1))) 00557 <= temp2 * (ascale * atol)) { 00558 ilazr2 = TRUE_; 00559 } 00560 } 00561 00562 /* If both tests pass (1 & 2), i.e., the leading diagonal 00563 element of B in the block is zero, split a 1x1 block off 00564 at the top. (I.e., at the J-th row/column) The leading 00565 diagonal element of the remainder can also be zero, so 00566 this may have to be done repeatedly. */ 00567 00568 if (ilazro || ilazr2) { 00569 i__3 = ilast - 1; 00570 for (jch = j; jch <= i__3; ++jch) { 00571 temp = a_ref(jch, jch); 00572 template_lapack_lartg(&temp, &a_ref(jch + 1, jch), &c__, &s, &a_ref( 00573 jch, jch)); 00574 a_ref(jch + 1, jch) = 0.; 00575 i__4 = ilastm - jch; 00576 template_blas_rot(&i__4, &a_ref(jch, jch + 1), lda, &a_ref(jch + 00577 1, jch + 1), lda, &c__, &s); 00578 i__4 = ilastm - jch; 00579 template_blas_rot(&i__4, &b_ref(jch, jch + 1), ldb, &b_ref(jch + 00580 1, jch + 1), ldb, &c__, &s); 00581 if (ilq) { 00582 template_blas_rot(n, &q_ref(1, jch), &c__1, &q_ref(1, jch + 1) 00583 , &c__1, &c__, &s); 00584 } 00585 if (ilazr2) { 00586 a_ref(jch, jch - 1) = a_ref(jch, jch - 1) * c__; 00587 } 00588 ilazr2 = FALSE_; 00589 if ((d__1 = b_ref(jch + 1, jch + 1), absMACRO(d__1)) >= 00590 btol) { 00591 if (jch + 1 >= ilast) { 00592 goto L80; 00593 } else { 00594 ifirst = jch + 1; 00595 goto L110; 00596 } 00597 } 00598 b_ref(jch + 1, jch + 1) = 0.; 00599 /* L40: */ 00600 } 00601 goto L70; 00602 } else { 00603 00604 /* Only test 2 passed -- chase the zero to B(ILAST,ILAST) 00605 Then process as in the case B(ILAST,ILAST)=0 */ 00606 00607 i__3 = ilast - 1; 00608 for (jch = j; jch <= i__3; ++jch) { 00609 temp = b_ref(jch, jch + 1); 00610 template_lapack_lartg(&temp, &b_ref(jch + 1, jch + 1), &c__, &s, & 00611 b_ref(jch, jch + 1)); 00612 b_ref(jch + 1, jch + 1) = 0.; 00613 if (jch < ilastm - 1) { 00614 i__4 = ilastm - jch - 1; 00615 template_blas_rot(&i__4, &b_ref(jch, jch + 2), ldb, &b_ref( 00616 jch + 1, jch + 2), ldb, &c__, &s); 00617 } 00618 i__4 = ilastm - jch + 2; 00619 template_blas_rot(&i__4, &a_ref(jch, jch - 1), lda, &a_ref(jch + 00620 1, jch - 1), lda, &c__, &s); 00621 if (ilq) { 00622 template_blas_rot(n, &q_ref(1, jch), &c__1, &q_ref(1, jch + 1) 00623 , &c__1, &c__, &s); 00624 } 00625 temp = a_ref(jch + 1, jch); 00626 template_lapack_lartg(&temp, &a_ref(jch + 1, jch - 1), &c__, &s, & 00627 a_ref(jch + 1, jch)); 00628 a_ref(jch + 1, jch - 1) = 0.; 00629 i__4 = jch + 1 - ifrstm; 00630 template_blas_rot(&i__4, &a_ref(ifrstm, jch), &c__1, &a_ref( 00631 ifrstm, jch - 1), &c__1, &c__, &s); 00632 i__4 = jch - ifrstm; 00633 template_blas_rot(&i__4, &b_ref(ifrstm, jch), &c__1, &b_ref( 00634 ifrstm, jch - 1), &c__1, &c__, &s); 00635 if (ilz) { 00636 template_blas_rot(n, &z___ref(1, jch), &c__1, &z___ref(1, jch 00637 - 1), &c__1, &c__, &s); 00638 } 00639 /* L50: */ 00640 } 00641 goto L70; 00642 } 00643 } else if (ilazro) { 00644 00645 /* Only test 1 passed -- work on J:ILAST */ 00646 00647 ifirst = j; 00648 goto L110; 00649 } 00650 00651 /* Neither test passed -- try next J 00652 00653 L60: */ 00654 } 00655 00656 /* (Drop-through is "impossible") */ 00657 00658 *info = *n + 1; 00659 goto L420; 00660 00661 /* B(ILAST,ILAST)=0 -- clear A(ILAST,ILAST-1) to split off a 00662 1x1 block. */ 00663 00664 L70: 00665 temp = a_ref(ilast, ilast); 00666 template_lapack_lartg(&temp, &a_ref(ilast, ilast - 1), &c__, &s, &a_ref(ilast, 00667 ilast)); 00668 a_ref(ilast, ilast - 1) = 0.; 00669 i__2 = ilast - ifrstm; 00670 template_blas_rot(&i__2, &a_ref(ifrstm, ilast), &c__1, &a_ref(ifrstm, ilast - 1), 00671 &c__1, &c__, &s); 00672 i__2 = ilast - ifrstm; 00673 template_blas_rot(&i__2, &b_ref(ifrstm, ilast), &c__1, &b_ref(ifrstm, ilast - 1), 00674 &c__1, &c__, &s); 00675 if (ilz) { 00676 template_blas_rot(n, &z___ref(1, ilast), &c__1, &z___ref(1, ilast - 1), &c__1, 00677 &c__, &s); 00678 } 00679 00680 /* A(ILAST,ILAST-1)=0 -- Standardize B, set ALPHAR, ALPHAI, 00681 and BETA */ 00682 00683 L80: 00684 if (b_ref(ilast, ilast) < 0.) { 00685 if (ilschr) { 00686 i__2 = ilast; 00687 for (j = ifrstm; j <= i__2; ++j) { 00688 a_ref(j, ilast) = -a_ref(j, ilast); 00689 b_ref(j, ilast) = -b_ref(j, ilast); 00690 /* L90: */ 00691 } 00692 } else { 00693 a_ref(ilast, ilast) = -a_ref(ilast, ilast); 00694 b_ref(ilast, ilast) = -b_ref(ilast, ilast); 00695 } 00696 if (ilz) { 00697 i__2 = *n; 00698 for (j = 1; j <= i__2; ++j) { 00699 z___ref(j, ilast) = -z___ref(j, ilast); 00700 /* L100: */ 00701 } 00702 } 00703 } 00704 alphar[ilast] = a_ref(ilast, ilast); 00705 alphai[ilast] = 0.; 00706 beta[ilast] = b_ref(ilast, ilast); 00707 00708 /* Go to next block -- exit if finished. */ 00709 00710 --ilast; 00711 if (ilast < *ilo) { 00712 goto L380; 00713 } 00714 00715 /* Reset counters */ 00716 00717 iiter = 0; 00718 eshift = 0.; 00719 if (! ilschr) { 00720 ilastm = ilast; 00721 if (ifrstm > ilast) { 00722 ifrstm = *ilo; 00723 } 00724 } 00725 goto L350; 00726 00727 /* QZ step 00728 00729 This iteration only involves rows/columns IFIRST:ILAST. We 00730 assume IFIRST < ILAST, and that the diagonal of B is non-zero. */ 00731 00732 L110: 00733 ++iiter; 00734 if (! ilschr) { 00735 ifrstm = ifirst; 00736 } 00737 00738 /* Compute single shifts. 00739 00740 At this point, IFIRST < ILAST, and the diagonal elements of 00741 B(IFIRST:ILAST,IFIRST,ILAST) are larger than BTOL (in 00742 magnitude) */ 00743 00744 if (iiter / 10 * 10 == iiter) { 00745 00746 /* Exceptional shift. Chosen for no particularly good reason. 00747 (Single shift only.) */ 00748 00749 if ((Treal) maxit * safmin * (d__1 = a_ref(ilast - 1, ilast), 00750 absMACRO(d__1)) < (d__2 = b_ref(ilast - 1, ilast - 1), absMACRO( 00751 d__2))) { 00752 eshift += a_ref(ilast - 1, ilast) / b_ref(ilast - 1, ilast - 00753 1); 00754 } else { 00755 eshift += 1. / (safmin * (Treal) maxit); 00756 } 00757 s1 = 1.; 00758 wr = eshift; 00759 00760 } else { 00761 00762 /* Shifts based on the generalized eigenvalues of the 00763 bottom-right 2x2 block of A and B. The first eigenvalue 00764 returned by DLAG2 is the Wilkinson shift (AEP p.512), */ 00765 00766 d__1 = safmin * 100.; 00767 template_lapack_lag2(&a_ref(ilast - 1, ilast - 1), lda, &b_ref(ilast - 1, ilast 00768 - 1), ldb, &d__1, &s1, &s2, &wr, &wr2, &wi); 00769 00770 /* Computing MAX 00771 Computing MAX */ 00772 d__3 = 1., d__4 = absMACRO(wr), d__3 = maxMACRO(d__3,d__4), d__4 = absMACRO(wi); 00773 d__1 = s1, d__2 = safmin * maxMACRO(d__3,d__4); 00774 temp = maxMACRO(d__1,d__2); 00775 if (wi != 0.) { 00776 goto L200; 00777 } 00778 } 00779 00780 /* Fiddle with shift to avoid overflow */ 00781 00782 temp = minMACRO(ascale,1.) * (safmax * .5); 00783 if (s1 > temp) { 00784 scale = temp / s1; 00785 } else { 00786 scale = 1.; 00787 } 00788 00789 temp = minMACRO(bscale,1.) * (safmax * .5); 00790 if (absMACRO(wr) > temp) { 00791 /* Computing MIN */ 00792 d__1 = scale, d__2 = temp / absMACRO(wr); 00793 scale = minMACRO(d__1,d__2); 00794 } 00795 s1 = scale * s1; 00796 wr = scale * wr; 00797 00798 /* Now check for two consecutive small subdiagonals. */ 00799 00800 i__2 = ifirst + 1; 00801 for (j = ilast - 1; j >= i__2; --j) { 00802 istart = j; 00803 temp = (d__1 = s1 * a_ref(j, j - 1), absMACRO(d__1)); 00804 temp2 = (d__1 = s1 * a_ref(j, j) - wr * b_ref(j, j), absMACRO(d__1)); 00805 tempr = maxMACRO(temp,temp2); 00806 if (tempr < 1. && tempr != 0.) { 00807 temp /= tempr; 00808 temp2 /= tempr; 00809 } 00810 if ((d__1 = ascale * a_ref(j + 1, j) * temp, absMACRO(d__1)) <= ascale 00811 * atol * temp2) { 00812 goto L130; 00813 } 00814 /* L120: */ 00815 } 00816 00817 istart = ifirst; 00818 L130: 00819 00820 /* Do an implicit single-shift QZ sweep. 00821 00822 Initial Q */ 00823 00824 temp = s1 * a_ref(istart, istart) - wr * b_ref(istart, istart); 00825 temp2 = s1 * a_ref(istart + 1, istart); 00826 template_lapack_lartg(&temp, &temp2, &c__, &s, &tempr); 00827 00828 /* Sweep */ 00829 00830 i__2 = ilast - 1; 00831 for (j = istart; j <= i__2; ++j) { 00832 if (j > istart) { 00833 temp = a_ref(j, j - 1); 00834 template_lapack_lartg(&temp, &a_ref(j + 1, j - 1), &c__, &s, &a_ref(j, j - 00835 1)); 00836 a_ref(j + 1, j - 1) = 0.; 00837 } 00838 00839 i__3 = ilastm; 00840 for (jc = j; jc <= i__3; ++jc) { 00841 temp = c__ * a_ref(j, jc) + s * a_ref(j + 1, jc); 00842 a_ref(j + 1, jc) = -s * a_ref(j, jc) + c__ * a_ref(j + 1, jc); 00843 a_ref(j, jc) = temp; 00844 temp2 = c__ * b_ref(j, jc) + s * b_ref(j + 1, jc); 00845 b_ref(j + 1, jc) = -s * b_ref(j, jc) + c__ * b_ref(j + 1, jc); 00846 b_ref(j, jc) = temp2; 00847 /* L140: */ 00848 } 00849 if (ilq) { 00850 i__3 = *n; 00851 for (jr = 1; jr <= i__3; ++jr) { 00852 temp = c__ * q_ref(jr, j) + s * q_ref(jr, j + 1); 00853 q_ref(jr, j + 1) = -s * q_ref(jr, j) + c__ * q_ref(jr, j 00854 + 1); 00855 q_ref(jr, j) = temp; 00856 /* L150: */ 00857 } 00858 } 00859 00860 temp = b_ref(j + 1, j + 1); 00861 template_lapack_lartg(&temp, &b_ref(j + 1, j), &c__, &s, &b_ref(j + 1, j + 1)); 00862 b_ref(j + 1, j) = 0.; 00863 00864 /* Computing MIN */ 00865 i__4 = j + 2; 00866 i__3 = minMACRO(i__4,ilast); 00867 for (jr = ifrstm; jr <= i__3; ++jr) { 00868 temp = c__ * a_ref(jr, j + 1) + s * a_ref(jr, j); 00869 a_ref(jr, j) = -s * a_ref(jr, j + 1) + c__ * a_ref(jr, j); 00870 a_ref(jr, j + 1) = temp; 00871 /* L160: */ 00872 } 00873 i__3 = j; 00874 for (jr = ifrstm; jr <= i__3; ++jr) { 00875 temp = c__ * b_ref(jr, j + 1) + s * b_ref(jr, j); 00876 b_ref(jr, j) = -s * b_ref(jr, j + 1) + c__ * b_ref(jr, j); 00877 b_ref(jr, j + 1) = temp; 00878 /* L170: */ 00879 } 00880 if (ilz) { 00881 i__3 = *n; 00882 for (jr = 1; jr <= i__3; ++jr) { 00883 temp = c__ * z___ref(jr, j + 1) + s * z___ref(jr, j); 00884 z___ref(jr, j) = -s * z___ref(jr, j + 1) + c__ * z___ref( 00885 jr, j); 00886 z___ref(jr, j + 1) = temp; 00887 /* L180: */ 00888 } 00889 } 00890 /* L190: */ 00891 } 00892 00893 goto L350; 00894 00895 /* Use Francis double-shift 00896 00897 Note: the Francis double-shift should work with real shifts, 00898 but only if the block is at least 3x3. 00899 This code may break if this point is reached with 00900 a 2x2 block with real eigenvalues. */ 00901 00902 L200: 00903 if (ifirst + 1 == ilast) { 00904 00905 /* Special case -- 2x2 block with complex eigenvectors 00906 00907 Step 1: Standardize, that is, rotate so that 00908 00909 ( B11 0 ) 00910 B = ( ) with B11 non-negative. 00911 ( 0 B22 ) */ 00912 00913 template_lapack_lasv2(&b_ref(ilast - 1, ilast - 1), &b_ref(ilast - 1, ilast), & 00914 b_ref(ilast, ilast), &b22, &b11, &sr, &cr, &sl, &cl); 00915 00916 if (b11 < 0.) { 00917 cr = -cr; 00918 sr = -sr; 00919 b11 = -b11; 00920 b22 = -b22; 00921 } 00922 00923 i__2 = ilastm + 1 - ifirst; 00924 template_blas_rot(&i__2, &a_ref(ilast - 1, ilast - 1), lda, &a_ref(ilast, 00925 ilast - 1), lda, &cl, &sl); 00926 i__2 = ilast + 1 - ifrstm; 00927 template_blas_rot(&i__2, &a_ref(ifrstm, ilast - 1), &c__1, &a_ref(ifrstm, 00928 ilast), &c__1, &cr, &sr); 00929 00930 if (ilast < ilastm) { 00931 i__2 = ilastm - ilast; 00932 template_blas_rot(&i__2, &b_ref(ilast - 1, ilast + 1), ldb, &b_ref(ilast, 00933 ilast + 1), lda, &cl, &sl); 00934 } 00935 if (ifrstm < ilast - 1) { 00936 i__2 = ifirst - ifrstm; 00937 template_blas_rot(&i__2, &b_ref(ifrstm, ilast - 1), &c__1, &b_ref(ifrstm, 00938 ilast), &c__1, &cr, &sr); 00939 } 00940 00941 if (ilq) { 00942 template_blas_rot(n, &q_ref(1, ilast - 1), &c__1, &q_ref(1, ilast), &c__1, 00943 &cl, &sl); 00944 } 00945 if (ilz) { 00946 template_blas_rot(n, &z___ref(1, ilast - 1), &c__1, &z___ref(1, ilast), & 00947 c__1, &cr, &sr); 00948 } 00949 00950 b_ref(ilast - 1, ilast - 1) = b11; 00951 b_ref(ilast - 1, ilast) = 0.; 00952 b_ref(ilast, ilast - 1) = 0.; 00953 b_ref(ilast, ilast) = b22; 00954 00955 /* If B22 is negative, negate column ILAST */ 00956 00957 if (b22 < 0.) { 00958 i__2 = ilast; 00959 for (j = ifrstm; j <= i__2; ++j) { 00960 a_ref(j, ilast) = -a_ref(j, ilast); 00961 b_ref(j, ilast) = -b_ref(j, ilast); 00962 /* L210: */ 00963 } 00964 00965 if (ilz) { 00966 i__2 = *n; 00967 for (j = 1; j <= i__2; ++j) { 00968 z___ref(j, ilast) = -z___ref(j, ilast); 00969 /* L220: */ 00970 } 00971 } 00972 } 00973 00974 /* Step 2: Compute ALPHAR, ALPHAI, and BETA (see refs.) 00975 00976 Recompute shift */ 00977 00978 d__1 = safmin * 100.; 00979 template_lapack_lag2(&a_ref(ilast - 1, ilast - 1), lda, &b_ref(ilast - 1, ilast 00980 - 1), ldb, &d__1, &s1, &temp, &wr, &temp2, &wi); 00981 00982 /* If standardization has perturbed the shift onto real line, 00983 do another (real single-shift) QR step. */ 00984 00985 if (wi == 0.) { 00986 goto L350; 00987 } 00988 s1inv = 1. / s1; 00989 00990 /* Do EISPACK (QZVAL) computation of alpha and beta */ 00991 00992 a11 = a_ref(ilast - 1, ilast - 1); 00993 a21 = a_ref(ilast, ilast - 1); 00994 a12 = a_ref(ilast - 1, ilast); 00995 a22 = a_ref(ilast, ilast); 00996 00997 /* Compute complex Givens rotation on right 00998 (Assume some element of C = (sA - wB) > unfl ) 00999 __ 01000 (sA - wB) ( CZ -SZ ) 01001 ( SZ CZ ) */ 01002 01003 c11r = s1 * a11 - wr * b11; 01004 c11i = -wi * b11; 01005 c12 = s1 * a12; 01006 c21 = s1 * a21; 01007 c22r = s1 * a22 - wr * b22; 01008 c22i = -wi * b22; 01009 01010 if (absMACRO(c11r) + absMACRO(c11i) + absMACRO(c12) > absMACRO(c21) + absMACRO(c22r) + absMACRO( 01011 c22i)) { 01012 t = template_lapack_lapy3(&c12, &c11r, &c11i); 01013 cz = c12 / t; 01014 szr = -c11r / t; 01015 szi = -c11i / t; 01016 } else { 01017 cz = template_lapack_lapy2(&c22r, &c22i); 01018 if (cz <= safmin) { 01019 cz = 0.; 01020 szr = 1.; 01021 szi = 0.; 01022 } else { 01023 tempr = c22r / cz; 01024 tempi = c22i / cz; 01025 t = template_lapack_lapy2(&cz, &c21); 01026 cz /= t; 01027 szr = -c21 * tempr / t; 01028 szi = c21 * tempi / t; 01029 } 01030 } 01031 01032 /* Compute Givens rotation on left 01033 01034 ( CQ SQ ) 01035 ( __ ) A or B 01036 ( -SQ CQ ) */ 01037 01038 an = absMACRO(a11) + absMACRO(a12) + absMACRO(a21) + absMACRO(a22); 01039 bn = absMACRO(b11) + absMACRO(b22); 01040 wabs = absMACRO(wr) + absMACRO(wi); 01041 if (s1 * an > wabs * bn) { 01042 cq = cz * b11; 01043 sqr = szr * b22; 01044 sqi = -szi * b22; 01045 } else { 01046 a1r = cz * a11 + szr * a12; 01047 a1i = szi * a12; 01048 a2r = cz * a21 + szr * a22; 01049 a2i = szi * a22; 01050 cq = template_lapack_lapy2(&a1r, &a1i); 01051 if (cq <= safmin) { 01052 cq = 0.; 01053 sqr = 1.; 01054 sqi = 0.; 01055 } else { 01056 tempr = a1r / cq; 01057 tempi = a1i / cq; 01058 sqr = tempr * a2r + tempi * a2i; 01059 sqi = tempi * a2r - tempr * a2i; 01060 } 01061 } 01062 t = template_lapack_lapy3(&cq, &sqr, &sqi); 01063 cq /= t; 01064 sqr /= t; 01065 sqi /= t; 01066 01067 /* Compute diagonal elements of QBZ */ 01068 01069 tempr = sqr * szr - sqi * szi; 01070 tempi = sqr * szi + sqi * szr; 01071 b1r = cq * cz * b11 + tempr * b22; 01072 b1i = tempi * b22; 01073 b1a = template_lapack_lapy2(&b1r, &b1i); 01074 b2r = cq * cz * b22 + tempr * b11; 01075 b2i = -tempi * b11; 01076 b2a = template_lapack_lapy2(&b2r, &b2i); 01077 01078 /* Normalize so beta > 0, and Im( alpha1 ) > 0 */ 01079 01080 beta[ilast - 1] = b1a; 01081 beta[ilast] = b2a; 01082 alphar[ilast - 1] = wr * b1a * s1inv; 01083 alphai[ilast - 1] = wi * b1a * s1inv; 01084 alphar[ilast] = wr * b2a * s1inv; 01085 alphai[ilast] = -(wi * b2a) * s1inv; 01086 01087 /* Step 3: Go to next block -- exit if finished. */ 01088 01089 ilast = ifirst - 1; 01090 if (ilast < *ilo) { 01091 goto L380; 01092 } 01093 01094 /* Reset counters */ 01095 01096 iiter = 0; 01097 eshift = 0.; 01098 if (! ilschr) { 01099 ilastm = ilast; 01100 if (ifrstm > ilast) { 01101 ifrstm = *ilo; 01102 } 01103 } 01104 goto L350; 01105 } else { 01106 01107 /* Usual case: 3x3 or larger block, using Francis implicit 01108 double-shift 01109 01110 2 01111 Eigenvalue equation is w - c w + d = 0, 01112 01113 -1 2 -1 01114 so compute 1st column of (A B ) - c A B + d 01115 using the formula in QZIT (from EISPACK) 01116 01117 We assume that the block is at least 3x3 */ 01118 01119 ad11 = ascale * a_ref(ilast - 1, ilast - 1) / (bscale * b_ref( 01120 ilast - 1, ilast - 1)); 01121 ad21 = ascale * a_ref(ilast, ilast - 1) / (bscale * b_ref(ilast - 01122 1, ilast - 1)); 01123 ad12 = ascale * a_ref(ilast - 1, ilast) / (bscale * b_ref(ilast, 01124 ilast)); 01125 ad22 = ascale * a_ref(ilast, ilast) / (bscale * b_ref(ilast, 01126 ilast)); 01127 u12 = b_ref(ilast - 1, ilast) / b_ref(ilast, ilast); 01128 ad11l = ascale * a_ref(ifirst, ifirst) / (bscale * b_ref(ifirst, 01129 ifirst)); 01130 ad21l = ascale * a_ref(ifirst + 1, ifirst) / (bscale * b_ref( 01131 ifirst, ifirst)); 01132 ad12l = ascale * a_ref(ifirst, ifirst + 1) / (bscale * b_ref( 01133 ifirst + 1, ifirst + 1)); 01134 ad22l = ascale * a_ref(ifirst + 1, ifirst + 1) / (bscale * b_ref( 01135 ifirst + 1, ifirst + 1)); 01136 ad32l = ascale * a_ref(ifirst + 2, ifirst + 1) / (bscale * b_ref( 01137 ifirst + 1, ifirst + 1)); 01138 u12l = b_ref(ifirst, ifirst + 1) / b_ref(ifirst + 1, ifirst + 1); 01139 01140 v[0] = (ad11 - ad11l) * (ad22 - ad11l) - ad12 * ad21 + ad21 * u12 01141 * ad11l + (ad12l - ad11l * u12l) * ad21l; 01142 v[1] = (ad22l - ad11l - ad21l * u12l - (ad11 - ad11l) - (ad22 - 01143 ad11l) + ad21 * u12) * ad21l; 01144 v[2] = ad32l * ad21l; 01145 01146 istart = ifirst; 01147 01148 template_lapack_larfg(&c__3, v, &v[1], &c__1, &tau); 01149 v[0] = 1.; 01150 01151 /* Sweep */ 01152 01153 i__2 = ilast - 2; 01154 for (j = istart; j <= i__2; ++j) { 01155 01156 /* All but last elements: use 3x3 Householder transforms. 01157 01158 Zero (j-1)st column of A */ 01159 01160 if (j > istart) { 01161 v[0] = a_ref(j, j - 1); 01162 v[1] = a_ref(j + 1, j - 1); 01163 v[2] = a_ref(j + 2, j - 1); 01164 01165 template_lapack_larfg(&c__3, &a_ref(j, j - 1), &v[1], &c__1, &tau); 01166 v[0] = 1.; 01167 a_ref(j + 1, j - 1) = 0.; 01168 a_ref(j + 2, j - 1) = 0.; 01169 } 01170 01171 i__3 = ilastm; 01172 for (jc = j; jc <= i__3; ++jc) { 01173 temp = tau * (a_ref(j, jc) + v[1] * a_ref(j + 1, jc) + v[ 01174 2] * a_ref(j + 2, jc)); 01175 a_ref(j, jc) = a_ref(j, jc) - temp; 01176 a_ref(j + 1, jc) = a_ref(j + 1, jc) - temp * v[1]; 01177 a_ref(j + 2, jc) = a_ref(j + 2, jc) - temp * v[2]; 01178 temp2 = tau * (b_ref(j, jc) + v[1] * b_ref(j + 1, jc) + v[ 01179 2] * b_ref(j + 2, jc)); 01180 b_ref(j, jc) = b_ref(j, jc) - temp2; 01181 b_ref(j + 1, jc) = b_ref(j + 1, jc) - temp2 * v[1]; 01182 b_ref(j + 2, jc) = b_ref(j + 2, jc) - temp2 * v[2]; 01183 /* L230: */ 01184 } 01185 if (ilq) { 01186 i__3 = *n; 01187 for (jr = 1; jr <= i__3; ++jr) { 01188 temp = tau * (q_ref(jr, j) + v[1] * q_ref(jr, j + 1) 01189 + v[2] * q_ref(jr, j + 2)); 01190 q_ref(jr, j) = q_ref(jr, j) - temp; 01191 q_ref(jr, j + 1) = q_ref(jr, j + 1) - temp * v[1]; 01192 q_ref(jr, j + 2) = q_ref(jr, j + 2) - temp * v[2]; 01193 /* L240: */ 01194 } 01195 } 01196 01197 /* Zero j-th column of B (see DLAGBC for details) 01198 01199 Swap rows to pivot */ 01200 01201 ilpivt = FALSE_; 01202 /* Computing MAX */ 01203 d__3 = (d__1 = b_ref(j + 1, j + 1), absMACRO(d__1)), d__4 = (d__2 = 01204 b_ref(j + 1, j + 2), absMACRO(d__2)); 01205 temp = maxMACRO(d__3,d__4); 01206 /* Computing MAX */ 01207 d__3 = (d__1 = b_ref(j + 2, j + 1), absMACRO(d__1)), d__4 = (d__2 = 01208 b_ref(j + 2, j + 2), absMACRO(d__2)); 01209 temp2 = maxMACRO(d__3,d__4); 01210 if (maxMACRO(temp,temp2) < safmin) { 01211 scale = 0.; 01212 u1 = 1.; 01213 u2 = 0.; 01214 goto L250; 01215 } else if (temp >= temp2) { 01216 w11 = b_ref(j + 1, j + 1); 01217 w21 = b_ref(j + 2, j + 1); 01218 w12 = b_ref(j + 1, j + 2); 01219 w22 = b_ref(j + 2, j + 2); 01220 u1 = b_ref(j + 1, j); 01221 u2 = b_ref(j + 2, j); 01222 } else { 01223 w21 = b_ref(j + 1, j + 1); 01224 w11 = b_ref(j + 2, j + 1); 01225 w22 = b_ref(j + 1, j + 2); 01226 w12 = b_ref(j + 2, j + 2); 01227 u2 = b_ref(j + 1, j); 01228 u1 = b_ref(j + 2, j); 01229 } 01230 01231 /* Swap columns if nec. */ 01232 01233 if (absMACRO(w12) > absMACRO(w11)) { 01234 ilpivt = TRUE_; 01235 temp = w12; 01236 temp2 = w22; 01237 w12 = w11; 01238 w22 = w21; 01239 w11 = temp; 01240 w21 = temp2; 01241 } 01242 01243 /* LU-factor */ 01244 01245 temp = w21 / w11; 01246 u2 -= temp * u1; 01247 w22 -= temp * w12; 01248 w21 = 0.; 01249 01250 /* Compute SCALE */ 01251 01252 scale = 1.; 01253 if (absMACRO(w22) < safmin) { 01254 scale = 0.; 01255 u2 = 1.; 01256 u1 = -w12 / w11; 01257 goto L250; 01258 } 01259 if (absMACRO(w22) < absMACRO(u2)) { 01260 scale = (d__1 = w22 / u2, absMACRO(d__1)); 01261 } 01262 if (absMACRO(w11) < absMACRO(u1)) { 01263 /* Computing MIN */ 01264 d__2 = scale, d__3 = (d__1 = w11 / u1, absMACRO(d__1)); 01265 scale = minMACRO(d__2,d__3); 01266 } 01267 01268 /* Solve */ 01269 01270 u2 = scale * u2 / w22; 01271 u1 = (scale * u1 - w12 * u2) / w11; 01272 01273 L250: 01274 if (ilpivt) { 01275 temp = u2; 01276 u2 = u1; 01277 u1 = temp; 01278 } 01279 01280 /* Compute Householder Vector 01281 01282 Computing 2nd power */ 01283 d__1 = scale; 01284 /* Computing 2nd power */ 01285 d__2 = u1; 01286 /* Computing 2nd power */ 01287 d__3 = u2; 01288 t = template_blas_sqrt(d__1 * d__1 + d__2 * d__2 + d__3 * d__3); 01289 tau = scale / t + 1.; 01290 vs = -1. / (scale + t); 01291 v[0] = 1.; 01292 v[1] = vs * u1; 01293 v[2] = vs * u2; 01294 01295 /* Apply transformations from the right. 01296 01297 Computing MIN */ 01298 i__4 = j + 3; 01299 i__3 = minMACRO(i__4,ilast); 01300 for (jr = ifrstm; jr <= i__3; ++jr) { 01301 temp = tau * (a_ref(jr, j) + v[1] * a_ref(jr, j + 1) + v[ 01302 2] * a_ref(jr, j + 2)); 01303 a_ref(jr, j) = a_ref(jr, j) - temp; 01304 a_ref(jr, j + 1) = a_ref(jr, j + 1) - temp * v[1]; 01305 a_ref(jr, j + 2) = a_ref(jr, j + 2) - temp * v[2]; 01306 /* L260: */ 01307 } 01308 i__3 = j + 2; 01309 for (jr = ifrstm; jr <= i__3; ++jr) { 01310 temp = tau * (b_ref(jr, j) + v[1] * b_ref(jr, j + 1) + v[ 01311 2] * b_ref(jr, j + 2)); 01312 b_ref(jr, j) = b_ref(jr, j) - temp; 01313 b_ref(jr, j + 1) = b_ref(jr, j + 1) - temp * v[1]; 01314 b_ref(jr, j + 2) = b_ref(jr, j + 2) - temp * v[2]; 01315 /* L270: */ 01316 } 01317 if (ilz) { 01318 i__3 = *n; 01319 for (jr = 1; jr <= i__3; ++jr) { 01320 temp = tau * (z___ref(jr, j) + v[1] * z___ref(jr, j + 01321 1) + v[2] * z___ref(jr, j + 2)); 01322 z___ref(jr, j) = z___ref(jr, j) - temp; 01323 z___ref(jr, j + 1) = z___ref(jr, j + 1) - temp * v[1]; 01324 z___ref(jr, j + 2) = z___ref(jr, j + 2) - temp * v[2]; 01325 /* L280: */ 01326 } 01327 } 01328 b_ref(j + 1, j) = 0.; 01329 b_ref(j + 2, j) = 0.; 01330 /* L290: */ 01331 } 01332 01333 /* Last elements: Use Givens rotations 01334 01335 Rotations from the left */ 01336 01337 j = ilast - 1; 01338 temp = a_ref(j, j - 1); 01339 template_lapack_lartg(&temp, &a_ref(j + 1, j - 1), &c__, &s, &a_ref(j, j - 1)); 01340 a_ref(j + 1, j - 1) = 0.; 01341 01342 i__2 = ilastm; 01343 for (jc = j; jc <= i__2; ++jc) { 01344 temp = c__ * a_ref(j, jc) + s * a_ref(j + 1, jc); 01345 a_ref(j + 1, jc) = -s * a_ref(j, jc) + c__ * a_ref(j + 1, jc); 01346 a_ref(j, jc) = temp; 01347 temp2 = c__ * b_ref(j, jc) + s * b_ref(j + 1, jc); 01348 b_ref(j + 1, jc) = -s * b_ref(j, jc) + c__ * b_ref(j + 1, jc); 01349 b_ref(j, jc) = temp2; 01350 /* L300: */ 01351 } 01352 if (ilq) { 01353 i__2 = *n; 01354 for (jr = 1; jr <= i__2; ++jr) { 01355 temp = c__ * q_ref(jr, j) + s * q_ref(jr, j + 1); 01356 q_ref(jr, j + 1) = -s * q_ref(jr, j) + c__ * q_ref(jr, j 01357 + 1); 01358 q_ref(jr, j) = temp; 01359 /* L310: */ 01360 } 01361 } 01362 01363 /* Rotations from the right. */ 01364 01365 temp = b_ref(j + 1, j + 1); 01366 template_lapack_lartg(&temp, &b_ref(j + 1, j), &c__, &s, &b_ref(j + 1, j + 1)); 01367 b_ref(j + 1, j) = 0.; 01368 01369 i__2 = ilast; 01370 for (jr = ifrstm; jr <= i__2; ++jr) { 01371 temp = c__ * a_ref(jr, j + 1) + s * a_ref(jr, j); 01372 a_ref(jr, j) = -s * a_ref(jr, j + 1) + c__ * a_ref(jr, j); 01373 a_ref(jr, j + 1) = temp; 01374 /* L320: */ 01375 } 01376 i__2 = ilast - 1; 01377 for (jr = ifrstm; jr <= i__2; ++jr) { 01378 temp = c__ * b_ref(jr, j + 1) + s * b_ref(jr, j); 01379 b_ref(jr, j) = -s * b_ref(jr, j + 1) + c__ * b_ref(jr, j); 01380 b_ref(jr, j + 1) = temp; 01381 /* L330: */ 01382 } 01383 if (ilz) { 01384 i__2 = *n; 01385 for (jr = 1; jr <= i__2; ++jr) { 01386 temp = c__ * z___ref(jr, j + 1) + s * z___ref(jr, j); 01387 z___ref(jr, j) = -s * z___ref(jr, j + 1) + c__ * z___ref( 01388 jr, j); 01389 z___ref(jr, j + 1) = temp; 01390 /* L340: */ 01391 } 01392 } 01393 01394 /* End of Double-Shift code */ 01395 01396 } 01397 01398 goto L350; 01399 01400 /* End of iteration loop */ 01401 01402 L350: 01403 /* L360: */ 01404 ; 01405 } 01406 01407 /* Drop-through = non-convergence 01408 01409 L370: */ 01410 *info = ilast; 01411 goto L420; 01412 01413 /* Successful completion of all QZ steps */ 01414 01415 L380: 01416 01417 /* Set Eigenvalues 1:ILO-1 */ 01418 01419 i__1 = *ilo - 1; 01420 for (j = 1; j <= i__1; ++j) { 01421 if (b_ref(j, j) < 0.) { 01422 if (ilschr) { 01423 i__2 = j; 01424 for (jr = 1; jr <= i__2; ++jr) { 01425 a_ref(jr, j) = -a_ref(jr, j); 01426 b_ref(jr, j) = -b_ref(jr, j); 01427 /* L390: */ 01428 } 01429 } else { 01430 a_ref(j, j) = -a_ref(j, j); 01431 b_ref(j, j) = -b_ref(j, j); 01432 } 01433 if (ilz) { 01434 i__2 = *n; 01435 for (jr = 1; jr <= i__2; ++jr) { 01436 z___ref(jr, j) = -z___ref(jr, j); 01437 /* L400: */ 01438 } 01439 } 01440 } 01441 alphar[j] = a_ref(j, j); 01442 alphai[j] = 0.; 01443 beta[j] = b_ref(j, j); 01444 /* L410: */ 01445 } 01446 01447 /* Normal Termination */ 01448 01449 *info = 0; 01450 01451 /* Exit (other than argument error) -- return optimal workspace size */ 01452 01453 L420: 01454 work[1] = (Treal) (*n); 01455 return 0; 01456 01457 /* End of DHGEQZ */ 01458 01459 } /* dhgeqz_ */ 01460 01461 #undef z___ref 01462 #undef q_ref 01463 #undef b_ref 01464 #undef a_ref 01465 01466 01467 #endif