ergo
template_lapack_hgeqz.h
Go to the documentation of this file.
1 /* Ergo, version 3.7, a program for linear scaling electronic structure
2  * calculations.
3  * Copyright (C) 2018 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 
41 template<class Treal>
42 int 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 
666 L70:
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 
685 L80:
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 
734 L110:
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;
820 L130:
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 
904 L200:
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 
1275 L250:
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 
1404 L350:
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 
1417 L380:
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 
1455 L420:
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
#define absMACRO(x)
Definition: template_blas_common.h:47
#define b_ref(a_1, a_2)
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
Treal template_lapack_lapy3(Treal *x, Treal *y, Treal *z__)
Definition: template_lapack_lapy3.h:42
int integer
Definition: template_blas_common.h:40
Treal template_lapack_lapy2(Treal *x, Treal *y)
Definition: template_lapack_lapy2.h:42
#define maxMACRO(a, b)
Definition: template_blas_common.h:45
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 minMACRO(a, b)
Definition: template_blas_common.h:46
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 z___ref(a_1, a_2)
int template_blas_erbla(const char *srname, integer *info)
Definition: template_blas_common.cc:146
#define a_ref(a_1, a_2)
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_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
int template_lapack_larfg(const integer *n, Treal *alpha, Treal *x, const integer *incx, Treal *tau)
Definition: template_lapack_larfg.h:43
Treal template_lapack_lamch(const char *cmach, Treal dummyReal)
Definition: template_lapack_lamch.h:202
bool logical
Definition: template_blas_common.h:41
#define TRUE_
Definition: template_lapack_common.h:42
#define FALSE_
Definition: template_lapack_common.h:43
#define q_ref(a_1, a_2)
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
Treal template_blas_sqrt(Treal x)
logical template_blas_lsame(const char *ca, const char *cb)
Definition: template_blas_common.cc:46
Treal dlanhs_(const char *norm, const integer *n, const Treal *a, const integer *lda, Treal *work)
Definition: template_lapack_lanhs.h:42
#define cr