24
22
integer *ldz, integer *info)
26
24
/* System generated locals */
27
integer h_dim1, h_offset, z_dim1, z_offset, i__1, i__2, i__3, i__4;
28
doublereal d__1, d__2;
25
integer h_dim1, h_offset, z_dim1, z_offset, i__1, i__2, i__3;
26
doublereal d__1, d__2, d__3, d__4;
30
28
/* Builtin functions */
31
double sqrt(doublereal), igraphd_sign(doublereal *, doublereal *);
29
double sqrt(doublereal);
33
31
/* Local variables */
34
static integer i__, j, k, l, m;
35
static doublereal s, v[3];
36
static integer i1, i2;
37
static doublereal t1, t2, t3, v1, v2, v3, h00, h10, h11, h12, h21, h22,
44
static doublereal ave, h33s, h44s;
45
static integer itn, its;
46
static doublereal ulp, sum, tst1, h43h34, disc, unfl, ovfl;
32
integer i__, j, k, l, m;
35
doublereal t1, t2, t3, v2, v3, aa, ab, ba, bb, h11, h12, h21, h22, cs;
43
doublereal ulp, sum, tst, rt1i, rt2i, rt1r, rt2r;
47
44
extern /* Subroutine */ int igraphdrot_(integer *, doublereal *, integer *,
48
doublereal *, integer *, doublereal *, doublereal *);
49
static doublereal work[1];
50
extern /* Subroutine */ int igraphdcopy_(integer *, doublereal *, integer *,
51
doublereal *, integer *), igraphdlanv2_(doublereal *, doublereal *,
45
doublereal *, integer *, doublereal *, doublereal *), igraphdcopy_(
46
integer *, doublereal *, integer *, doublereal *, integer *),
47
igraphdlanv2_(doublereal *, doublereal *, doublereal *, doublereal *,
52
48
doublereal *, doublereal *, doublereal *, doublereal *,
53
doublereal *, doublereal *, doublereal *, doublereal *), igraphdlabad_(
54
doublereal *, doublereal *);
49
doublereal *, doublereal *), igraphdlabad_(doublereal *, doublereal *);
55
50
extern doublereal igraphdlamch_(char *);
56
51
extern /* Subroutine */ int igraphdlarfg_(integer *, doublereal *, doublereal *,
57
52
integer *, doublereal *);
58
extern doublereal igraphdlanhs_(char *, integer *, doublereal *, integer *,
60
static doublereal smlnum;
63
/* -- LAPACK auxiliary routine (version 3.0) -- */
64
/* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., */
65
/* Courant Institute, Argonne National Lab, and Rice University */
68
/* .. Scalar Arguments .. */
70
/* .. Array Arguments .. */
76
/* DLAHQR is an auxiliary routine called by DHSEQR to update the */
77
/* eigenvalues and Schur decomposition already computed by DHSEQR, by */
78
/* dealing with the Hessenberg submatrix in rows and columns ILO to IHI. */
83
/* WANTT (input) LOGICAL */
84
/* = .TRUE. : the full Schur form T is required; */
85
/* = .FALSE.: only eigenvalues are required. */
87
/* WANTZ (input) LOGICAL */
88
/* = .TRUE. : the matrix of Schur vectors Z is required; */
89
/* = .FALSE.: Schur vectors are not required. */
91
/* N (input) INTEGER */
92
/* The order of the matrix H. N >= 0. */
94
/* ILO (input) INTEGER */
95
/* IHI (input) INTEGER */
96
/* It is assumed that H is already upper quasi-triangular in */
97
/* rows and columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless */
98
/* ILO = 1). DLAHQR works primarily with the Hessenberg */
99
/* submatrix in rows and columns ILO to IHI, but applies */
100
/* transformations to all of H if WANTT is .TRUE.. */
101
/* 1 <= ILO <= max(1,IHI); IHI <= N. */
103
/* H (input/output) DOUBLE PRECISION array, dimension (LDH,N) */
104
/* On entry, the upper Hessenberg matrix H. */
105
/* On exit, if WANTT is .TRUE., H is upper quasi-triangular in */
106
/* rows and columns ILO:IHI, with any 2-by-2 diagonal blocks in */
107
/* standard form. If WANTT is .FALSE., the contents of H are */
108
/* unspecified on exit. */
110
/* LDH (input) INTEGER */
111
/* The leading dimension of the array H. LDH >= max(1,N). */
113
/* WR (output) DOUBLE PRECISION array, dimension (N) */
114
/* WI (output) DOUBLE PRECISION array, dimension (N) */
115
/* The real and imaginary parts, respectively, of the computed */
116
/* eigenvalues ILO to IHI are stored in the corresponding */
117
/* elements of WR and WI. If two eigenvalues are computed as a */
118
/* complex conjugate pair, they are stored in consecutive */
119
/* elements of WR and WI, say the i-th and (i+1)th, with */
120
/* WI(i) > 0 and WI(i+1) < 0. If WANTT is .TRUE., the */
121
/* eigenvalues are stored in the same order as on the diagonal */
122
/* of the Schur form returned in H, with WR(i) = H(i,i), and, if */
123
/* H(i:i+1,i:i+1) is a 2-by-2 diagonal block, */
124
/* WI(i) = sqrt(H(i+1,i)*H(i,i+1)) and WI(i+1) = -WI(i). */
126
/* ILOZ (input) INTEGER */
127
/* IHIZ (input) INTEGER */
128
/* Specify the rows of Z to which transformations must be */
129
/* applied if WANTZ is .TRUE.. */
130
/* 1 <= ILOZ <= ILO; IHI <= IHIZ <= N. */
132
/* Z (input/output) DOUBLE PRECISION array, dimension (LDZ,N) */
133
/* If WANTZ is .TRUE., on entry Z must contain the current */
134
/* matrix Z of transformations accumulated by DHSEQR, and on */
135
/* exit Z has been updated; transformations are applied only to */
136
/* the submatrix Z(ILOZ:IHIZ,ILO:IHI). */
137
/* If WANTZ is .FALSE., Z is not referenced. */
139
/* LDZ (input) INTEGER */
140
/* The leading dimension of the array Z. LDZ >= max(1,N). */
142
/* INFO (output) INTEGER */
143
/* = 0: successful exit */
144
/* > 0: DLAHQR failed to compute all the eigenvalues ILO to IHI */
145
/* in a total of 30*(IHI-ILO+1) iterations; if INFO = i, */
146
/* elements i+1:ihi of WR and WI contain those eigenvalues */
147
/* which have been successfully computed. */
149
/* Further Details */
150
/* =============== */
152
/* 2-96 Based on modifications by */
153
/* David Day, Sandia National Laboratory, USA */
155
/* ===================================================================== */
157
/* .. Parameters .. */
159
/* .. Local Scalars .. */
161
/* .. Local Arrays .. */
163
/* .. External Functions .. */
165
/* .. External Subroutines .. */
167
/* .. Intrinsic Functions .. */
169
/* .. Executable Statements .. */
171
/* Parameter adjustments */
53
doublereal safmin, safmax, rtdisc, smlnum;
56
/* -- LAPACK auxiliary routine (version 3.2) --
57
Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..
64
DLAHQR is an auxiliary routine called by DHSEQR to update the
65
eigenvalues and Schur decomposition already computed by DHSEQR, by
66
dealing with the Hessenberg submatrix in rows and columns ILO to
73
= .TRUE. : the full Schur form T is required;
74
= .FALSE.: only eigenvalues are required.
77
= .TRUE. : the matrix of Schur vectors Z is required;
78
= .FALSE.: Schur vectors are not required.
81
The order of the matrix H. N >= 0.
85
It is assumed that H is already upper quasi-triangular in
86
rows and columns IHI+1:N, and that H(ILO,ILO-1) = 0 (unless
87
ILO = 1). DLAHQR works primarily with the Hessenberg
88
submatrix in rows and columns ILO to IHI, but applies
89
transformations to all of H if WANTT is .TRUE..
90
1 <= ILO <= max(1,IHI); IHI <= N.
92
H (input/output) DOUBLE PRECISION array, dimension (LDH,N)
93
On entry, the upper Hessenberg matrix H.
94
On exit, if INFO is zero and if WANTT is .TRUE., H is upper
95
quasi-triangular in rows and columns ILO:IHI, with any
96
2-by-2 diagonal blocks in standard form. If INFO is zero
97
and WANTT is .FALSE., the contents of H are unspecified on
98
exit. The output state of H if INFO is nonzero is given
99
below under the description of INFO.
102
The leading dimension of the array H. LDH >= max(1,N).
104
WR (output) DOUBLE PRECISION array, dimension (N)
105
WI (output) DOUBLE PRECISION array, dimension (N)
106
The real and imaginary parts, respectively, of the computed
107
eigenvalues ILO to IHI are stored in the corresponding
108
elements of WR and WI. If two eigenvalues are computed as a
109
complex conjugate pair, they are stored in consecutive
110
elements of WR and WI, say the i-th and (i+1)th, with
111
WI(i) > 0 and WI(i+1) < 0. If WANTT is .TRUE., the
112
eigenvalues are stored in the same order as on the diagonal
113
of the Schur form returned in H, with WR(i) = H(i,i), and, if
114
H(i:i+1,i:i+1) is a 2-by-2 diagonal block,
115
WI(i) = sqrt(H(i+1,i)*H(i,i+1)) and WI(i+1) = -WI(i).
119
Specify the rows of Z to which transformations must be
120
applied if WANTZ is .TRUE..
121
1 <= ILOZ <= ILO; IHI <= IHIZ <= N.
123
Z (input/output) DOUBLE PRECISION array, dimension (LDZ,N)
124
If WANTZ is .TRUE., on entry Z must contain the current
125
matrix Z of transformations accumulated by DHSEQR, and on
126
exit Z has been updated; transformations are applied only to
127
the submatrix Z(ILOZ:IHIZ,ILO:IHI).
128
If WANTZ is .FALSE., Z is not referenced.
131
The leading dimension of the array Z. LDZ >= max(1,N).
133
INFO (output) INTEGER
135
.GT. 0: If INFO = i, DLAHQR failed to compute all the
136
eigenvalues ILO to IHI in a total of 30 iterations
137
per eigenvalue; elements i+1:ihi of WR and WI
138
contain those eigenvalues which have been
139
successfully computed.
141
If INFO .GT. 0 and WANTT is .FALSE., then on exit,
142
the remaining unconverged eigenvalues are the
143
eigenvalues of the upper Hessenberg matrix rows
144
and columns ILO thorugh INFO of the final, output
147
If INFO .GT. 0 and WANTT is .TRUE., then on exit
148
(*) (initial value of H)*U = U*(final value of H)
149
where U is an orthognal matrix. The final
150
value of H is upper Hessenberg and triangular in
151
rows and columns INFO+1 through IHI.
153
If INFO .GT. 0 and WANTZ is .TRUE., then on exit
154
(final value of Z) = (initial value of Z)*U
155
where U is the orthogonal matrix in (*)
156
(regardless of the value of WANTT.)
161
02-96 Based on modifications by
162
David Day, Sandia National Laboratory, USA
164
12-04 Further modifications by
165
Ralph Byers, University of Kansas, USA
166
This is a modified version of DLAHQR from LAPACK version 3.0.
167
It is (1) more robust against overflow and underflow and
168
(2) adopts the more conservative Ahues & Tisseur stopping
169
criterion (LAWN 122, 1997).
171
=========================================================
174
Parameter adjustments */
173
176
h_offset = 1 + h_dim1;
198
/* ==== clear out the trash ==== */
200
for (j = *ilo; j <= i__1; ++j) {
201
h__[j + 2 + j * h_dim1] = 0.;
202
h__[j + 3 + j * h_dim1] = 0.;
205
if (*ilo <= *ihi - 2) {
206
h__[*ihi + (*ihi - 2) * h_dim1] = 0.;
195
209
nh = *ihi - *ilo + 1;
196
210
nz = *ihiz - *iloz + 1;
198
212
/* Set machine-dependent constants for the stopping criterion. */
199
/* If norm(H) <= sqrt(OVFL), overflow should not occur. */
201
unfl = igraphdlamch_("Safe minimum");
203
igraphdlabad_(&unfl, &ovfl);
204
ulp = igraphdlamch_("Precision");
205
smlnum = unfl * (nh / ulp);
207
/* I1 and I2 are the indices of the first row and last column of H */
208
/* to which transformations must be applied. If eigenvalues only are */
209
/* being computed, I1 and I2 are set inside the main loop. */
214
safmin = igraphdlamch_("SAFE MINIMUM");
215
safmax = 1. / safmin;
216
igraphdlabad_(&safmin, &safmax);
217
ulp = igraphdlamch_("PRECISION");
218
smlnum = safmin * ((doublereal) nh / ulp);
220
/* I1 and I2 are the indices of the first row and last column of H
221
to which transformations must be applied. If eigenvalues only are
222
being computed, I1 and I2 are set inside the main loop. */
216
/* ITN is the total number of QR iterations allowed. */
220
/* The main loop begins here. I is the loop index and decreases from */
221
/* IHI to ILO in steps of 1 or 2. Each iteration of the loop works */
222
/* with the active submatrix in rows and columns L to I. */
223
/* Eigenvalues I+1 to IHI have already converged. Either L = ILO or */
224
/* H(L,L-1) is negligible so that the matrix splits. */
229
/* The main loop begins here. I is the loop index and decreases from
230
IHI to ILO in steps of 1 or 2. Each iteration of the loop works
231
with the active submatrix in rows and columns L to I.
232
Eigenvalues I+1 to IHI have already converged. Either L = ILO or
233
H(L,L-1) is negligible so that the matrix splits. */
229
238
if (i__ < *ilo) {
233
/* Perform QR iterations on rows and columns ILO to I until a */
234
/* submatrix of order 1 or 2 splits off at the bottom because a */
235
/* subdiagonal element has become negligible. */
242
/* Perform QR iterations on rows and columns ILO to I until a
243
submatrix of order 1 or 2 splits off at the bottom because a
244
subdiagonal element has become negligible. */
238
for (its = 0; its <= i__1; ++its) {
246
for (its = 0; its <= 30; ++its) {
240
248
/* Look for a single small subdiagonal element. */
243
for (k = i__; k >= i__2; --k) {
244
tst1 = (d__1 = h__[k - 1 + (k - 1) * h_dim1], abs(d__1)) + (d__2 =
245
h__[k + k * h_dim1], abs(d__2));
248
tst1 = igraphdlanhs_("1", &i__3, &h__[l + l * h_dim1], ldh, work);
252
if ((d__1 = h__[k + (k - 1) * h_dim1], abs(d__1)) <= max(d__2,
251
for (k = i__; k >= i__1; --k) {
252
if ((d__1 = h__[k + (k - 1) * h_dim1], abs(d__1)) <= smlnum) {
255
tst = (d__1 = h__[k - 1 + (k - 1) * h_dim1], abs(d__1)) + (d__2 =
256
h__[k + k * h_dim1], abs(d__2));
259
tst += (d__1 = h__[k - 1 + (k - 2) * h_dim1], abs(d__1));
262
tst += (d__1 = h__[k + 1 + k * h_dim1], abs(d__1));
265
/* ==== The following is a conservative small subdiagonal
266
. deflation criterion due to Ahues & Tisseur (LAWN 122,
267
. 1997). It has better mathematical foundation and
268
. improves accuracy in some cases. ==== */
269
if ((d__1 = h__[k + (k - 1) * h_dim1], abs(d__1)) <= ulp * tst) {
271
d__3 = (d__1 = h__[k + (k - 1) * h_dim1], abs(d__1)), d__4 = (
272
d__2 = h__[k - 1 + k * h_dim1], abs(d__2));
275
d__3 = (d__1 = h__[k + (k - 1) * h_dim1], abs(d__1)), d__4 = (
276
d__2 = h__[k - 1 + k * h_dim1], abs(d__2));
279
d__3 = (d__1 = h__[k + k * h_dim1], abs(d__1)), d__4 = (d__2 =
280
h__[k - 1 + (k - 1) * h_dim1] - h__[k + k * h_dim1],
284
d__3 = (d__1 = h__[k + k * h_dim1], abs(d__1)), d__4 = (d__2 =
285
h__[k - 1 + (k - 1) * h_dim1] - h__[k + k * h_dim1],
290
d__1 = smlnum, d__2 = ulp * (bb * (aa / s));
291
if (ba * (ab / s) <= max(d__1,d__2)) {
267
306
/* Exit from loop if a submatrix of order 1 or 2 has split off. */
269
308
if (l >= i__ - 1) {
273
/* Now the active submatrix is in rows and columns L to I. If */
274
/* eigenvalues only are being computed, only the active submatrix */
275
/* need be transformed. */
312
/* Now the active submatrix is in rows and columns L to I. If
313
eigenvalues only are being computed, only the active submatrix
314
need be transformed. */
277
316
if (! (*wantt)) {
282
if (its == 10 || its == 20) {
323
/* Exceptional shift. */
325
s = (d__1 = h__[l + 1 + l * h_dim1], abs(d__1)) + (d__2 = h__[l +
326
2 + (l + 1) * h_dim1], abs(d__2));
327
h11 = s * .75 + h__[l + l * h_dim1];
331
} else if (its == 20) {
284
333
/* Exceptional shift. */
286
335
s = (d__1 = h__[i__ + (i__ - 1) * h_dim1], abs(d__1)) + (d__2 =
287
336
h__[i__ - 1 + (i__ - 2) * h_dim1], abs(d__2));
288
h44 = s * .75 + h__[i__ + i__ * h_dim1];
290
h43h34 = s * -.4375 * s;
293
/* Prepare to use Francis' double shift */
294
/* (i.e. 2nd degree generalized Rayleigh quotient) */
296
h44 = h__[i__ + i__ * h_dim1];
297
h33 = h__[i__ - 1 + (i__ - 1) * h_dim1];
298
h43h34 = h__[i__ + (i__ - 1) * h_dim1] * h__[i__ - 1 + i__ *
300
s = h__[i__ - 1 + (i__ - 2) * h_dim1] * h__[i__ - 1 + (i__ - 2) *
302
disc = (h33 - h44) * .5;
303
disc = disc * disc + h43h34;
306
/* Real roots: use Wilkinson's shift twice */
309
ave = (h33 + h44) * .5;
310
if (abs(h33) - abs(h44) > 0.) {
311
h33 = h33 * h44 - h43h34;
312
h44 = h33 / (igraphd_sign(&disc, &ave) + ave);
337
h11 = s * .75 + h__[i__ + i__ * h_dim1];
343
/* Prepare to use Francis' double shift
344
(i.e. 2nd degree generalized Rayleigh quotient) */
346
h11 = h__[i__ - 1 + (i__ - 1) * h_dim1];
347
h21 = h__[i__ + (i__ - 1) * h_dim1];
348
h12 = h__[i__ - 1 + i__ * h_dim1];
349
h22 = h__[i__ + i__ * h_dim1];
351
s = abs(h11) + abs(h12) + abs(h21) + abs(h22);
362
tr = (h11 + h22) / 2.;
363
det = (h11 - tr) * (h22 - tr) - h12 * h21;
364
rtdisc = sqrt((abs(det)));
367
/* ==== complex conjugate shifts ==== */
375
/* ==== real shifts (use only one of them) ==== */
379
if ((d__1 = rt1r - h22, abs(d__1)) <= (d__2 = rt2r - h22, abs(
314
h44 = igraphd_sign(&disc, &ave) + ave;
321
392
/* Look for two consecutive small subdiagonal elements. */
324
for (m = i__ - 2; m >= i__2; --m) {
325
/* Determine the effect of starting the double-shift QR */
326
/* iteration at row M, and see if this would make H(M,M-1) */
395
for (m = i__ - 2; m >= i__1; --m) {
396
/* Determine the effect of starting the double-shift QR
397
iteration at row M, and see if this would make H(M,M-1)
398
negligible. (The following uses scaling to avoid
399
overflows and most underflows.) */
329
h11 = h__[m + m * h_dim1];
330
h22 = h__[m + 1 + (m + 1) * h_dim1];
331
h21 = h__[m + 1 + m * h_dim1];
332
h12 = h__[m + (m + 1) * h_dim1];
335
v1 = (h33s * h44s - h43h34) / h21 + h12;
336
v2 = h22 - h11 - h33s - h44s;
337
v3 = h__[m + 2 + (m + 1) * h_dim1];
338
s = abs(v1) + abs(v2) + abs(v3);
401
h21s = h__[m + 1 + m * h_dim1];
402
s = (d__1 = h__[m + m * h_dim1] - rt2r, abs(d__1)) + abs(rt2i) +
404
h21s = h__[m + 1 + m * h_dim1] / s;
405
v[0] = h21s * h__[m + (m + 1) * h_dim1] + (h__[m + m * h_dim1] -
406
rt1r) * ((h__[m + m * h_dim1] - rt2r) / s) - rt1i * (rt2i
408
v[1] = h21s * (h__[m + m * h_dim1] + h__[m + 1 + (m + 1) * h_dim1]
410
v[2] = h21s * h__[m + 2 + (m + 1) * h_dim1];
411
s = abs(v[0]) + abs(v[1]) + abs(v[2]);
348
h00 = h__[m - 1 + (m - 1) * h_dim1];
349
h10 = h__[m + (m - 1) * h_dim1];
350
tst1 = abs(v1) * (abs(h00) + abs(h11) + abs(h22));
351
if (abs(h10) * (abs(v2) + abs(v3)) <= ulp * tst1) {
418
if ((d__1 = h__[m + (m - 1) * h_dim1], abs(d__1)) * (abs(v[1]) +
419
abs(v[2])) <= ulp * abs(v[0]) * ((d__2 = h__[m - 1 + (m -
420
1) * h_dim1], abs(d__2)) + (d__3 = h__[m + m * h_dim1],
421
abs(d__3)) + (d__4 = h__[m + 1 + (m + 1) * h_dim1], abs(
358
429
/* Double-shift QR step */
361
for (k = m; k <= i__2; ++k) {
363
/* The first iteration of this loop determines a reflection G */
364
/* from the vector V and applies it from left and right to H, */
365
/* thus creating a nonzero bulge below the subdiagonal. */
367
/* Each subsequent iteration determines a reflection G to */
368
/* restore the Hessenberg form in the (K-1)th column, and thus */
369
/* chases the bulge one step toward the bottom of the active */
370
/* submatrix. NR is the order of G. */
373
i__3 = 3, i__4 = i__ - k + 1;
432
for (k = m; k <= i__1; ++k) {
434
/* The first iteration of this loop determines a reflection G
435
from the vector V and applies it from left and right to H,
436
thus creating a nonzero bulge below the subdiagonal.
438
Each subsequent iteration determines a reflection G to
439
restore the Hessenberg form in the (K-1)th column, and thus
440
chases the bulge one step toward the bottom of the active
441
submatrix. NR is the order of G.
444
i__2 = 3, i__3 = i__ - k + 1;
376
447
igraphdcopy_(&nr, &h__[k + (k - 1) * h_dim1], &c__1, v, &c__1);
394
/* Apply G from the left to transform the rows of the matrix */
395
/* in columns K to I2. */
469
/* Apply G from the left to transform the rows of the matrix
470
in columns K to I2. */
398
for (j = k; j <= i__3; ++j) {
473
for (j = k; j <= i__2; ++j) {
399
474
sum = h__[k + j * h_dim1] + v2 * h__[k + 1 + j * h_dim1]
400
475
+ v3 * h__[k + 2 + j * h_dim1];
401
476
h__[k + j * h_dim1] -= sum * t1;
402
477
h__[k + 1 + j * h_dim1] -= sum * t2;
403
478
h__[k + 2 + j * h_dim1] -= sum * t3;
407
/* Apply G from the right to transform the columns of the */
408
/* matrix in rows I1 to min(K+3,I). */
482
/* Apply G from the right to transform the columns of the
483
matrix in rows I1 to min(K+3,I).
412
i__3 = min(i__4,i__);
413
for (j = i1; j <= i__3; ++j) {
487
i__2 = min(i__3,i__);
488
for (j = i1; j <= i__2; ++j) {
414
489
sum = h__[j + k * h_dim1] + v2 * h__[j + (k + 1) * h_dim1]
415
490
+ v3 * h__[j + (k + 2) * h_dim1];
416
491
h__[j + k * h_dim1] -= sum * t1;
417
492
h__[j + (k + 1) * h_dim1] -= sum * t2;
418
493
h__[j + (k + 2) * h_dim1] -= sum * t3;
424
499
/* Accumulate transformations in the matrix Z */
427
for (j = *iloz; j <= i__3; ++j) {
502
for (j = *iloz; j <= i__2; ++j) {
428
503
sum = z__[j + k * z_dim1] + v2 * z__[j + (k + 1) *
429
504
z_dim1] + v3 * z__[j + (k + 2) * z_dim1];
430
505
z__[j + k * z_dim1] -= sum * t1;
431
506
z__[j + (k + 1) * z_dim1] -= sum * t2;
432
507
z__[j + (k + 2) * z_dim1] -= sum * t3;
436
511
} else if (nr == 2) {
438
/* Apply G from the left to transform the rows of the matrix */
439
/* in columns K to I2. */
513
/* Apply G from the left to transform the rows of the matrix
514
in columns K to I2. */
442
for (j = k; j <= i__3; ++j) {
517
for (j = k; j <= i__2; ++j) {
443
518
sum = h__[k + j * h_dim1] + v2 * h__[k + 1 + j * h_dim1];
444
519
h__[k + j * h_dim1] -= sum * t1;
445
520
h__[k + 1 + j * h_dim1] -= sum * t2;
449
/* Apply G from the right to transform the columns of the */
450
/* matrix in rows I1 to min(K+3,I). */
524
/* Apply G from the right to transform the columns of the
525
matrix in rows I1 to min(K+3,I). */
453
for (j = i1; j <= i__3; ++j) {
528
for (j = i1; j <= i__2; ++j) {
454
529
sum = h__[j + k * h_dim1] + v2 * h__[j + (k + 1) * h_dim1]
456
531
h__[j + k * h_dim1] -= sum * t1;
457
532
h__[j + (k + 1) * h_dim1] -= sum * t2;
463
538
/* Accumulate transformations in the matrix Z */
466
for (j = *iloz; j <= i__3; ++j) {
541
for (j = *iloz; j <= i__2; ++j) {
467
542
sum = z__[j + k * z_dim1] + v2 * z__[j + (k + 1) *
469
544
z__[j + k * z_dim1] -= sum * t1;
470
545
z__[j + (k + 1) * z_dim1] -= sum * t2;
481
556
/* Failure to converge in remaining number of iterations */
added
removed
src/lapack/dlahqr.c
