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