vendor/gonum.org/v1/gonum/lapack/gonum/dlahqr.go

   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.
   4
   5 package gonum
   6
   7 import (
   8         "math"
   9
  10         "gonum.org/v1/gonum/blas/blas64"
  11 )
  12
  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.
  15 //
  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.
  20 //
  21 // It must hold that
  22 //  0 <= ilo <= max(0,ihi), and ihi < n,
  23 // and that
  24 //  H[ilo,ilo-1] == 0,  if ilo > 0,
  25 // otherwise Dlahqr will panic.
  26 //
  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].
  35 //
  36 // wr and wi must have length ihi+1.
  37 //
  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.
  42 //
  43 // unconverged indicates whether Dlahqr computed all the eigenvalues ilo to ihi
  44 // in a total of 30 iterations per eigenvalue.
  45 //
  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].
  48 //
  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.
  52 //
  53 // If unconverged is zero and if wantt is false, the contents of h on return is
  54 // unspecified.
  55 //
  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.
  59 //
  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.
  64 //
  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].
  68 //
  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.
  72 //
  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)
  76         switch {
  77         case ilo < 0 || max(0, ihi) < ilo:
  78                 panic(badIlo)
  79         case n <= ihi:
  80                 panic(badIhi)
  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")
  87         }
  88         if wantz {
  89                 checkMatrix(n, n, z, ldz)
  90                 switch {
  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")
  95                 }
  96         }
  97
  98         // Quick return if possible.
  99         if n == 0 {
 100                 return 0
 101         }
 102         if ilo == ihi {
 103                 wr[ilo] = h[ilo*ldh+ilo]
 104                 wi[ilo] = 0
 105                 return 0
 106         }
 107
 108         // Clear out the trash.
 109         for j := ilo; j < ihi-2; j++ {
 110                 h[(j+2)*ldh+j] = 0
 111                 h[(j+3)*ldh+j] = 0
 112         }
 113         if ilo <= ihi-2 {
 114                 h[ihi*ldh+ihi-2] = 0
 115         }
 116
 117         nh := ihi - ilo + 1
 118         nz := ihiz - iloz + 1
 119
 120         // Set machine-dependent constants for the stopping criterion.
 121         ulp := dlamchP
 122         smlnum := float64(nh) / ulp * dlamchS
 123
 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.
 127         var i1, i2 int
 128         if wantt {
 129                 i1 = 0
 130                 i2 = n - 1
 131         }
 132
 133         itmax := 30 * max(10, nh) // Total number of QR iterations allowed.
 134
 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()
 141         i := ihi
 142         for i >= ilo {
 143                 l := ilo
 144
 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.
 148                 converged := false
 149                 for its := 0; its <= itmax; its++ {
 150                         // Look for a single small subdiagonal element.
 151                         var k int
 152                         for k = i; k > l; k-- {
 153                                 if math.Abs(h[k*ldh+k-1]) <= smlnum {
 154                                         break
 155                                 }
 156                                 tst := math.Abs(h[(k-1)*ldh+k-1]) + math.Abs(h[k*ldh+k])
 157                                 if tst == 0 {
 158                                         if k-2 >= ilo {
 159                                                 tst += math.Abs(h[(k-1)*ldh+k-2])
 160                                         }
 161                                         if k+1 <= ihi {
 162                                                 tst += math.Abs(h[(k+1)*ldh+k])
 163                                         }
 164                                 }
 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
 169                                 // in some cases.
 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]))
 175                                         s := aa + ab
 176                                         if ab/s*ba <= math.Max(smlnum, aa/s*bb*ulp) {
 177                                                 break
 178                                         }
 179                                 }
 180                         }
 181                         l = k
 182                         if l > ilo {
 183                                 // H[l,l-1] is negligible.
 184                                 h[l*ldh+l-1] = 0
 185                         }
 186                         if l >= i-1 {
 187                                 // Break the loop because a submatrix of order 1
 188                                 // or 2 has split off.
 189                                 converged = true
 190                                 break
 191                         }
 192
 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.
 196                         if !wantt {
 197                                 i1 = l
 198                                 i2 = i
 199                         }
 200
 201                         const (
 202                                 dat1 = 3.0
 203                                 dat2 = -0.4375
 204                         )
 205                         var h11, h21, h12, h22 float64
 206                         switch its {
 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]
 210                                 h12 = dat2 * s
 211                                 h21 = s
 212                                 h22 = h11
 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]
 216                                 h12 = dat2 * s
 217                                 h21 = s
 218                                 h22 = h11
 219                         default: // Prepare to use Francis' double shift (i.e.,
 220                                 // 2nd degree generalized Rayleigh quotient).
 221                                 h11 = h[(i-1)*ldh+i-1]
 222                                 h21 = h[i*ldh+i-1]
 223                                 h12 = h[(i-1)*ldh+i]
 224                                 h22 = h[i*ldh+i]
 225                         }
 226                         s := math.Abs(h11) + math.Abs(h12) + math.Abs(h21) + math.Abs(h22)
 227                         var (
 228                                 rt1r, rt1i float64
 229                                 rt2r, rt2i float64
 230                         )
 231                         if s != 0 {
 232                                 h11 /= s
 233                                 h21 /= s
 234                                 h12 /= s
 235                                 h22 /= s
 236                                 tr := (h11 + h22) / 2
 237                                 det := (h11-tr)*(h22-tr) - h12*h21
 238                                 rtdisc := math.Sqrt(math.Abs(det))
 239                                 if det >= 0 {
 240                                         // Complex conjugate shifts.
 241                                         rt1r = tr * s
 242                                         rt2r = rt1r
 243                                         rt1i = rtdisc * s
 244                                         rt2i = -rt1i
 245                                 } else {
 246                                         // Real shifts (use only one of them).
 247                                         rt1r = tr + rtdisc
 248                                         rt2r = tr - rtdisc
 249                                         if math.Abs(rt1r-h22) <= math.Abs(rt2r-h22) {
 250                                                 rt1r *= s
 251                                                 rt2r = rt1r
 252                                         } else {
 253                                                 rt2r *= s
 254                                                 rt1r = rt2r
 255                                         }
 256                                         rt1i = 0
 257                                         rt2i = 0
 258                                 }
 259                         }
 260
 261                         // Look for two consecutive small subdiagonal elements.
 262                         var m int
 263                         var v [3]float64
 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
 269                                 // most underflows.
 270                                 h21s := h[(m+1)*ldh+m]
 271                                 s := math.Abs(h[m*ldh+m]-rt2r) + math.Abs(rt2i) + math.Abs(h21s)
 272                                 h21s /= s
 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])
 277                                 v[0] /= s
 278                                 v[1] /= s
 279                                 v[2] /= s
 280                                 if m == l {
 281                                         break
 282                                 }
 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 {
 285                                         break
 286                                 }
 287                         }
 288
 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.
 295                                 //
 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.
 301
 302                                 nr := min(3, i-k+1)
 303                                 if k > m {
 304                                         bi.Dcopy(nr, h[k*ldh+k-1:], ldh, v[:], 1)
 305                                 }
 306                                 var t0 float64
 307                                 v[0], t0 = impl.Dlarfg(nr, v[0], v[1:], 1)
 308                                 if k > m {
 309                                         h[k*ldh+k-1] = v[0]
 310                                         h[(k+1)*ldh+k-1] = 0
 311                                         if k < i-1 {
 312                                                 h[(k+2)*ldh+k-1] = 0
 313                                         }
 314                                 } else if m > l {
 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
 318                                 }
 319                                 t1 := t0 * v[1]
 320                                 if nr == 3 {
 321                                         t2 := t0 * v[2]
 322
 323                                         // Apply G from the left to transform
 324                                         // the rows of the matrix in columns k
 325                                         // to i2.
 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
 331                                         }
 332
 333                                         // Apply G from the right to transform
 334                                         // the columns of the matrix in rows i1
 335                                         // to min(k+3,i).
 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
 341                                         }
 342
 343                                         if wantz {
 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
 350                                                 }
 351                                         }
 352                                 } else if nr == 2 {
 353                                         // Apply G from the left to transform
 354                                         // the rows of the matrix in columns k
 355                                         // to i2.
 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
 360                                         }
 361
 362                                         // Apply G from the right to transform
 363                                         // the columns of the matrix in rows i1
 364                                         // to min(k+3,i).
 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
 369                                         }
 370
 371                                         if wantz {
 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
 377                                                 }
 378                                         }
 379                                 }
 380                         }
 381                 }
 382
 383                 if !converged {
 384                         // The QR iteration finished without splitting off a
 385                         // submatrix of order 1 or 2.
 386                         return i + 1
 387                 }
 388
 389                 if l == i {
 390                         // H[i,i-1] is negligible: one eigenvalue has converged.
 391                         wr[i] = h[i*ldh+i]
 392                         wi[i] = 0
 393                 } else if l == i-1 {
 394                         // H[i-1,i-2] is negligible: a pair of eigenvalues have converged.
 395
 396                         // Transform the 2×2 submatrix to standard Schur form,
 397                         // and compute and store the eigenvalues.
 398                         var cs, sn float64
 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
 404
 405                         if wantt {
 406                                 // Apply the transformation to the rest of H.
 407                                 if i2 > i {
 408                                         bi.Drot(i2-i, h[(i-1)*ldh+i+1:], 1, h[i*ldh+i+1:], 1, cs, sn)
 409                                 }
 410                                 bi.Drot(i-i1-1, h[i1*ldh+i-1:], ldh, h[i1*ldh+i:], ldh, cs, sn)
 411                         }
 412
 413                         if wantz {
 414                                 // Apply the transformation to Z.
 415                                 bi.Drot(nz, z[iloz*ldz+i-1:], ldz, z[iloz*ldz+i:], ldz, cs, sn)
 416                         }
 417                 }
 418
 419                 // Return to start of the main loop with new value of i.
 420                 i = l - 1
 421         }
 422         return 0
 423 }