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