ergo
template_lapack_hgeqz.h
Go to the documentation of this file.
1/* Ergo, version 3.8.2, a program for linear scaling electronic structure
2 * calculations.
3 * Copyright (C) 2023 Elias Rudberg, Emanuel H. Rubensson, Pawel Salek,
4 * and Anastasia Kruchinina.
5 *
6 * This program is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program. If not, see <http://www.gnu.org/licenses/>.
18 *
19 * Primary academic reference:
20 * Ergo: An open-source program for linear-scaling electronic structure
21 * calculations,
22 * Elias Rudberg, Emanuel H. Rubensson, Pawel Salek, and Anastasia
23 * Kruchinina,
24 * SoftwareX 7, 107 (2018),
25 * <http://dx.doi.org/10.1016/j.softx.2018.03.005>
26 *
27 * For further information about Ergo, see <http://www.ergoscf.org>.
28 */
29
30 /* This file belongs to the template_lapack part of the Ergo source
31 * code. The source files in the template_lapack directory are modified
32 * versions of files originally distributed as CLAPACK, see the
33 * Copyright/license notice in the file template_lapack/COPYING.
34 */
35
36
37#ifndef TEMPLATE_LAPACK_HGEQZ_HEADER
38#define TEMPLATE_LAPACK_HGEQZ_HEADER
39
40
41template<class Treal>
42int template_lapack_hgeqz(const char *job, const char *compq, const char *compz, const integer *n,
43 const integer *ilo, const integer *ihi, Treal *a, const integer *lda, Treal *
44 b, const integer *ldb, Treal *alphar, Treal *alphai, Treal *
45 beta, Treal *q, const integer *ldq, Treal *z__, const integer *ldz,
46 Treal *work, const integer *lwork, integer *info)
47{
48/* -- LAPACK routine (version 3.0) --
49 Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
50 Courant Institute, Argonne National Lab, and Rice University
51 June 30, 1999
52
53
54 Purpose
55 =======
56
57 DHGEQZ implements a single-/double-shift version of the QZ method for
58 finding the generalized eigenvalues
59
60 w(j)=(ALPHAR(j) + i*ALPHAI(j))/BETAR(j) of the equation
61
62 det( A - w(i) B ) = 0
63
64 In addition, the pair A,B may be reduced to generalized Schur form:
65 B is upper triangular, and A is block upper triangular, where the
66 diagonal blocks are either 1-by-1 or 2-by-2, the 2-by-2 blocks having
67 complex generalized eigenvalues (see the description of the argument
68 JOB.)
69
70 If JOB='S', then the pair (A,B) is simultaneously reduced to Schur
71 form by applying one orthogonal tranformation (usually called Q) on
72 the left and another (usually called Z) on the right. The 2-by-2
73 upper-triangular diagonal blocks of B corresponding to 2-by-2 blocks
74 of A will be reduced to positive diagonal matrices. (I.e.,
75 if A(j+1,j) is non-zero, then B(j+1,j)=B(j,j+1)=0 and B(j,j) and
76 B(j+1,j+1) will be positive.)
77
78 If JOB='E', then at each iteration, the same transformations
79 are computed, but they are only applied to those parts of A and B
80 which are needed to compute ALPHAR, ALPHAI, and BETAR.
81
82 If JOB='S' and COMPQ and COMPZ are 'V' or 'I', then the orthogonal
83 transformations used to reduce (A,B) are accumulated into the arrays
84 Q and Z s.t.:
85
86 Q(in) A(in) Z(in)* = Q(out) A(out) Z(out)*
87 Q(in) B(in) Z(in)* = Q(out) B(out) Z(out)*
88
89 Ref: C.B. Moler & G.W. Stewart, "An Algorithm for Generalized Matrix
90 Eigenvalue Problems", SIAM J. Numer. Anal., 10(1973),
91 pp. 241--256.
92
93 Arguments
94 =========
95
96 JOB (input) CHARACTER*1
97 = 'E': compute only ALPHAR, ALPHAI, and BETA. A and B will
98 not necessarily be put into generalized Schur form.
99 = 'S': put A and B into generalized Schur form, as well
100 as computing ALPHAR, ALPHAI, and BETA.
101
102 COMPQ (input) CHARACTER*1
103 = 'N': do not modify Q.
104 = 'V': multiply the array Q on the right by the transpose of
105 the orthogonal tranformation that is applied to the
106 left side of A and B to reduce them to Schur form.
107 = 'I': like COMPQ='V', except that Q will be initialized to
108 the identity first.
109
110 COMPZ (input) CHARACTER*1
111 = 'N': do not modify Z.
112 = 'V': multiply the array Z on the right by the orthogonal
113 tranformation that is applied to the right side of
114 A and B to reduce them to Schur form.
115 = 'I': like COMPZ='V', except that Z will be initialized to
116 the identity first.
117
118 N (input) INTEGER
119 The order of the matrices A, B, Q, and Z. N >= 0.
120
121 ILO (input) INTEGER
122 IHI (input) INTEGER
123 It is assumed that A is already upper triangular in rows and
124 columns 1:ILO-1 and IHI+1:N.
125 1 <= ILO <= IHI <= N, if N > 0; ILO=1 and IHI=0, if N=0.
126
127 A (input/output) DOUBLE PRECISION array, dimension (LDA, N)
128 On entry, the N-by-N upper Hessenberg matrix A. Elements
129 below the subdiagonal must be zero.
130 If JOB='S', then on exit A and B will have been
131 simultaneously reduced to generalized Schur form.
132 If JOB='E', then on exit A will have been destroyed.
133 The diagonal blocks will be correct, but the off-diagonal
134 portion will be meaningless.
135
136 LDA (input) INTEGER
137 The leading dimension of the array A. LDA >= max( 1, N ).
138
139 B (input/output) DOUBLE PRECISION array, dimension (LDB, N)
140 On entry, the N-by-N upper triangular matrix B. Elements
141 below the diagonal must be zero. 2-by-2 blocks in B
142 corresponding to 2-by-2 blocks in A will be reduced to
143 positive diagonal form. (I.e., if A(j+1,j) is non-zero,
144 then B(j+1,j)=B(j,j+1)=0 and B(j,j) and B(j+1,j+1) will be
145 positive.)
146 If JOB='S', then on exit A and B will have been
147 simultaneously reduced to Schur form.
148 If JOB='E', then on exit B will have been destroyed.
149 Elements corresponding to diagonal blocks of A will be
150 correct, but the off-diagonal portion will be meaningless.
151
152 LDB (input) INTEGER
153 The leading dimension of the array B. LDB >= max( 1, N ).
154
155 ALPHAR (output) DOUBLE PRECISION array, dimension (N)
156 ALPHAR(1:N) will be set to real parts of the diagonal
157 elements of A that would result from reducing A and B to
158 Schur form and then further reducing them both to triangular
159 form using unitary transformations s.t. the diagonal of B
160 was non-negative real. Thus, if A(j,j) is in a 1-by-1 block
161 (i.e., A(j+1,j)=A(j,j+1)=0), then ALPHAR(j)=A(j,j).
162 Note that the (real or complex) values
163 (ALPHAR(j) + i*ALPHAI(j))/BETA(j), j=1,...,N, are the
164 generalized eigenvalues of the matrix pencil A - wB.
165
166 ALPHAI (output) DOUBLE PRECISION array, dimension (N)
167 ALPHAI(1:N) will be set to imaginary parts of the diagonal
168 elements of A that would result from reducing A and B to
169 Schur form and then further reducing them both to triangular
170 form using unitary transformations s.t. the diagonal of B
171 was non-negative real. Thus, if A(j,j) is in a 1-by-1 block
172 (i.e., A(j+1,j)=A(j,j+1)=0), then ALPHAR(j)=0.
173 Note that the (real or complex) values
174 (ALPHAR(j) + i*ALPHAI(j))/BETA(j), j=1,...,N, are the
175 generalized eigenvalues of the matrix pencil A - wB.
176
177 BETA (output) DOUBLE PRECISION array, dimension (N)
178 BETA(1:N) will be set to the (real) diagonal elements of B
179 that would result from reducing A and B to Schur form and
180 then further reducing them both to triangular form using
181 unitary transformations s.t. the diagonal of B was
182 non-negative real. Thus, if A(j,j) is in a 1-by-1 block
183 (i.e., A(j+1,j)=A(j,j+1)=0), then BETA(j)=B(j,j).
184 Note that the (real or complex) values
185 (ALPHAR(j) + i*ALPHAI(j))/BETA(j), j=1,...,N, are the
186 generalized eigenvalues of the matrix pencil A - wB.
187 (Note that BETA(1:N) will always be non-negative, and no
188 BETAI is necessary.)
189
190 Q (input/output) DOUBLE PRECISION array, dimension (LDQ, N)
191 If COMPQ='N', then Q will not be referenced.
192 If COMPQ='V' or 'I', then the transpose of the orthogonal
193 transformations which are applied to A and B on the left
194 will be applied to the array Q on the right.
195
196 LDQ (input) INTEGER
197 The leading dimension of the array Q. LDQ >= 1.
198 If COMPQ='V' or 'I', then LDQ >= N.
199
200 Z (input/output) DOUBLE PRECISION array, dimension (LDZ, N)
201 If COMPZ='N', then Z will not be referenced.
202 If COMPZ='V' or 'I', then the orthogonal transformations
203 which are applied to A and B on the right will be applied
204 to the array Z on the right.
205
206 LDZ (input) INTEGER
207 The leading dimension of the array Z. LDZ >= 1.
208 If COMPZ='V' or 'I', then LDZ >= N.
209
210 WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
211 On exit, if INFO >= 0, WORK(1) returns the optimal LWORK.
212
213 LWORK (input) INTEGER
214 The dimension of the array WORK. LWORK >= max(1,N).
215
216 If LWORK = -1, then a workspace query is assumed; the routine
217 only calculates the optimal size of the WORK array, returns
218 this value as the first entry of the WORK array, and no error
219 message related to LWORK is issued by XERBLA.
220
221 INFO (output) INTEGER
222 = 0: successful exit
223 < 0: if INFO = -i, the i-th argument had an illegal value
224 = 1,...,N: the QZ iteration did not converge. (A,B) is not
225 in Schur form, but ALPHAR(i), ALPHAI(i), and
226 BETA(i), i=INFO+1,...,N should be correct.
227 = N+1,...,2*N: the shift calculation failed. (A,B) is not
228 in Schur form, but ALPHAR(i), ALPHAI(i), and
229 BETA(i), i=INFO-N+1,...,N should be correct.
230 > 2*N: various "impossible" errors.
231
232 Further Details
233 ===============
234
235 Iteration counters:
236
237 JITER -- counts iterations.
238 IITER -- counts iterations run since ILAST was last
239 changed. This is therefore reset only when a 1-by-1 or
240 2-by-2 block deflates off the bottom.
241
242 =====================================================================
243
244 $ SAFETY = 1.0E+0 )
245
246 Decode JOB, COMPQ, COMPZ
247
248 Parameter adjustments */
249 /* Table of constant values */
250 Treal c_b12 = 0.;
251 Treal c_b13 = 1.;
252 integer c__1 = 1;
253 integer c__3 = 3;
254
255 /* System generated locals */
256 integer a_dim1, a_offset, b_dim1, b_offset, q_dim1, q_offset, z_dim1,
257 z_offset, i__1, i__2, i__3, i__4;
258 Treal d__1, d__2, d__3, d__4;
259 /* Local variables */
260 Treal ad11l, ad12l, ad21l, ad22l, ad32l, wabs, atol, btol,
261 temp;
262 Treal temp2, s1inv, c__;
263 integer j;
264 Treal s, t, v[3], scale;
265 integer iiter, ilast, jiter;
266 Treal anorm, bnorm;
267 integer maxit;
268 Treal tempi, tempr, s1, s2, u1, u2;
269 logical ilazr2;
270 Treal a11, a12, a21, a22, b11, b22, c12, c21;
271 integer jc;
272 Treal an, bn, cl, cq, cr;
273 integer in;
274 Treal ascale, bscale, u12, w11;
275 integer jr;
276 Treal cz, sl, w12, w21, w22, wi;
277 Treal sr;
278 Treal vs, wr;
279 Treal safmin;
280 Treal safmax;
281 Treal eshift;
282 logical ilschr;
283 Treal b1a, b2a;
284 integer icompq, ilastm;
285 Treal a1i;
286 integer ischur;
287 Treal a2i, b1i;
288 logical ilazro;
289 integer icompz, ifirst;
290 Treal b2i;
291 integer ifrstm;
292 Treal a1r;
293 integer istart;
294 logical ilpivt;
295 Treal a2r, b1r, b2r;
296 logical lquery;
297 Treal wr2, ad11, ad12, ad21, ad22, c11i, c22i;
298 integer jch;
299 Treal c11r, c22r, u12l;
300 logical ilq;
301 Treal tau, sqi;
302 logical ilz;
303 Treal ulp, sqr, szi, szr;
304#define a_ref(a_1,a_2) a[(a_2)*a_dim1 + a_1]
305#define b_ref(a_1,a_2) b[(a_2)*b_dim1 + a_1]
306#define q_ref(a_1,a_2) q[(a_2)*q_dim1 + a_1]
307#define z___ref(a_1,a_2) z__[(a_2)*z_dim1 + a_1]
308
309
310 a_dim1 = *lda;
311 a_offset = 1 + a_dim1 * 1;
312 a -= a_offset;
313 b_dim1 = *ldb;
314 b_offset = 1 + b_dim1 * 1;
315 b -= b_offset;
316 --alphar;
317 --alphai;
318 --beta;
319 q_dim1 = *ldq;
320 q_offset = 1 + q_dim1 * 1;
321 q -= q_offset;
322 z_dim1 = *ldz;
323 z_offset = 1 + z_dim1 * 1;
324 z__ -= z_offset;
325 --work;
326
327 /* Initialization added by Elias to get rid of compiler warnings. */
328 ilschr = ilq = ilz = 0;
329 /* Function Body */
330 if (template_blas_lsame(job, "E")) {
331 ilschr = FALSE_;
332 ischur = 1;
333 } else if (template_blas_lsame(job, "S")) {
334 ilschr = TRUE_;
335 ischur = 2;
336 } else {
337 ischur = 0;
338 }
339
340 if (template_blas_lsame(compq, "N")) {
341 ilq = FALSE_;
342 icompq = 1;
343 } else if (template_blas_lsame(compq, "V")) {
344 ilq = TRUE_;
345 icompq = 2;
346 } else if (template_blas_lsame(compq, "I")) {
347 ilq = TRUE_;
348 icompq = 3;
349 } else {
350 icompq = 0;
351 }
352
353 if (template_blas_lsame(compz, "N")) {
354 ilz = FALSE_;
355 icompz = 1;
356 } else if (template_blas_lsame(compz, "V")) {
357 ilz = TRUE_;
358 icompz = 2;
359 } else if (template_blas_lsame(compz, "I")) {
360 ilz = TRUE_;
361 icompz = 3;
362 } else {
363 icompz = 0;
364 }
365
366/* Check Argument Values */
367
368 *info = 0;
369 work[1] = (Treal) maxMACRO(1,*n);
370 lquery = *lwork == -1;
371 if (ischur == 0) {
372 *info = -1;
373 } else if (icompq == 0) {
374 *info = -2;
375 } else if (icompz == 0) {
376 *info = -3;
377 } else if (*n < 0) {
378 *info = -4;
379 } else if (*ilo < 1) {
380 *info = -5;
381 } else if (*ihi > *n || *ihi < *ilo - 1) {
382 *info = -6;
383 } else if (*lda < *n) {
384 *info = -8;
385 } else if (*ldb < *n) {
386 *info = -10;
387 } else if (*ldq < 1 || ( ilq && *ldq < *n ) ) {
388 *info = -15;
389 } else if (*ldz < 1 || ( ilz && *ldz < *n ) ) {
390 *info = -17;
391 } else if (*lwork < maxMACRO(1,*n) && ! lquery) {
392 *info = -19;
393 }
394 if (*info != 0) {
395 i__1 = -(*info);
396 template_blas_erbla("HGEQZ ", &i__1);
397 return 0;
398 } else if (lquery) {
399 return 0;
400 }
401
402/* Quick return if possible */
403
404 if (*n <= 0) {
405 work[1] = 1.;
406 return 0;
407 }
408
409/* Initialize Q and Z */
410
411 if (icompq == 3) {
412 template_lapack_laset("Full", n, n, &c_b12, &c_b13, &q[q_offset], ldq);
413 }
414 if (icompz == 3) {
415 template_lapack_laset("Full", n, n, &c_b12, &c_b13, &z__[z_offset], ldz);
416 }
417
418/* Machine Constants */
419
420 in = *ihi + 1 - *ilo;
421 safmin = template_lapack_lamch("S", (Treal)0);
422 safmax = 1. / safmin;
423 ulp = template_lapack_lamch("E", (Treal)0) * template_lapack_lamch("B", (Treal)0);
424 anorm = dlanhs_("F", &in, &a_ref(*ilo, *ilo), lda, &work[1]);
425 bnorm = dlanhs_("F", &in, &b_ref(*ilo, *ilo), ldb, &work[1]);
426/* Computing MAX */
427 d__1 = safmin, d__2 = ulp * anorm;
428 atol = maxMACRO(d__1,d__2);
429/* Computing MAX */
430 d__1 = safmin, d__2 = ulp * bnorm;
431 btol = maxMACRO(d__1,d__2);
432 ascale = 1. / maxMACRO(safmin,anorm);
433 bscale = 1. / maxMACRO(safmin,bnorm);
434
435/* Set Eigenvalues IHI+1:N */
436
437 i__1 = *n;
438 for (j = *ihi + 1; j <= i__1; ++j) {
439 if (b_ref(j, j) < 0.) {
440 if (ilschr) {
441 i__2 = j;
442 for (jr = 1; jr <= i__2; ++jr) {
443 a_ref(jr, j) = -a_ref(jr, j);
444 b_ref(jr, j) = -b_ref(jr, j);
445/* L10: */
446 }
447 } else {
448 a_ref(j, j) = -a_ref(j, j);
449 b_ref(j, j) = -b_ref(j, j);
450 }
451 if (ilz) {
452 i__2 = *n;
453 for (jr = 1; jr <= i__2; ++jr) {
454 z___ref(jr, j) = -z___ref(jr, j);
455/* L20: */
456 }
457 }
458 }
459 alphar[j] = a_ref(j, j);
460 alphai[j] = 0.;
461 beta[j] = b_ref(j, j);
462/* L30: */
463 }
464
465/* If IHI < ILO, skip QZ steps */
466
467 if (*ihi < *ilo) {
468 goto L380;
469 }
470
471/* MAIN QZ ITERATION LOOP
472
473 Initialize dynamic indices
474
475 Eigenvalues ILAST+1:N have been found.
476 Column operations modify rows IFRSTM:whatever.
477 Row operations modify columns whatever:ILASTM.
478
479 If only eigenvalues are being computed, then
480 IFRSTM is the row of the last splitting row above row ILAST;
481 this is always at least ILO.
482 IITER counts iterations since the last eigenvalue was found,
483 to tell when to use an extraordinary shift.
484 MAXIT is the maximum number of QZ sweeps allowed. */
485
486 ilast = *ihi;
487 if (ilschr) {
488 ifrstm = 1;
489 ilastm = *n;
490 } else {
491 ifrstm = *ilo;
492 ilastm = *ihi;
493 }
494 iiter = 0;
495 eshift = 0.;
496 maxit = (*ihi - *ilo + 1) * 30;
497
498 i__1 = maxit;
499 for (jiter = 1; jiter <= i__1; ++jiter) {
500
501/* Split the matrix if possible.
502
503 Two tests:
504 1: A(j,j-1)=0 or j=ILO
505 2: B(j,j)=0 */
506
507 if (ilast == *ilo) {
508
509/* Special case: j=ILAST */
510
511 goto L80;
512 } else {
513 if ((d__1 = a_ref(ilast, ilast - 1), absMACRO(d__1)) <= atol) {
514 a_ref(ilast, ilast - 1) = 0.;
515 goto L80;
516 }
517 }
518
519 if ((d__1 = b_ref(ilast, ilast), absMACRO(d__1)) <= btol) {
520 b_ref(ilast, ilast) = 0.;
521 goto L70;
522 }
523
524/* General case: j<ILAST */
525
526 i__2 = *ilo;
527 for (j = ilast - 1; j >= i__2; --j) {
528
529/* Test 1: for A(j,j-1)=0 or j=ILO */
530
531 if (j == *ilo) {
532 ilazro = TRUE_;
533 } else {
534 if ((d__1 = a_ref(j, j - 1), absMACRO(d__1)) <= atol) {
535 a_ref(j, j - 1) = 0.;
536 ilazro = TRUE_;
537 } else {
538 ilazro = FALSE_;
539 }
540 }
541
542/* Test 2: for B(j,j)=0 */
543
544 if ((d__1 = b_ref(j, j), absMACRO(d__1)) < btol) {
545 b_ref(j, j) = 0.;
546
547/* Test 1a: Check for 2 consecutive small subdiagonals in A */
548
549 ilazr2 = FALSE_;
550 if (! ilazro) {
551 temp = (d__1 = a_ref(j, j - 1), absMACRO(d__1));
552 temp2 = (d__1 = a_ref(j, j), absMACRO(d__1));
553 tempr = maxMACRO(temp,temp2);
554 if (tempr < 1. && tempr != 0.) {
555 temp /= tempr;
556 temp2 /= tempr;
557 }
558 if (temp * (ascale * (d__1 = a_ref(j + 1, j), absMACRO(d__1)))
559 <= temp2 * (ascale * atol)) {
560 ilazr2 = TRUE_;
561 }
562 }
563
564/* If both tests pass (1 & 2), i.e., the leading diagonal
565 element of B in the block is zero, split a 1x1 block off
566 at the top. (I.e., at the J-th row/column) The leading
567 diagonal element of the remainder can also be zero, so
568 this may have to be done repeatedly. */
569
570 if (ilazro || ilazr2) {
571 i__3 = ilast - 1;
572 for (jch = j; jch <= i__3; ++jch) {
573 temp = a_ref(jch, jch);
574 template_lapack_lartg(&temp, &a_ref(jch + 1, jch), &c__, &s, &a_ref(
575 jch, jch));
576 a_ref(jch + 1, jch) = 0.;
577 i__4 = ilastm - jch;
578 template_blas_rot(&i__4, &a_ref(jch, jch + 1), lda, &a_ref(jch +
579 1, jch + 1), lda, &c__, &s);
580 i__4 = ilastm - jch;
581 template_blas_rot(&i__4, &b_ref(jch, jch + 1), ldb, &b_ref(jch +
582 1, jch + 1), ldb, &c__, &s);
583 if (ilq) {
584 template_blas_rot(n, &q_ref(1, jch), &c__1, &q_ref(1, jch + 1)
585 , &c__1, &c__, &s);
586 }
587 if (ilazr2) {
588 a_ref(jch, jch - 1) = a_ref(jch, jch - 1) * c__;
589 }
590 ilazr2 = FALSE_;
591 if ((d__1 = b_ref(jch + 1, jch + 1), absMACRO(d__1)) >=
592 btol) {
593 if (jch + 1 >= ilast) {
594 goto L80;
595 } else {
596 ifirst = jch + 1;
597 goto L110;
598 }
599 }
600 b_ref(jch + 1, jch + 1) = 0.;
601/* L40: */
602 }
603 goto L70;
604 } else {
605
606/* Only test 2 passed -- chase the zero to B(ILAST,ILAST)
607 Then process as in the case B(ILAST,ILAST)=0 */
608
609 i__3 = ilast - 1;
610 for (jch = j; jch <= i__3; ++jch) {
611 temp = b_ref(jch, jch + 1);
612 template_lapack_lartg(&temp, &b_ref(jch + 1, jch + 1), &c__, &s, &
613 b_ref(jch, jch + 1));
614 b_ref(jch + 1, jch + 1) = 0.;
615 if (jch < ilastm - 1) {
616 i__4 = ilastm - jch - 1;
617 template_blas_rot(&i__4, &b_ref(jch, jch + 2), ldb, &b_ref(
618 jch + 1, jch + 2), ldb, &c__, &s);
619 }
620 i__4 = ilastm - jch + 2;
621 template_blas_rot(&i__4, &a_ref(jch, jch - 1), lda, &a_ref(jch +
622 1, jch - 1), lda, &c__, &s);
623 if (ilq) {
624 template_blas_rot(n, &q_ref(1, jch), &c__1, &q_ref(1, jch + 1)
625 , &c__1, &c__, &s);
626 }
627 temp = a_ref(jch + 1, jch);
628 template_lapack_lartg(&temp, &a_ref(jch + 1, jch - 1), &c__, &s, &
629 a_ref(jch + 1, jch));
630 a_ref(jch + 1, jch - 1) = 0.;
631 i__4 = jch + 1 - ifrstm;
632 template_blas_rot(&i__4, &a_ref(ifrstm, jch), &c__1, &a_ref(
633 ifrstm, jch - 1), &c__1, &c__, &s);
634 i__4 = jch - ifrstm;
635 template_blas_rot(&i__4, &b_ref(ifrstm, jch), &c__1, &b_ref(
636 ifrstm, jch - 1), &c__1, &c__, &s);
637 if (ilz) {
638 template_blas_rot(n, &z___ref(1, jch), &c__1, &z___ref(1, jch
639 - 1), &c__1, &c__, &s);
640 }
641/* L50: */
642 }
643 goto L70;
644 }
645 } else if (ilazro) {
646
647/* Only test 1 passed -- work on J:ILAST */
648
649 ifirst = j;
650 goto L110;
651 }
652
653/* Neither test passed -- try next J
654
655 L60: */
656 }
657
658/* (Drop-through is "impossible") */
659
660 *info = *n + 1;
661 goto L420;
662
663/* B(ILAST,ILAST)=0 -- clear A(ILAST,ILAST-1) to split off a
664 1x1 block. */
665
666L70:
667 temp = a_ref(ilast, ilast);
668 template_lapack_lartg(&temp, &a_ref(ilast, ilast - 1), &c__, &s, &a_ref(ilast,
669 ilast));
670 a_ref(ilast, ilast - 1) = 0.;
671 i__2 = ilast - ifrstm;
672 template_blas_rot(&i__2, &a_ref(ifrstm, ilast), &c__1, &a_ref(ifrstm, ilast - 1),
673 &c__1, &c__, &s);
674 i__2 = ilast - ifrstm;
675 template_blas_rot(&i__2, &b_ref(ifrstm, ilast), &c__1, &b_ref(ifrstm, ilast - 1),
676 &c__1, &c__, &s);
677 if (ilz) {
678 template_blas_rot(n, &z___ref(1, ilast), &c__1, &z___ref(1, ilast - 1), &c__1,
679 &c__, &s);
680 }
681
682/* A(ILAST,ILAST-1)=0 -- Standardize B, set ALPHAR, ALPHAI,
683 and BETA */
684
685L80:
686 if (b_ref(ilast, ilast) < 0.) {
687 if (ilschr) {
688 i__2 = ilast;
689 for (j = ifrstm; j <= i__2; ++j) {
690 a_ref(j, ilast) = -a_ref(j, ilast);
691 b_ref(j, ilast) = -b_ref(j, ilast);
692/* L90: */
693 }
694 } else {
695 a_ref(ilast, ilast) = -a_ref(ilast, ilast);
696 b_ref(ilast, ilast) = -b_ref(ilast, ilast);
697 }
698 if (ilz) {
699 i__2 = *n;
700 for (j = 1; j <= i__2; ++j) {
701 z___ref(j, ilast) = -z___ref(j, ilast);
702/* L100: */
703 }
704 }
705 }
706 alphar[ilast] = a_ref(ilast, ilast);
707 alphai[ilast] = 0.;
708 beta[ilast] = b_ref(ilast, ilast);
709
710/* Go to next block -- exit if finished. */
711
712 --ilast;
713 if (ilast < *ilo) {
714 goto L380;
715 }
716
717/* Reset counters */
718
719 iiter = 0;
720 eshift = 0.;
721 if (! ilschr) {
722 ilastm = ilast;
723 if (ifrstm > ilast) {
724 ifrstm = *ilo;
725 }
726 }
727 goto L350;
728
729/* QZ step
730
731 This iteration only involves rows/columns IFIRST:ILAST. We
732 assume IFIRST < ILAST, and that the diagonal of B is non-zero. */
733
734L110:
735 ++iiter;
736 if (! ilschr) {
737 ifrstm = ifirst;
738 }
739
740/* Compute single shifts.
741
742 At this point, IFIRST < ILAST, and the diagonal elements of
743 B(IFIRST:ILAST,IFIRST,ILAST) are larger than BTOL (in
744 magnitude) */
745
746 if (iiter / 10 * 10 == iiter) {
747
748/* Exceptional shift. Chosen for no particularly good reason.
749 (Single shift only.) */
750
751 if ((Treal) maxit * safmin * (d__1 = a_ref(ilast - 1, ilast),
752 absMACRO(d__1)) < (d__2 = b_ref(ilast - 1, ilast - 1), absMACRO(
753 d__2))) {
754 eshift += a_ref(ilast - 1, ilast) / b_ref(ilast - 1, ilast -
755 1);
756 } else {
757 eshift += 1. / (safmin * (Treal) maxit);
758 }
759 s1 = 1.;
760 wr = eshift;
761
762 } else {
763
764/* Shifts based on the generalized eigenvalues of the
765 bottom-right 2x2 block of A and B. The first eigenvalue
766 returned by DLAG2 is the Wilkinson shift (AEP p.512), */
767
768 d__1 = safmin * 100.;
769 template_lapack_lag2(&a_ref(ilast - 1, ilast - 1), lda, &b_ref(ilast - 1, ilast
770 - 1), ldb, &d__1, &s1, &s2, &wr, &wr2, &wi);
771
772/* Computing MAX
773 Computing MAX */
774 d__3 = 1., d__4 = absMACRO(wr), d__3 = maxMACRO(d__3,d__4), d__4 = absMACRO(wi);
775 d__1 = s1, d__2 = safmin * maxMACRO(d__3,d__4);
776 temp = maxMACRO(d__1,d__2);
777 if (wi != 0.) {
778 goto L200;
779 }
780 }
781
782/* Fiddle with shift to avoid overflow */
783
784 temp = minMACRO(ascale,1.) * (safmax * .5);
785 if (s1 > temp) {
786 scale = temp / s1;
787 } else {
788 scale = 1.;
789 }
790
791 temp = minMACRO(bscale,1.) * (safmax * .5);
792 if (absMACRO(wr) > temp) {
793/* Computing MIN */
794 d__1 = scale, d__2 = temp / absMACRO(wr);
795 scale = minMACRO(d__1,d__2);
796 }
797 s1 = scale * s1;
798 wr = scale * wr;
799
800/* Now check for two consecutive small subdiagonals. */
801
802 i__2 = ifirst + 1;
803 for (j = ilast - 1; j >= i__2; --j) {
804 istart = j;
805 temp = (d__1 = s1 * a_ref(j, j - 1), absMACRO(d__1));
806 temp2 = (d__1 = s1 * a_ref(j, j) - wr * b_ref(j, j), absMACRO(d__1));
807 tempr = maxMACRO(temp,temp2);
808 if (tempr < 1. && tempr != 0.) {
809 temp /= tempr;
810 temp2 /= tempr;
811 }
812 if ((d__1 = ascale * a_ref(j + 1, j) * temp, absMACRO(d__1)) <= ascale
813 * atol * temp2) {
814 goto L130;
815 }
816/* L120: */
817 }
818
819 istart = ifirst;
820L130:
821
822/* Do an implicit single-shift QZ sweep.
823
824 Initial Q */
825
826 temp = s1 * a_ref(istart, istart) - wr * b_ref(istart, istart);
827 temp2 = s1 * a_ref(istart + 1, istart);
828 template_lapack_lartg(&temp, &temp2, &c__, &s, &tempr);
829
830/* Sweep */
831
832 i__2 = ilast - 1;
833 for (j = istart; j <= i__2; ++j) {
834 if (j > istart) {
835 temp = a_ref(j, j - 1);
836 template_lapack_lartg(&temp, &a_ref(j + 1, j - 1), &c__, &s, &a_ref(j, j -
837 1));
838 a_ref(j + 1, j - 1) = 0.;
839 }
840
841 i__3 = ilastm;
842 for (jc = j; jc <= i__3; ++jc) {
843 temp = c__ * a_ref(j, jc) + s * a_ref(j + 1, jc);
844 a_ref(j + 1, jc) = -s * a_ref(j, jc) + c__ * a_ref(j + 1, jc);
845 a_ref(j, jc) = temp;
846 temp2 = c__ * b_ref(j, jc) + s * b_ref(j + 1, jc);
847 b_ref(j + 1, jc) = -s * b_ref(j, jc) + c__ * b_ref(j + 1, jc);
848 b_ref(j, jc) = temp2;
849/* L140: */
850 }
851 if (ilq) {
852 i__3 = *n;
853 for (jr = 1; jr <= i__3; ++jr) {
854 temp = c__ * q_ref(jr, j) + s * q_ref(jr, j + 1);
855 q_ref(jr, j + 1) = -s * q_ref(jr, j) + c__ * q_ref(jr, j
856 + 1);
857 q_ref(jr, j) = temp;
858/* L150: */
859 }
860 }
861
862 temp = b_ref(j + 1, j + 1);
863 template_lapack_lartg(&temp, &b_ref(j + 1, j), &c__, &s, &b_ref(j + 1, j + 1));
864 b_ref(j + 1, j) = 0.;
865
866/* Computing MIN */
867 i__4 = j + 2;
868 i__3 = minMACRO(i__4,ilast);
869 for (jr = ifrstm; jr <= i__3; ++jr) {
870 temp = c__ * a_ref(jr, j + 1) + s * a_ref(jr, j);
871 a_ref(jr, j) = -s * a_ref(jr, j + 1) + c__ * a_ref(jr, j);
872 a_ref(jr, j + 1) = temp;
873/* L160: */
874 }
875 i__3 = j;
876 for (jr = ifrstm; jr <= i__3; ++jr) {
877 temp = c__ * b_ref(jr, j + 1) + s * b_ref(jr, j);
878 b_ref(jr, j) = -s * b_ref(jr, j + 1) + c__ * b_ref(jr, j);
879 b_ref(jr, j + 1) = temp;
880/* L170: */
881 }
882 if (ilz) {
883 i__3 = *n;
884 for (jr = 1; jr <= i__3; ++jr) {
885 temp = c__ * z___ref(jr, j + 1) + s * z___ref(jr, j);
886 z___ref(jr, j) = -s * z___ref(jr, j + 1) + c__ * z___ref(
887 jr, j);
888 z___ref(jr, j + 1) = temp;
889/* L180: */
890 }
891 }
892/* L190: */
893 }
894
895 goto L350;
896
897/* Use Francis double-shift
898
899 Note: the Francis double-shift should work with real shifts,
900 but only if the block is at least 3x3.
901 This code may break if this point is reached with
902 a 2x2 block with real eigenvalues. */
903
904L200:
905 if (ifirst + 1 == ilast) {
906
907/* Special case -- 2x2 block with complex eigenvectors
908
909 Step 1: Standardize, that is, rotate so that
910
911 ( B11 0 )
912 B = ( ) with B11 non-negative.
913 ( 0 B22 ) */
914
915 template_lapack_lasv2(&b_ref(ilast - 1, ilast - 1), &b_ref(ilast - 1, ilast), &
916 b_ref(ilast, ilast), &b22, &b11, &sr, &cr, &sl, &cl);
917
918 if (b11 < 0.) {
919 cr = -cr;
920 sr = -sr;
921 b11 = -b11;
922 b22 = -b22;
923 }
924
925 i__2 = ilastm + 1 - ifirst;
926 template_blas_rot(&i__2, &a_ref(ilast - 1, ilast - 1), lda, &a_ref(ilast,
927 ilast - 1), lda, &cl, &sl);
928 i__2 = ilast + 1 - ifrstm;
929 template_blas_rot(&i__2, &a_ref(ifrstm, ilast - 1), &c__1, &a_ref(ifrstm,
930 ilast), &c__1, &cr, &sr);
931
932 if (ilast < ilastm) {
933 i__2 = ilastm - ilast;
934 template_blas_rot(&i__2, &b_ref(ilast - 1, ilast + 1), ldb, &b_ref(ilast,
935 ilast + 1), lda, &cl, &sl);
936 }
937 if (ifrstm < ilast - 1) {
938 i__2 = ifirst - ifrstm;
939 template_blas_rot(&i__2, &b_ref(ifrstm, ilast - 1), &c__1, &b_ref(ifrstm,
940 ilast), &c__1, &cr, &sr);
941 }
942
943 if (ilq) {
944 template_blas_rot(n, &q_ref(1, ilast - 1), &c__1, &q_ref(1, ilast), &c__1,
945 &cl, &sl);
946 }
947 if (ilz) {
948 template_blas_rot(n, &z___ref(1, ilast - 1), &c__1, &z___ref(1, ilast), &
949 c__1, &cr, &sr);
950 }
951
952 b_ref(ilast - 1, ilast - 1) = b11;
953 b_ref(ilast - 1, ilast) = 0.;
954 b_ref(ilast, ilast - 1) = 0.;
955 b_ref(ilast, ilast) = b22;
956
957/* If B22 is negative, negate column ILAST */
958
959 if (b22 < 0.) {
960 i__2 = ilast;
961 for (j = ifrstm; j <= i__2; ++j) {
962 a_ref(j, ilast) = -a_ref(j, ilast);
963 b_ref(j, ilast) = -b_ref(j, ilast);
964/* L210: */
965 }
966
967 if (ilz) {
968 i__2 = *n;
969 for (j = 1; j <= i__2; ++j) {
970 z___ref(j, ilast) = -z___ref(j, ilast);
971/* L220: */
972 }
973 }
974 }
975
976/* Step 2: Compute ALPHAR, ALPHAI, and BETA (see refs.)
977
978 Recompute shift */
979
980 d__1 = safmin * 100.;
981 template_lapack_lag2(&a_ref(ilast - 1, ilast - 1), lda, &b_ref(ilast - 1, ilast
982 - 1), ldb, &d__1, &s1, &temp, &wr, &temp2, &wi);
983
984/* If standardization has perturbed the shift onto real line,
985 do another (real single-shift) QR step. */
986
987 if (wi == 0.) {
988 goto L350;
989 }
990 s1inv = 1. / s1;
991
992/* Do EISPACK (QZVAL) computation of alpha and beta */
993
994 a11 = a_ref(ilast - 1, ilast - 1);
995 a21 = a_ref(ilast, ilast - 1);
996 a12 = a_ref(ilast - 1, ilast);
997 a22 = a_ref(ilast, ilast);
998
999/* Compute complex Givens rotation on right
1000 (Assume some element of C = (sA - wB) > unfl )
1001 __
1002 (sA - wB) ( CZ -SZ )
1003 ( SZ CZ ) */
1004
1005 c11r = s1 * a11 - wr * b11;
1006 c11i = -wi * b11;
1007 c12 = s1 * a12;
1008 c21 = s1 * a21;
1009 c22r = s1 * a22 - wr * b22;
1010 c22i = -wi * b22;
1011
1012 if (absMACRO(c11r) + absMACRO(c11i) + absMACRO(c12) > absMACRO(c21) + absMACRO(c22r) + absMACRO(
1013 c22i)) {
1014 t = template_lapack_lapy3(&c12, &c11r, &c11i);
1015 cz = c12 / t;
1016 szr = -c11r / t;
1017 szi = -c11i / t;
1018 } else {
1019 cz = template_lapack_lapy2(&c22r, &c22i);
1020 if (cz <= safmin) {
1021 cz = 0.;
1022 szr = 1.;
1023 szi = 0.;
1024 } else {
1025 tempr = c22r / cz;
1026 tempi = c22i / cz;
1027 t = template_lapack_lapy2(&cz, &c21);
1028 cz /= t;
1029 szr = -c21 * tempr / t;
1030 szi = c21 * tempi / t;
1031 }
1032 }
1033
1034/* Compute Givens rotation on left
1035
1036 ( CQ SQ )
1037 ( __ ) A or B
1038 ( -SQ CQ ) */
1039
1040 an = absMACRO(a11) + absMACRO(a12) + absMACRO(a21) + absMACRO(a22);
1041 bn = absMACRO(b11) + absMACRO(b22);
1042 wabs = absMACRO(wr) + absMACRO(wi);
1043 if (s1 * an > wabs * bn) {
1044 cq = cz * b11;
1045 sqr = szr * b22;
1046 sqi = -szi * b22;
1047 } else {
1048 a1r = cz * a11 + szr * a12;
1049 a1i = szi * a12;
1050 a2r = cz * a21 + szr * a22;
1051 a2i = szi * a22;
1052 cq = template_lapack_lapy2(&a1r, &a1i);
1053 if (cq <= safmin) {
1054 cq = 0.;
1055 sqr = 1.;
1056 sqi = 0.;
1057 } else {
1058 tempr = a1r / cq;
1059 tempi = a1i / cq;
1060 sqr = tempr * a2r + tempi * a2i;
1061 sqi = tempi * a2r - tempr * a2i;
1062 }
1063 }
1064 t = template_lapack_lapy3(&cq, &sqr, &sqi);
1065 cq /= t;
1066 sqr /= t;
1067 sqi /= t;
1068
1069/* Compute diagonal elements of QBZ */
1070
1071 tempr = sqr * szr - sqi * szi;
1072 tempi = sqr * szi + sqi * szr;
1073 b1r = cq * cz * b11 + tempr * b22;
1074 b1i = tempi * b22;
1075 b1a = template_lapack_lapy2(&b1r, &b1i);
1076 b2r = cq * cz * b22 + tempr * b11;
1077 b2i = -tempi * b11;
1078 b2a = template_lapack_lapy2(&b2r, &b2i);
1079
1080/* Normalize so beta > 0, and Im( alpha1 ) > 0 */
1081
1082 beta[ilast - 1] = b1a;
1083 beta[ilast] = b2a;
1084 alphar[ilast - 1] = wr * b1a * s1inv;
1085 alphai[ilast - 1] = wi * b1a * s1inv;
1086 alphar[ilast] = wr * b2a * s1inv;
1087 alphai[ilast] = -(wi * b2a) * s1inv;
1088
1089/* Step 3: Go to next block -- exit if finished. */
1090
1091 ilast = ifirst - 1;
1092 if (ilast < *ilo) {
1093 goto L380;
1094 }
1095
1096/* Reset counters */
1097
1098 iiter = 0;
1099 eshift = 0.;
1100 if (! ilschr) {
1101 ilastm = ilast;
1102 if (ifrstm > ilast) {
1103 ifrstm = *ilo;
1104 }
1105 }
1106 goto L350;
1107 } else {
1108
1109/* Usual case: 3x3 or larger block, using Francis implicit
1110 double-shift
1111
1112 2
1113 Eigenvalue equation is w - c w + d = 0,
1114
1115 -1 2 -1
1116 so compute 1st column of (A B ) - c A B + d
1117 using the formula in QZIT (from EISPACK)
1118
1119 We assume that the block is at least 3x3 */
1120
1121 ad11 = ascale * a_ref(ilast - 1, ilast - 1) / (bscale * b_ref(
1122 ilast - 1, ilast - 1));
1123 ad21 = ascale * a_ref(ilast, ilast - 1) / (bscale * b_ref(ilast -
1124 1, ilast - 1));
1125 ad12 = ascale * a_ref(ilast - 1, ilast) / (bscale * b_ref(ilast,
1126 ilast));
1127 ad22 = ascale * a_ref(ilast, ilast) / (bscale * b_ref(ilast,
1128 ilast));
1129 u12 = b_ref(ilast - 1, ilast) / b_ref(ilast, ilast);
1130 ad11l = ascale * a_ref(ifirst, ifirst) / (bscale * b_ref(ifirst,
1131 ifirst));
1132 ad21l = ascale * a_ref(ifirst + 1, ifirst) / (bscale * b_ref(
1133 ifirst, ifirst));
1134 ad12l = ascale * a_ref(ifirst, ifirst + 1) / (bscale * b_ref(
1135 ifirst + 1, ifirst + 1));
1136 ad22l = ascale * a_ref(ifirst + 1, ifirst + 1) / (bscale * b_ref(
1137 ifirst + 1, ifirst + 1));
1138 ad32l = ascale * a_ref(ifirst + 2, ifirst + 1) / (bscale * b_ref(
1139 ifirst + 1, ifirst + 1));
1140 u12l = b_ref(ifirst, ifirst + 1) / b_ref(ifirst + 1, ifirst + 1);
1141
1142 v[0] = (ad11 - ad11l) * (ad22 - ad11l) - ad12 * ad21 + ad21 * u12
1143 * ad11l + (ad12l - ad11l * u12l) * ad21l;
1144 v[1] = (ad22l - ad11l - ad21l * u12l - (ad11 - ad11l) - (ad22 -
1145 ad11l) + ad21 * u12) * ad21l;
1146 v[2] = ad32l * ad21l;
1147
1148 istart = ifirst;
1149
1150 template_lapack_larfg(&c__3, v, &v[1], &c__1, &tau);
1151 v[0] = 1.;
1152
1153/* Sweep */
1154
1155 i__2 = ilast - 2;
1156 for (j = istart; j <= i__2; ++j) {
1157
1158/* All but last elements: use 3x3 Householder transforms.
1159
1160 Zero (j-1)st column of A */
1161
1162 if (j > istart) {
1163 v[0] = a_ref(j, j - 1);
1164 v[1] = a_ref(j + 1, j - 1);
1165 v[2] = a_ref(j + 2, j - 1);
1166
1167 template_lapack_larfg(&c__3, &a_ref(j, j - 1), &v[1], &c__1, &tau);
1168 v[0] = 1.;
1169 a_ref(j + 1, j - 1) = 0.;
1170 a_ref(j + 2, j - 1) = 0.;
1171 }
1172
1173 i__3 = ilastm;
1174 for (jc = j; jc <= i__3; ++jc) {
1175 temp = tau * (a_ref(j, jc) + v[1] * a_ref(j + 1, jc) + v[
1176 2] * a_ref(j + 2, jc));
1177 a_ref(j, jc) = a_ref(j, jc) - temp;
1178 a_ref(j + 1, jc) = a_ref(j + 1, jc) - temp * v[1];
1179 a_ref(j + 2, jc) = a_ref(j + 2, jc) - temp * v[2];
1180 temp2 = tau * (b_ref(j, jc) + v[1] * b_ref(j + 1, jc) + v[
1181 2] * b_ref(j + 2, jc));
1182 b_ref(j, jc) = b_ref(j, jc) - temp2;
1183 b_ref(j + 1, jc) = b_ref(j + 1, jc) - temp2 * v[1];
1184 b_ref(j + 2, jc) = b_ref(j + 2, jc) - temp2 * v[2];
1185/* L230: */
1186 }
1187 if (ilq) {
1188 i__3 = *n;
1189 for (jr = 1; jr <= i__3; ++jr) {
1190 temp = tau * (q_ref(jr, j) + v[1] * q_ref(jr, j + 1)
1191 + v[2] * q_ref(jr, j + 2));
1192 q_ref(jr, j) = q_ref(jr, j) - temp;
1193 q_ref(jr, j + 1) = q_ref(jr, j + 1) - temp * v[1];
1194 q_ref(jr, j + 2) = q_ref(jr, j + 2) - temp * v[2];
1195/* L240: */
1196 }
1197 }
1198
1199/* Zero j-th column of B (see DLAGBC for details)
1200
1201 Swap rows to pivot */
1202
1203 ilpivt = FALSE_;
1204/* Computing MAX */
1205 d__3 = (d__1 = b_ref(j + 1, j + 1), absMACRO(d__1)), d__4 = (d__2 =
1206 b_ref(j + 1, j + 2), absMACRO(d__2));
1207 temp = maxMACRO(d__3,d__4);
1208/* Computing MAX */
1209 d__3 = (d__1 = b_ref(j + 2, j + 1), absMACRO(d__1)), d__4 = (d__2 =
1210 b_ref(j + 2, j + 2), absMACRO(d__2));
1211 temp2 = maxMACRO(d__3,d__4);
1212 if (maxMACRO(temp,temp2) < safmin) {
1213 scale = 0.;
1214 u1 = 1.;
1215 u2 = 0.;
1216 goto L250;
1217 } else if (temp >= temp2) {
1218 w11 = b_ref(j + 1, j + 1);
1219 w21 = b_ref(j + 2, j + 1);
1220 w12 = b_ref(j + 1, j + 2);
1221 w22 = b_ref(j + 2, j + 2);
1222 u1 = b_ref(j + 1, j);
1223 u2 = b_ref(j + 2, j);
1224 } else {
1225 w21 = b_ref(j + 1, j + 1);
1226 w11 = b_ref(j + 2, j + 1);
1227 w22 = b_ref(j + 1, j + 2);
1228 w12 = b_ref(j + 2, j + 2);
1229 u2 = b_ref(j + 1, j);
1230 u1 = b_ref(j + 2, j);
1231 }
1232
1233/* Swap columns if nec. */
1234
1235 if (absMACRO(w12) > absMACRO(w11)) {
1236 ilpivt = TRUE_;
1237 temp = w12;
1238 temp2 = w22;
1239 w12 = w11;
1240 w22 = w21;
1241 w11 = temp;
1242 w21 = temp2;
1243 }
1244
1245/* LU-factor */
1246
1247 temp = w21 / w11;
1248 u2 -= temp * u1;
1249 w22 -= temp * w12;
1250 w21 = 0.;
1251
1252/* Compute SCALE */
1253
1254 scale = 1.;
1255 if (absMACRO(w22) < safmin) {
1256 scale = 0.;
1257 u2 = 1.;
1258 u1 = -w12 / w11;
1259 goto L250;
1260 }
1261 if (absMACRO(w22) < absMACRO(u2)) {
1262 scale = (d__1 = w22 / u2, absMACRO(d__1));
1263 }
1264 if (absMACRO(w11) < absMACRO(u1)) {
1265/* Computing MIN */
1266 d__2 = scale, d__3 = (d__1 = w11 / u1, absMACRO(d__1));
1267 scale = minMACRO(d__2,d__3);
1268 }
1269
1270/* Solve */
1271
1272 u2 = scale * u2 / w22;
1273 u1 = (scale * u1 - w12 * u2) / w11;
1274
1275L250:
1276 if (ilpivt) {
1277 temp = u2;
1278 u2 = u1;
1279 u1 = temp;
1280 }
1281
1282/* Compute Householder Vector
1283
1284 Computing 2nd power */
1285 d__1 = scale;
1286/* Computing 2nd power */
1287 d__2 = u1;
1288/* Computing 2nd power */
1289 d__3 = u2;
1290 t = template_blas_sqrt(d__1 * d__1 + d__2 * d__2 + d__3 * d__3);
1291 tau = scale / t + 1.;
1292 vs = -1. / (scale + t);
1293 v[0] = 1.;
1294 v[1] = vs * u1;
1295 v[2] = vs * u2;
1296
1297/* Apply transformations from the right.
1298
1299 Computing MIN */
1300 i__4 = j + 3;
1301 i__3 = minMACRO(i__4,ilast);
1302 for (jr = ifrstm; jr <= i__3; ++jr) {
1303 temp = tau * (a_ref(jr, j) + v[1] * a_ref(jr, j + 1) + v[
1304 2] * a_ref(jr, j + 2));
1305 a_ref(jr, j) = a_ref(jr, j) - temp;
1306 a_ref(jr, j + 1) = a_ref(jr, j + 1) - temp * v[1];
1307 a_ref(jr, j + 2) = a_ref(jr, j + 2) - temp * v[2];
1308/* L260: */
1309 }
1310 i__3 = j + 2;
1311 for (jr = ifrstm; jr <= i__3; ++jr) {
1312 temp = tau * (b_ref(jr, j) + v[1] * b_ref(jr, j + 1) + v[
1313 2] * b_ref(jr, j + 2));
1314 b_ref(jr, j) = b_ref(jr, j) - temp;
1315 b_ref(jr, j + 1) = b_ref(jr, j + 1) - temp * v[1];
1316 b_ref(jr, j + 2) = b_ref(jr, j + 2) - temp * v[2];
1317/* L270: */
1318 }
1319 if (ilz) {
1320 i__3 = *n;
1321 for (jr = 1; jr <= i__3; ++jr) {
1322 temp = tau * (z___ref(jr, j) + v[1] * z___ref(jr, j +
1323 1) + v[2] * z___ref(jr, j + 2));
1324 z___ref(jr, j) = z___ref(jr, j) - temp;
1325 z___ref(jr, j + 1) = z___ref(jr, j + 1) - temp * v[1];
1326 z___ref(jr, j + 2) = z___ref(jr, j + 2) - temp * v[2];
1327/* L280: */
1328 }
1329 }
1330 b_ref(j + 1, j) = 0.;
1331 b_ref(j + 2, j) = 0.;
1332/* L290: */
1333 }
1334
1335/* Last elements: Use Givens rotations
1336
1337 Rotations from the left */
1338
1339 j = ilast - 1;
1340 temp = a_ref(j, j - 1);
1341 template_lapack_lartg(&temp, &a_ref(j + 1, j - 1), &c__, &s, &a_ref(j, j - 1));
1342 a_ref(j + 1, j - 1) = 0.;
1343
1344 i__2 = ilastm;
1345 for (jc = j; jc <= i__2; ++jc) {
1346 temp = c__ * a_ref(j, jc) + s * a_ref(j + 1, jc);
1347 a_ref(j + 1, jc) = -s * a_ref(j, jc) + c__ * a_ref(j + 1, jc);
1348 a_ref(j, jc) = temp;
1349 temp2 = c__ * b_ref(j, jc) + s * b_ref(j + 1, jc);
1350 b_ref(j + 1, jc) = -s * b_ref(j, jc) + c__ * b_ref(j + 1, jc);
1351 b_ref(j, jc) = temp2;
1352/* L300: */
1353 }
1354 if (ilq) {
1355 i__2 = *n;
1356 for (jr = 1; jr <= i__2; ++jr) {
1357 temp = c__ * q_ref(jr, j) + s * q_ref(jr, j + 1);
1358 q_ref(jr, j + 1) = -s * q_ref(jr, j) + c__ * q_ref(jr, j
1359 + 1);
1360 q_ref(jr, j) = temp;
1361/* L310: */
1362 }
1363 }
1364
1365/* Rotations from the right. */
1366
1367 temp = b_ref(j + 1, j + 1);
1368 template_lapack_lartg(&temp, &b_ref(j + 1, j), &c__, &s, &b_ref(j + 1, j + 1));
1369 b_ref(j + 1, j) = 0.;
1370
1371 i__2 = ilast;
1372 for (jr = ifrstm; jr <= i__2; ++jr) {
1373 temp = c__ * a_ref(jr, j + 1) + s * a_ref(jr, j);
1374 a_ref(jr, j) = -s * a_ref(jr, j + 1) + c__ * a_ref(jr, j);
1375 a_ref(jr, j + 1) = temp;
1376/* L320: */
1377 }
1378 i__2 = ilast - 1;
1379 for (jr = ifrstm; jr <= i__2; ++jr) {
1380 temp = c__ * b_ref(jr, j + 1) + s * b_ref(jr, j);
1381 b_ref(jr, j) = -s * b_ref(jr, j + 1) + c__ * b_ref(jr, j);
1382 b_ref(jr, j + 1) = temp;
1383/* L330: */
1384 }
1385 if (ilz) {
1386 i__2 = *n;
1387 for (jr = 1; jr <= i__2; ++jr) {
1388 temp = c__ * z___ref(jr, j + 1) + s * z___ref(jr, j);
1389 z___ref(jr, j) = -s * z___ref(jr, j + 1) + c__ * z___ref(
1390 jr, j);
1391 z___ref(jr, j + 1) = temp;
1392/* L340: */
1393 }
1394 }
1395
1396/* End of Double-Shift code */
1397
1398 }
1399
1400 goto L350;
1401
1402/* End of iteration loop */
1403
1404L350:
1405/* L360: */
1406 ;
1407 }
1408
1409/* Drop-through = non-convergence
1410
1411 L370: */
1412 *info = ilast;
1413 goto L420;
1414
1415/* Successful completion of all QZ steps */
1416
1417L380:
1418
1419/* Set Eigenvalues 1:ILO-1 */
1420
1421 i__1 = *ilo - 1;
1422 for (j = 1; j <= i__1; ++j) {
1423 if (b_ref(j, j) < 0.) {
1424 if (ilschr) {
1425 i__2 = j;
1426 for (jr = 1; jr <= i__2; ++jr) {
1427 a_ref(jr, j) = -a_ref(jr, j);
1428 b_ref(jr, j) = -b_ref(jr, j);
1429/* L390: */
1430 }
1431 } else {
1432 a_ref(j, j) = -a_ref(j, j);
1433 b_ref(j, j) = -b_ref(j, j);
1434 }
1435 if (ilz) {
1436 i__2 = *n;
1437 for (jr = 1; jr <= i__2; ++jr) {
1438 z___ref(jr, j) = -z___ref(jr, j);
1439/* L400: */
1440 }
1441 }
1442 }
1443 alphar[j] = a_ref(j, j);
1444 alphai[j] = 0.;
1445 beta[j] = b_ref(j, j);
1446/* L410: */
1447 }
1448
1449/* Normal Termination */
1450
1451 *info = 0;
1452
1453/* Exit (other than argument error) -- return optimal workspace size */
1454
1455L420:
1456 work[1] = (Treal) (*n);
1457 return 0;
1458
1459/* End of DHGEQZ */
1460
1461} /* dhgeqz_ */
1462
1463#undef z___ref
1464#undef q_ref
1465#undef b_ref
1466#undef a_ref
1467
1468
1469#endif
Treal template_blas_sqrt(Treal x)
int template_blas_erbla(const char *srname, integer *info)
Definition template_blas_common.cc:146
logical template_blas_lsame(const char *ca, const char *cb)
Definition template_blas_common.cc:46
int integer
Definition template_blas_common.h:40
#define absMACRO(x)
Definition template_blas_common.h:47
#define minMACRO(a, b)
Definition template_blas_common.h:46
#define maxMACRO(a, b)
Definition template_blas_common.h:45
bool logical
Definition template_blas_common.h:41
int template_blas_rot(const integer *n, Treal *dx, const integer *incx, Treal *dy, const integer *incy, const Treal *c__, const Treal *s)
Definition template_blas_rot.h:42
#define TRUE_
Definition template_lapack_common.h:42
#define FALSE_
Definition template_lapack_common.h:43
int template_lapack_hgeqz(const char *job, const char *compq, const char *compz, const integer *n, const integer *ilo, const integer *ihi, Treal *a, const integer *lda, Treal *b, const integer *ldb, Treal *alphar, Treal *alphai, Treal *beta, Treal *q, const integer *ldq, Treal *z__, const integer *ldz, Treal *work, const integer *lwork, integer *info)
Definition template_lapack_hgeqz.h:42
#define a_ref(a_1, a_2)
#define b_ref(a_1, a_2)
#define z___ref(a_1, a_2)
#define q_ref(a_1, a_2)
int template_lapack_lag2(const Treal *a, const integer *lda, const Treal *b, const integer *ldb, const Treal *safmin, Treal *scale1, Treal *scale2, Treal *wr1, Treal *wr2, Treal *wi)
Definition template_lapack_lag2.h:42
#define cr
Treal template_lapack_lamch(const char *cmach, Treal dummyReal)
Definition template_lapack_lamch.h:202
Treal dlanhs_(const char *norm, const integer *n, const Treal *a, const integer *lda, Treal *work)
Definition template_lapack_lanhs.h:42
Treal template_lapack_lapy2(Treal *x, Treal *y)
Definition template_lapack_lapy2.h:42
Treal template_lapack_lapy3(Treal *x, Treal *y, Treal *z__)
Definition template_lapack_lapy3.h:42
int template_lapack_larfg(const integer *n, Treal *alpha, Treal *x, const integer *incx, Treal *tau)
Definition template_lapack_larfg.h:43
int template_lapack_lartg(const Treal *f, const Treal *g, Treal *cs, Treal *sn, Treal *r__)
Definition template_lapack_lartg.h:42
int template_lapack_laset(const char *uplo, const integer *m, const integer *n, const Treal *alpha, const Treal *beta, Treal *a, const integer *lda)
Definition template_lapack_laset.h:42
int template_lapack_lasv2(const Treal *f, const Treal *g, const Treal *h__, Treal *ssmin, Treal *ssmax, Treal *snr, Treal *csr, Treal *snl, Treal *csl)
Definition template_lapack_lasv2.h:42