1 /* dlahqr.f -- translated by f2c (version 20061008).
2 You must link the resulting object file with libf2c:
3 on Microsoft Windows system, link with libf2c.lib;
4 on Linux or Unix systems, link with .../path/to/libf2c.a -lm
5 or, if you install libf2c.a in a standard place, with -lf2c -lm
6 -- in that order, at the end of the command line, as in
7 cc *.o -lf2c -lm
8 Source for libf2c is in /netlib/f2c/libf2c.zip, e.g.,
9
10 http://www.netlib.org/f2c/libf2c.zip
11 */
12
13 #include "f2c.h"
14 #include "blaswrap.h"
15
16 /* Table of constant values */
17
18 static integer c__1 = 1;
19
dlahqr_(logical * wantt,logical * wantz,integer * n,integer * ilo,integer * ihi,doublereal * h__,integer * ldh,doublereal * wr,doublereal * wi,integer * iloz,integer * ihiz,doublereal * z__,integer * ldz,integer * info)20 /* Subroutine */ int dlahqr_(logical *wantt, logical *wantz, integer *n,
21 integer *ilo, integer *ihi, doublereal *h__, integer *ldh, doublereal
22 *wr, doublereal *wi, integer *iloz, integer *ihiz, doublereal *z__,
23 integer *ldz, integer *info)
24 {
25 /* System generated locals */
26 integer h_dim1, h_offset, z_dim1, z_offset, i__1, i__2, i__3;
27 doublereal d__1, d__2, d__3, d__4;
28
29 /* Builtin functions */
30 double sqrt(doublereal);
31
32 /* Local variables */
33 integer i__, j, k, l, m;
34 doublereal s, v[3];
35 integer i1, i2;
36 doublereal t1, t2, t3, v2, v3, aa, ab, ba, bb, h11, h12, h21, h22, cs;
37 integer nh;
38 doublereal sn;
39 integer nr;
40 doublereal tr;
41 integer nz;
42 doublereal det, h21s;
43 integer its;
44 doublereal ulp, sum, tst, rt1i, rt2i, rt1r, rt2r;
45 extern /* Subroutine */ int drot_(integer *, doublereal *, integer *,
46 doublereal *, integer *, doublereal *, doublereal *), dcopy_(
47 integer *, doublereal *, integer *, doublereal *, integer *),
48 dlanv2_(doublereal *, doublereal *, doublereal *, doublereal *,
49 doublereal *, doublereal *, doublereal *, doublereal *,
50 doublereal *, doublereal *), dlabad_(doublereal *, doublereal *);
51 extern doublereal dlamch_(char *);
52 extern /* Subroutine */ int dlarfg_(integer *, doublereal *, doublereal *,
53 integer *, doublereal *);
54 doublereal safmin, safmax, rtdisc, smlnum;
55
56
57 /* -- LAPACK auxiliary routine (version 3.2) -- */
58 /* Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd.. */
59 /* November 2006 */
60
61 /* .. Scalar Arguments .. */
62 /* .. */
63 /* .. Array Arguments .. */
64 /* .. */
65
66 /* Purpose */
67 /* ======= */
68
69 /* DLAHQR is an auxiliary routine called by DHSEQR to update the */
70 /* eigenvalues and Schur decomposition already computed by DHSEQR, by */
71 /* dealing with the Hessenberg submatrix in rows and columns ILO to */
72 /* IHI. */
73
74 /* Arguments */
75 /* ========= */
76
77 /* WANTT (input) LOGICAL */
78 /* = .TRUE. : the full Schur form T is required; */
79 /* = .FALSE.: only eigenvalues are required. */
80
81 /* WANTZ (input) LOGICAL */
82 /* = .TRUE. : the matrix of Schur vectors Z is required; */
83 /* = .FALSE.: Schur vectors are not required. */
84
85 /* N (input) INTEGER */
86 /* The order of the matrix H. N >= 0. */
87
88 /* ILO (input) INTEGER */
89 /* IHI (input) INTEGER */
90 /* It is assumed that H is already upper quasi-triangular in */
91 /* rows and columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless */
92 /* ILO = 1). DLAHQR works primarily with the Hessenberg */
93 /* submatrix in rows and columns ILO to IHI, but applies */
94 /* transformations to all of H if WANTT is .TRUE.. */
95 /* 1 <= ILO <= max(1,IHI); IHI <= N. */
96
97 /* H (input/output) DOUBLE PRECISION array, dimension (LDH,N) */
98 /* On entry, the upper Hessenberg matrix H. */
99 /* On exit, if INFO is zero and if WANTT is .TRUE., H is upper */
100 /* quasi-triangular in rows and columns ILO:IHI, with any */
101 /* 2-by-2 diagonal blocks in standard form. If INFO is zero */
102 /* and WANTT is .FALSE., the contents of H are unspecified on */
103 /* exit. The output state of H if INFO is nonzero is given */
104 /* below under the description of INFO. */
105
106 /* LDH (input) INTEGER */
107 /* The leading dimension of the array H. LDH >= max(1,N). */
108
109 /* WR (output) DOUBLE PRECISION array, dimension (N) */
110 /* WI (output) DOUBLE PRECISION array, dimension (N) */
111 /* The real and imaginary parts, respectively, of the computed */
112 /* eigenvalues ILO to IHI are stored in the corresponding */
113 /* elements of WR and WI. If two eigenvalues are computed as a */
114 /* complex conjugate pair, they are stored in consecutive */
115 /* elements of WR and WI, say the i-th and (i+1)th, with */
116 /* WI(i) > 0 and WI(i+1) < 0. If WANTT is .TRUE., the */
117 /* eigenvalues are stored in the same order as on the diagonal */
118 /* of the Schur form returned in H, with WR(i) = H(i,i), and, if */
119 /* H(i:i+1,i:i+1) is a 2-by-2 diagonal block, */
120 /* WI(i) = sqrt(H(i+1,i)*H(i,i+1)) and WI(i+1) = -WI(i). */
121
122 /* ILOZ (input) INTEGER */
123 /* IHIZ (input) INTEGER */
124 /* Specify the rows of Z to which transformations must be */
125 /* applied if WANTZ is .TRUE.. */
126 /* 1 <= ILOZ <= ILO; IHI <= IHIZ <= N. */
127
128 /* Z (input/output) DOUBLE PRECISION array, dimension (LDZ,N) */
129 /* If WANTZ is .TRUE., on entry Z must contain the current */
130 /* matrix Z of transformations accumulated by DHSEQR, and on */
131 /* exit Z has been updated; transformations are applied only to */
132 /* the submatrix Z(ILOZ:IHIZ,ILO:IHI). */
133 /* If WANTZ is .FALSE., Z is not referenced. */
134
135 /* LDZ (input) INTEGER */
136 /* The leading dimension of the array Z. LDZ >= max(1,N). */
137
138 /* INFO (output) INTEGER */
139 /* = 0: successful exit */
140 /* .GT. 0: If INFO = i, DLAHQR failed to compute all the */
141 /* eigenvalues ILO to IHI in a total of 30 iterations */
142 /* per eigenvalue; elements i+1:ihi of WR and WI */
143 /* contain those eigenvalues which have been */
144 /* successfully computed. */
145
146 /* If INFO .GT. 0 and WANTT is .FALSE., then on exit, */
147 /* the remaining unconverged eigenvalues are the */
148 /* eigenvalues of the upper Hessenberg matrix rows */
149 /* and columns ILO thorugh INFO of the final, output */
150 /* value of H. */
151
152 /* If INFO .GT. 0 and WANTT is .TRUE., then on exit */
153 /* (*) (initial value of H)*U = U*(final value of H) */
154 /* where U is an orthognal matrix. The final */
155 /* value of H is upper Hessenberg and triangular in */
156 /* rows and columns INFO+1 through IHI. */
157
158 /* If INFO .GT. 0 and WANTZ is .TRUE., then on exit */
159 /* (final value of Z) = (initial value of Z)*U */
160 /* where U is the orthogonal matrix in (*) */
161 /* (regardless of the value of WANTT.) */
162
163 /* Further Details */
164 /* =============== */
165
166 /* 02-96 Based on modifications by */
167 /* David Day, Sandia National Laboratory, USA */
168
169 /* 12-04 Further modifications by */
170 /* Ralph Byers, University of Kansas, USA */
171 /* This is a modified version of DLAHQR from LAPACK version 3.0. */
172 /* It is (1) more robust against overflow and underflow and */
173 /* (2) adopts the more conservative Ahues & Tisseur stopping */
174 /* criterion (LAWN 122, 1997). */
175
176 /* ========================================================= */
177
178 /* .. Parameters .. */
179 /* .. */
180 /* .. Local Scalars .. */
181 /* .. */
182 /* .. Local Arrays .. */
183 /* .. */
184 /* .. External Functions .. */
185 /* .. */
186 /* .. External Subroutines .. */
187 /* .. */
188 /* .. Intrinsic Functions .. */
189 /* .. */
190 /* .. Executable Statements .. */
191
192 /* Parameter adjustments */
193 h_dim1 = *ldh;
194 h_offset = 1 + h_dim1;
195 h__ -= h_offset;
196 --wr;
197 --wi;
198 z_dim1 = *ldz;
199 z_offset = 1 + z_dim1;
200 z__ -= z_offset;
201
202 /* Function Body */
203 *info = 0;
204
205 /* Quick return if possible */
206
207 if (*n == 0) {
208 return 0;
209 }
210 if (*ilo == *ihi) {
211 wr[*ilo] = h__[*ilo + *ilo * h_dim1];
212 wi[*ilo] = 0.;
213 return 0;
214 }
215
216 /* ==== clear out the trash ==== */
217 i__1 = *ihi - 3;
218 for (j = *ilo; j <= i__1; ++j) {
219 h__[j + 2 + j * h_dim1] = 0.;
220 h__[j + 3 + j * h_dim1] = 0.;
221 /* L10: */
222 }
223 if (*ilo <= *ihi - 2) {
224 h__[*ihi + (*ihi - 2) * h_dim1] = 0.;
225 }
226
227 nh = *ihi - *ilo + 1;
228 nz = *ihiz - *iloz + 1;
229
230 /* Set machine-dependent constants for the stopping criterion. */
231
232 safmin = dlamch_("SAFE MINIMUM");
233 safmax = 1. / safmin;
234 dlabad_(&safmin, &safmax);
235 ulp = dlamch_("PRECISION");
236 smlnum = safmin * ((doublereal) nh / ulp);
237
238 /* I1 and I2 are the indices of the first row and last column of H */
239 /* to which transformations must be applied. If eigenvalues only are */
240 /* being computed, I1 and I2 are set inside the main loop. */
241
242 if (*wantt) {
243 i1 = 1;
244 i2 = *n;
245 }
246
247 /* The main loop begins here. I is the loop index and decreases from */
248 /* IHI to ILO in steps of 1 or 2. Each iteration of the loop works */
249 /* with the active submatrix in rows and columns L to I. */
250 /* Eigenvalues I+1 to IHI have already converged. Either L = ILO or */
251 /* H(L,L-1) is negligible so that the matrix splits. */
252
253 i__ = *ihi;
254 L20:
255 l = *ilo;
256 if (i__ < *ilo) {
257 goto L160;
258 }
259
260 /* Perform QR iterations on rows and columns ILO to I until a */
261 /* submatrix of order 1 or 2 splits off at the bottom because a */
262 /* subdiagonal element has become negligible. */
263
264 for (its = 0; its <= 30; ++its) {
265
266 /* Look for a single small subdiagonal element. */
267
268 i__1 = l + 1;
269 for (k = i__; k >= i__1; --k) {
270 if ((d__1 = h__[k + (k - 1) * h_dim1], abs(d__1)) <= smlnum) {
271 goto L40;
272 }
273 tst = (d__1 = h__[k - 1 + (k - 1) * h_dim1], abs(d__1)) + (d__2 =
274 h__[k + k * h_dim1], abs(d__2));
275 if (tst == 0.) {
276 if (k - 2 >= *ilo) {
277 tst += (d__1 = h__[k - 1 + (k - 2) * h_dim1], abs(d__1));
278 }
279 if (k + 1 <= *ihi) {
280 tst += (d__1 = h__[k + 1 + k * h_dim1], abs(d__1));
281 }
282 }
283 /* ==== The following is a conservative small subdiagonal */
284 /* . deflation criterion due to Ahues & Tisseur (LAWN 122, */
285 /* . 1997). It has better mathematical foundation and */
286 /* . improves accuracy in some cases. ==== */
287 if ((d__1 = h__[k + (k - 1) * h_dim1], abs(d__1)) <= ulp * tst) {
288 /* Computing MAX */
289 d__3 = (d__1 = h__[k + (k - 1) * h_dim1], abs(d__1)), d__4 = (
290 d__2 = h__[k - 1 + k * h_dim1], abs(d__2));
291 ab = max(d__3,d__4);
292 /* Computing MIN */
293 d__3 = (d__1 = h__[k + (k - 1) * h_dim1], abs(d__1)), d__4 = (
294 d__2 = h__[k - 1 + k * h_dim1], abs(d__2));
295 ba = min(d__3,d__4);
296 /* Computing MAX */
297 d__3 = (d__1 = h__[k + k * h_dim1], abs(d__1)), d__4 = (d__2 =
298 h__[k - 1 + (k - 1) * h_dim1] - h__[k + k * h_dim1],
299 abs(d__2));
300 aa = max(d__3,d__4);
301 /* Computing MIN */
302 d__3 = (d__1 = h__[k + k * h_dim1], abs(d__1)), d__4 = (d__2 =
303 h__[k - 1 + (k - 1) * h_dim1] - h__[k + k * h_dim1],
304 abs(d__2));
305 bb = min(d__3,d__4);
306 s = aa + ab;
307 /* Computing MAX */
308 d__1 = smlnum, d__2 = ulp * (bb * (aa / s));
309 if (ba * (ab / s) <= max(d__1,d__2)) {
310 goto L40;
311 }
312 }
313 /* L30: */
314 }
315 L40:
316 l = k;
317 if (l > *ilo) {
318
319 /* H(L,L-1) is negligible */
320
321 h__[l + (l - 1) * h_dim1] = 0.;
322 }
323
324 /* Exit from loop if a submatrix of order 1 or 2 has split off. */
325
326 if (l >= i__ - 1) {
327 goto L150;
328 }
329
330 /* Now the active submatrix is in rows and columns L to I. If */
331 /* eigenvalues only are being computed, only the active submatrix */
332 /* need be transformed. */
333
334 if (! (*wantt)) {
335 i1 = l;
336 i2 = i__;
337 }
338
339 if (its == 10) {
340
341 /* Exceptional shift. */
342
343 s = (d__1 = h__[l + 1 + l * h_dim1], abs(d__1)) + (d__2 = h__[l +
344 2 + (l + 1) * h_dim1], abs(d__2));
345 h11 = s * .75 + h__[l + l * h_dim1];
346 h12 = s * -.4375;
347 h21 = s;
348 h22 = h11;
349 } else if (its == 20) {
350
351 /* Exceptional shift. */
352
353 s = (d__1 = h__[i__ + (i__ - 1) * h_dim1], abs(d__1)) + (d__2 =
354 h__[i__ - 1 + (i__ - 2) * h_dim1], abs(d__2));
355 h11 = s * .75 + h__[i__ + i__ * h_dim1];
356 h12 = s * -.4375;
357 h21 = s;
358 h22 = h11;
359 } else {
360
361 /* Prepare to use Francis' double shift */
362 /* (i.e. 2nd degree generalized Rayleigh quotient) */
363
364 h11 = h__[i__ - 1 + (i__ - 1) * h_dim1];
365 h21 = h__[i__ + (i__ - 1) * h_dim1];
366 h12 = h__[i__ - 1 + i__ * h_dim1];
367 h22 = h__[i__ + i__ * h_dim1];
368 }
369 s = abs(h11) + abs(h12) + abs(h21) + abs(h22);
370 if (s == 0.) {
371 rt1r = 0.;
372 rt1i = 0.;
373 rt2r = 0.;
374 rt2i = 0.;
375 } else {
376 h11 /= s;
377 h21 /= s;
378 h12 /= s;
379 h22 /= s;
380 tr = (h11 + h22) / 2.;
381 det = (h11 - tr) * (h22 - tr) - h12 * h21;
382 rtdisc = sqrt((abs(det)));
383 if (det >= 0.) {
384
385 /* ==== complex conjugate shifts ==== */
386
387 rt1r = tr * s;
388 rt2r = rt1r;
389 rt1i = rtdisc * s;
390 rt2i = -rt1i;
391 } else {
392
393 /* ==== real shifts (use only one of them) ==== */
394
395 rt1r = tr + rtdisc;
396 rt2r = tr - rtdisc;
397 if ((d__1 = rt1r - h22, abs(d__1)) <= (d__2 = rt2r - h22, abs(
398 d__2))) {
399 rt1r *= s;
400 rt2r = rt1r;
401 } else {
402 rt2r *= s;
403 rt1r = rt2r;
404 }
405 rt1i = 0.;
406 rt2i = 0.;
407 }
408 }
409
410 /* Look for two consecutive small subdiagonal elements. */
411
412 i__1 = l;
413 for (m = i__ - 2; m >= i__1; --m) {
414 /* Determine the effect of starting the double-shift QR */
415 /* iteration at row M, and see if this would make H(M,M-1) */
416 /* negligible. (The following uses scaling to avoid */
417 /* overflows and most underflows.) */
418
419 h21s = h__[m + 1 + m * h_dim1];
420 s = (d__1 = h__[m + m * h_dim1] - rt2r, abs(d__1)) + abs(rt2i) +
421 abs(h21s);
422 h21s = h__[m + 1 + m * h_dim1] / s;
423 v[0] = h21s * h__[m + (m + 1) * h_dim1] + (h__[m + m * h_dim1] -
424 rt1r) * ((h__[m + m * h_dim1] - rt2r) / s) - rt1i * (rt2i
425 / s);
426 v[1] = h21s * (h__[m + m * h_dim1] + h__[m + 1 + (m + 1) * h_dim1]
427 - rt1r - rt2r);
428 v[2] = h21s * h__[m + 2 + (m + 1) * h_dim1];
429 s = abs(v[0]) + abs(v[1]) + abs(v[2]);
430 v[0] /= s;
431 v[1] /= s;
432 v[2] /= s;
433 if (m == l) {
434 goto L60;
435 }
436 if ((d__1 = h__[m + (m - 1) * h_dim1], abs(d__1)) * (abs(v[1]) +
437 abs(v[2])) <= ulp * abs(v[0]) * ((d__2 = h__[m - 1 + (m -
438 1) * h_dim1], abs(d__2)) + (d__3 = h__[m + m * h_dim1],
439 abs(d__3)) + (d__4 = h__[m + 1 + (m + 1) * h_dim1], abs(
440 d__4)))) {
441 goto L60;
442 }
443 /* L50: */
444 }
445 L60:
446
447 /* Double-shift QR step */
448
449 i__1 = i__ - 1;
450 for (k = m; k <= i__1; ++k) {
451
452 /* The first iteration of this loop determines a reflection G */
453 /* from the vector V and applies it from left and right to H, */
454 /* thus creating a nonzero bulge below the subdiagonal. */
455
456 /* Each subsequent iteration determines a reflection G to */
457 /* restore the Hessenberg form in the (K-1)th column, and thus */
458 /* chases the bulge one step toward the bottom of the active */
459 /* submatrix. NR is the order of G. */
460
461 /* Computing MIN */
462 i__2 = 3, i__3 = i__ - k + 1;
463 nr = min(i__2,i__3);
464 if (k > m) {
465 dcopy_(&nr, &h__[k + (k - 1) * h_dim1], &c__1, v, &c__1);
466 }
467 dlarfg_(&nr, v, &v[1], &c__1, &t1);
468 if (k > m) {
469 h__[k + (k - 1) * h_dim1] = v[0];
470 h__[k + 1 + (k - 1) * h_dim1] = 0.;
471 if (k < i__ - 1) {
472 h__[k + 2 + (k - 1) * h_dim1] = 0.;
473 }
474 } else if (m > l) {
475 /* ==== Use the following instead of */
476 /* . H( K, K-1 ) = -H( K, K-1 ) to */
477 /* . avoid a bug when v(2) and v(3) */
478 /* . underflow. ==== */
479 h__[k + (k - 1) * h_dim1] *= 1. - t1;
480 }
481 v2 = v[1];
482 t2 = t1 * v2;
483 if (nr == 3) {
484 v3 = v[2];
485 t3 = t1 * v3;
486
487 /* Apply G from the left to transform the rows of the matrix */
488 /* in columns K to I2. */
489
490 i__2 = i2;
491 for (j = k; j <= i__2; ++j) {
492 sum = h__[k + j * h_dim1] + v2 * h__[k + 1 + j * h_dim1]
493 + v3 * h__[k + 2 + j * h_dim1];
494 h__[k + j * h_dim1] -= sum * t1;
495 h__[k + 1 + j * h_dim1] -= sum * t2;
496 h__[k + 2 + j * h_dim1] -= sum * t3;
497 /* L70: */
498 }
499
500 /* Apply G from the right to transform the columns of the */
501 /* matrix in rows I1 to min(K+3,I). */
502
503 /* Computing MIN */
504 i__3 = k + 3;
505 i__2 = min(i__3,i__);
506 for (j = i1; j <= i__2; ++j) {
507 sum = h__[j + k * h_dim1] + v2 * h__[j + (k + 1) * h_dim1]
508 + v3 * h__[j + (k + 2) * h_dim1];
509 h__[j + k * h_dim1] -= sum * t1;
510 h__[j + (k + 1) * h_dim1] -= sum * t2;
511 h__[j + (k + 2) * h_dim1] -= sum * t3;
512 /* L80: */
513 }
514
515 if (*wantz) {
516
517 /* Accumulate transformations in the matrix Z */
518
519 i__2 = *ihiz;
520 for (j = *iloz; j <= i__2; ++j) {
521 sum = z__[j + k * z_dim1] + v2 * z__[j + (k + 1) *
522 z_dim1] + v3 * z__[j + (k + 2) * z_dim1];
523 z__[j + k * z_dim1] -= sum * t1;
524 z__[j + (k + 1) * z_dim1] -= sum * t2;
525 z__[j + (k + 2) * z_dim1] -= sum * t3;
526 /* L90: */
527 }
528 }
529 } else if (nr == 2) {
530
531 /* Apply G from the left to transform the rows of the matrix */
532 /* in columns K to I2. */
533
534 i__2 = i2;
535 for (j = k; j <= i__2; ++j) {
536 sum = h__[k + j * h_dim1] + v2 * h__[k + 1 + j * h_dim1];
537 h__[k + j * h_dim1] -= sum * t1;
538 h__[k + 1 + j * h_dim1] -= sum * t2;
539 /* L100: */
540 }
541
542 /* Apply G from the right to transform the columns of the */
543 /* matrix in rows I1 to min(K+3,I). */
544
545 i__2 = i__;
546 for (j = i1; j <= i__2; ++j) {
547 sum = h__[j + k * h_dim1] + v2 * h__[j + (k + 1) * h_dim1]
548 ;
549 h__[j + k * h_dim1] -= sum * t1;
550 h__[j + (k + 1) * h_dim1] -= sum * t2;
551 /* L110: */
552 }
553
554 if (*wantz) {
555
556 /* Accumulate transformations in the matrix Z */
557
558 i__2 = *ihiz;
559 for (j = *iloz; j <= i__2; ++j) {
560 sum = z__[j + k * z_dim1] + v2 * z__[j + (k + 1) *
561 z_dim1];
562 z__[j + k * z_dim1] -= sum * t1;
563 z__[j + (k + 1) * z_dim1] -= sum * t2;
564 /* L120: */
565 }
566 }
567 }
568 /* L130: */
569 }
570
571 /* L140: */
572 }
573
574 /* Failure to converge in remaining number of iterations */
575
576 *info = i__;
577 return 0;
578
579 L150:
580
581 if (l == i__) {
582
583 /* H(I,I-1) is negligible: one eigenvalue has converged. */
584
585 wr[i__] = h__[i__ + i__ * h_dim1];
586 wi[i__] = 0.;
587 } else if (l == i__ - 1) {
588
589 /* H(I-1,I-2) is negligible: a pair of eigenvalues have converged. */
590
591 /* Transform the 2-by-2 submatrix to standard Schur form, */
592 /* and compute and store the eigenvalues. */
593
594 dlanv2_(&h__[i__ - 1 + (i__ - 1) * h_dim1], &h__[i__ - 1 + i__ *
595 h_dim1], &h__[i__ + (i__ - 1) * h_dim1], &h__[i__ + i__ *
596 h_dim1], &wr[i__ - 1], &wi[i__ - 1], &wr[i__], &wi[i__], &cs,
597 &sn);
598
599 if (*wantt) {
600
601 /* Apply the transformation to the rest of H. */
602
603 if (i2 > i__) {
604 i__1 = i2 - i__;
605 drot_(&i__1, &h__[i__ - 1 + (i__ + 1) * h_dim1], ldh, &h__[
606 i__ + (i__ + 1) * h_dim1], ldh, &cs, &sn);
607 }
608 i__1 = i__ - i1 - 1;
609 drot_(&i__1, &h__[i1 + (i__ - 1) * h_dim1], &c__1, &h__[i1 + i__ *
610 h_dim1], &c__1, &cs, &sn);
611 }
612 if (*wantz) {
613
614 /* Apply the transformation to Z. */
615
616 drot_(&nz, &z__[*iloz + (i__ - 1) * z_dim1], &c__1, &z__[*iloz +
617 i__ * z_dim1], &c__1, &cs, &sn);
618 }
619 }
620
621 /* return to start of the main loop with new value of I. */
622
623 i__ = l - 1;
624 goto L20;
625
626 L160:
627 return 0;
628
629 /* End of DLAHQR */
630
631 } /* dlahqr_ */
632