1 // Copyright ©2016 The Gonum Authors. All rights reserved.
2 // Use of this source code is governed by a BSD-style
3 // license that can be found in the LICENSE file.
10 "gonum.org/v1/gonum/blas/blas64"
13 // Dlahqr computes the eigenvalues and Schur factorization of a block of an n×n
14 // upper Hessenberg matrix H, using the double-shift/single-shift QR algorithm.
16 // h and ldh represent the matrix H. Dlahqr works primarily with the Hessenberg
17 // submatrix H[ilo:ihi+1,ilo:ihi+1], but applies transformations to all of H if
18 // wantt is true. It is assumed that H[ihi+1:n,ihi+1:n] is already upper
19 // quasi-triangular, although this is not checked.
22 // 0 <= ilo <= max(0,ihi), and ihi < n,
24 // H[ilo,ilo-1] == 0, if ilo > 0,
25 // otherwise Dlahqr will panic.
27 // If unconverged is zero on return, wr[ilo:ihi+1] and wi[ilo:ihi+1] will contain
28 // respectively the real and imaginary parts of the computed eigenvalues ilo
29 // to ihi. If two eigenvalues are computed as a complex conjugate pair, they are
30 // stored in consecutive elements of wr and wi, say the i-th and (i+1)th, with
31 // wi[i] > 0 and wi[i+1] < 0. If wantt is true, the eigenvalues are stored in
32 // the same order as on the diagonal of the Schur form returned in H, with
33 // wr[i] = H[i,i], and, if H[i:i+2,i:i+2] is a 2×2 diagonal block,
34 // wi[i] = sqrt(abs(H[i+1,i]*H[i,i+1])) and wi[i+1] = -wi[i].
36 // wr and wi must have length ihi+1.
38 // z and ldz represent an n×n matrix Z. If wantz is true, the transformations
39 // will be applied to the submatrix Z[iloz:ihiz+1,ilo:ihi+1] and it must hold that
40 // 0 <= iloz <= ilo, and ihi <= ihiz < n.
41 // If wantz is false, z is not referenced.
43 // unconverged indicates whether Dlahqr computed all the eigenvalues ilo to ihi
44 // in a total of 30 iterations per eigenvalue.
46 // If unconverged is zero, all the eigenvalues ilo to ihi have been computed and
47 // will be stored on return in wr[ilo:ihi+1] and wi[ilo:ihi+1].
49 // If unconverged is zero and wantt is true, H[ilo:ihi+1,ilo:ihi+1] will be
50 // overwritten on return by upper quasi-triangular full Schur form with any
51 // 2×2 diagonal blocks in standard form.
53 // If unconverged is zero and if wantt is false, the contents of h on return is
56 // If unconverged is positive, some eigenvalues have not converged, and
57 // wr[unconverged:ihi+1] and wi[unconverged:ihi+1] contain those eigenvalues
58 // which have been successfully computed.
60 // If unconverged is positive and wantt is true, then on return
61 // (initial H)*U = U*(final H), (*)
62 // where U is an orthogonal matrix. The final H is upper Hessenberg and
63 // H[unconverged:ihi+1,unconverged:ihi+1] is upper quasi-triangular.
65 // If unconverged is positive and wantt is false, on return the remaining
66 // unconverged eigenvalues are the eigenvalues of the upper Hessenberg matrix
67 // H[ilo:unconverged,ilo:unconverged].
69 // If unconverged is positive and wantz is true, then on return
70 // (final Z) = (initial Z)*U,
71 // where U is the orthogonal matrix in (*) regardless of the value of wantt.
73 // Dlahqr is an internal routine. It is exported for testing purposes.
74 func (impl Implementation) Dlahqr(wantt, wantz bool, n, ilo, ihi int, h []float64, ldh int, wr, wi []float64, iloz, ihiz int, z []float64, ldz int) (unconverged int) {
75 checkMatrix(n, n, h, ldh)
77 case ilo < 0 || max(0, ihi) < ilo:
81 case len(wr) != ihi+1:
82 panic("lapack: bad length of wr")
83 case len(wi) != ihi+1:
84 panic("lapack: bad length of wi")
85 case ilo > 0 && h[ilo*ldh+ilo-1] != 0:
86 panic("lapack: block is not isolated")
89 checkMatrix(n, n, z, ldz)
91 case iloz < 0 || ilo < iloz:
92 panic("lapack: iloz out of range")
93 case ihiz < ihi || n <= ihiz:
94 panic("lapack: ihiz out of range")
98 // Quick return if possible.
103 wr[ilo] = h[ilo*ldh+ilo]
108 // Clear out the trash.
109 for j := ilo; j < ihi-2; j++ {
118 nz := ihiz - iloz + 1
120 // Set machine-dependent constants for the stopping criterion.
122 smlnum := float64(nh) / ulp * dlamchS
124 // i1 and i2 are the indices of the first row and last column of H to
125 // which transformations must be applied. If eigenvalues only are being
126 // computed, i1 and i2 are set inside the main loop.
133 itmax := 30 * max(10, nh) // Total number of QR iterations allowed.
135 // The main loop begins here. i is the loop index and decreases from ihi
136 // to ilo in steps of 1 or 2. Each iteration of the loop works with the
137 // active submatrix in rows and columns l to i. Eigenvalues i+1 to ihi
138 // have already converged. Either l = ilo or H[l,l-1] is negligible so
139 // that the matrix splits.
140 bi := blas64.Implementation()
145 // Perform QR iterations on rows and columns ilo to i until a
146 // submatrix of order 1 or 2 splits off at the bottom because a
147 // subdiagonal element has become negligible.
149 for its := 0; its <= itmax; its++ {
150 // Look for a single small subdiagonal element.
152 for k = i; k > l; k-- {
153 if math.Abs(h[k*ldh+k-1]) <= smlnum {
156 tst := math.Abs(h[(k-1)*ldh+k-1]) + math.Abs(h[k*ldh+k])
159 tst += math.Abs(h[(k-1)*ldh+k-2])
162 tst += math.Abs(h[(k+1)*ldh+k])
165 // The following is a conservative small
166 // subdiagonal deflation criterion due to Ahues
167 // & Tisseur (LAWN 122, 1997). It has better
168 // mathematical foundation and improves accuracy
170 if math.Abs(h[k*ldh+k-1]) <= ulp*tst {
171 ab := math.Max(math.Abs(h[k*ldh+k-1]), math.Abs(h[(k-1)*ldh+k]))
172 ba := math.Min(math.Abs(h[k*ldh+k-1]), math.Abs(h[(k-1)*ldh+k]))
173 aa := math.Max(math.Abs(h[k*ldh+k]), math.Abs(h[(k-1)*ldh+k-1]-h[k*ldh+k]))
174 bb := math.Min(math.Abs(h[k*ldh+k]), math.Abs(h[(k-1)*ldh+k-1]-h[k*ldh+k]))
176 if ab/s*ba <= math.Max(smlnum, aa/s*bb*ulp) {
183 // H[l,l-1] is negligible.
187 // Break the loop because a submatrix of order 1
188 // or 2 has split off.
193 // Now the active submatrix is in rows and columns l to
194 // i. If eigenvalues only are being computed, only the
195 // active submatrix need be transformed.
205 var h11, h21, h12, h22 float64
207 case 10: // Exceptional shift.
208 s := math.Abs(h[(l+1)*ldh+l]) + math.Abs(h[(l+2)*ldh+l+1])
209 h11 = dat1*s + h[l*ldh+l]
213 case 20: // Exceptional shift.
214 s := math.Abs(h[i*ldh+i-1]) + math.Abs(h[(i-1)*ldh+i-2])
215 h11 = dat1*s + h[i*ldh+i]
219 default: // Prepare to use Francis' double shift (i.e.,
220 // 2nd degree generalized Rayleigh quotient).
221 h11 = h[(i-1)*ldh+i-1]
226 s := math.Abs(h11) + math.Abs(h12) + math.Abs(h21) + math.Abs(h22)
236 tr := (h11 + h22) / 2
237 det := (h11-tr)*(h22-tr) - h12*h21
238 rtdisc := math.Sqrt(math.Abs(det))
240 // Complex conjugate shifts.
246 // Real shifts (use only one of them).
249 if math.Abs(rt1r-h22) <= math.Abs(rt2r-h22) {
261 // Look for two consecutive small subdiagonal elements.
264 for m = i - 2; m >= l; m-- {
265 // Determine the effect of starting the
266 // double-shift QR iteration at row m, and see
267 // if this would make H[m,m-1] negligible. The
268 // following uses scaling to avoid overflows and
270 h21s := h[(m+1)*ldh+m]
271 s := math.Abs(h[m*ldh+m]-rt2r) + math.Abs(rt2i) + math.Abs(h21s)
273 v[0] = h21s*h[m*ldh+m+1] + (h[m*ldh+m]-rt1r)*((h[m*ldh+m]-rt2r)/s) - rt2i/s*rt1i
274 v[1] = h21s * (h[m*ldh+m] + h[(m+1)*ldh+m+1] - rt1r - rt2r)
275 v[2] = h21s * h[(m+2)*ldh+m+1]
276 s = math.Abs(v[0]) + math.Abs(v[1]) + math.Abs(v[2])
283 dsum := math.Abs(h[(m-1)*ldh+m-1]) + math.Abs(h[m*ldh+m]) + math.Abs(h[(m+1)*ldh+m+1])
284 if math.Abs(h[m*ldh+m-1])*(math.Abs(v[1])+math.Abs(v[2])) <= ulp*math.Abs(v[0])*dsum {
289 // Double-shift QR step.
290 for k := m; k < i; k++ {
291 // The first iteration of this loop determines a
292 // reflection G from the vector V and applies it
293 // from left and right to H, thus creating a
294 // non-zero bulge below the subdiagonal.
296 // Each subsequent iteration determines a
297 // reflection G to restore the Hessenberg form
298 // in the (k-1)th column, and thus chases the
299 // bulge one step toward the bottom of the
300 // active submatrix. nr is the order of G.
304 bi.Dcopy(nr, h[k*ldh+k-1:], ldh, v[:], 1)
307 v[0], t0 = impl.Dlarfg(nr, v[0], v[1:], 1)
315 // Use the following instead of H[k,k-1] = -H[k,k-1]
316 // to avoid a bug when v[1] and v[2] underflow.
317 h[k*ldh+k-1] *= 1 - t0
323 // Apply G from the left to transform
324 // the rows of the matrix in columns k
326 for j := k; j <= i2; j++ {
327 sum := h[k*ldh+j] + v[1]*h[(k+1)*ldh+j] + v[2]*h[(k+2)*ldh+j]
328 h[k*ldh+j] -= sum * t0
329 h[(k+1)*ldh+j] -= sum * t1
330 h[(k+2)*ldh+j] -= sum * t2
333 // Apply G from the right to transform
334 // the columns of the matrix in rows i1
336 for j := i1; j <= min(k+3, i); j++ {
337 sum := h[j*ldh+k] + v[1]*h[j*ldh+k+1] + v[2]*h[j*ldh+k+2]
338 h[j*ldh+k] -= sum * t0
339 h[j*ldh+k+1] -= sum * t1
340 h[j*ldh+k+2] -= sum * t2
344 // Accumulate transformations in the matrix Z.
345 for j := iloz; j <= ihiz; j++ {
346 sum := z[j*ldz+k] + v[1]*z[j*ldz+k+1] + v[2]*z[j*ldz+k+2]
347 z[j*ldz+k] -= sum * t0
348 z[j*ldz+k+1] -= sum * t1
349 z[j*ldz+k+2] -= sum * t2
353 // Apply G from the left to transform
354 // the rows of the matrix in columns k
356 for j := k; j <= i2; j++ {
357 sum := h[k*ldh+j] + v[1]*h[(k+1)*ldh+j]
358 h[k*ldh+j] -= sum * t0
359 h[(k+1)*ldh+j] -= sum * t1
362 // Apply G from the right to transform
363 // the columns of the matrix in rows i1
365 for j := i1; j <= i; j++ {
366 sum := h[j*ldh+k] + v[1]*h[j*ldh+k+1]
367 h[j*ldh+k] -= sum * t0
368 h[j*ldh+k+1] -= sum * t1
372 // Accumulate transformations in the matrix Z.
373 for j := iloz; j <= ihiz; j++ {
374 sum := z[j*ldz+k] + v[1]*z[j*ldz+k+1]
375 z[j*ldz+k] -= sum * t0
376 z[j*ldz+k+1] -= sum * t1
384 // The QR iteration finished without splitting off a
385 // submatrix of order 1 or 2.
390 // H[i,i-1] is negligible: one eigenvalue has converged.
394 // H[i-1,i-2] is negligible: a pair of eigenvalues have converged.
396 // Transform the 2×2 submatrix to standard Schur form,
397 // and compute and store the eigenvalues.
399 a, b := h[(i-1)*ldh+i-1], h[(i-1)*ldh+i]
400 c, d := h[i*ldh+i-1], h[i*ldh+i]
401 a, b, c, d, wr[i-1], wi[i-1], wr[i], wi[i], cs, sn = impl.Dlanv2(a, b, c, d)
402 h[(i-1)*ldh+i-1], h[(i-1)*ldh+i] = a, b
403 h[i*ldh+i-1], h[i*ldh+i] = c, d
406 // Apply the transformation to the rest of H.
408 bi.Drot(i2-i, h[(i-1)*ldh+i+1:], 1, h[i*ldh+i+1:], 1, cs, sn)
410 bi.Drot(i-i1-1, h[i1*ldh+i-1:], ldh, h[i1*ldh+i:], ldh, cs, sn)
414 // Apply the transformation to Z.
415 bi.Drot(nz, z[iloz*ldz+i-1:], ldz, z[iloz*ldz+i:], ldz, cs, sn)
419 // Return to start of the main loop with new value of i.