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