vendor/gonum.org/v1/gonum/lapack/gonum/dsterf.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/lapack"
  11 )
  12
  13 // Dsterf computes all eigenvalues of a symmetric tridiagonal matrix using the
  14 // Pal-Walker-Kahan variant of the QL or QR algorithm.
  15 //
  16 // d contains the diagonal elements of the tridiagonal matrix on entry, and
  17 // contains the eigenvalues in ascending order on exit. d must have length at
  18 // least n, or Dsterf will panic.
  19 //
  20 // e contains the off-diagonal elements of the tridiagonal matrix on entry, and is
  21 // overwritten during the call to Dsterf. e must have length of at least n-1 or
  22 // Dsterf will panic.
  23 //
  24 // Dsterf is an internal routine. It is exported for testing purposes.
  25 func (impl Implementation) Dsterf(n int, d, e []float64) (ok bool) {
  26         if n < 0 {
  27                 panic(nLT0)
  28         }
  29         if n == 0 {
  30                 return true
  31         }
  32         if len(d) < n {
  33                 panic(badD)
  34         }
  35         if len(e) < n-1 {
  36                 panic(badE)
  37         }
  38
  39         const (
  40                 none = 0 // The values are not scaled.
  41                 down = 1 // The values are scaled below ssfmax threshold.
  42                 up   = 2 // The values are scaled below ssfmin threshold.
  43         )
  44
  45         // Determine the unit roundoff for this environment.
  46         eps := dlamchE
  47         eps2 := eps * eps
  48         safmin := dlamchS
  49         safmax := 1 / safmin
  50         ssfmax := math.Sqrt(safmax) / 3
  51         ssfmin := math.Sqrt(safmin) / eps2
  52
  53         // Compute the eigenvalues of the tridiagonal matrix.
  54         maxit := 30
  55         nmaxit := n * maxit
  56         jtot := 0
  57
  58         l1 := 0
  59
  60         for {
  61                 if l1 > n-1 {
  62                         impl.Dlasrt(lapack.SortIncreasing, n, d)
  63                         return true
  64                 }
  65                 if l1 > 0 {
  66                         e[l1-1] = 0
  67                 }
  68                 var m int
  69                 for m = l1; m < n-1; m++ {
  70                         if math.Abs(e[m]) <= math.Sqrt(math.Abs(d[m]))*math.Sqrt(math.Abs(d[m+1]))*eps {
  71                                 e[m] = 0
  72                                 break
  73                         }
  74                 }
  75
  76                 l := l1
  77                 lsv := l
  78                 lend := m
  79                 lendsv := lend
  80                 l1 = m + 1
  81                 if lend == 0 {
  82                         continue
  83                 }
  84
  85                 // Scale submatrix in rows and columns l to lend.
  86                 anorm := impl.Dlanst(lapack.MaxAbs, lend-l+1, d[l:], e[l:])
  87                 iscale := none
  88                 if anorm == 0 {
  89                         continue
  90                 }
  91                 if anorm > ssfmax {
  92                         iscale = down
  93                         impl.Dlascl(lapack.General, 0, 0, anorm, ssfmax, lend-l+1, 1, d[l:], n)
  94                         impl.Dlascl(lapack.General, 0, 0, anorm, ssfmax, lend-l, 1, e[l:], n)
  95                 } else if anorm < ssfmin {
  96                         iscale = up
  97                         impl.Dlascl(lapack.General, 0, 0, anorm, ssfmin, lend-l+1, 1, d[l:], n)
  98                         impl.Dlascl(lapack.General, 0, 0, anorm, ssfmin, lend-l, 1, e[l:], n)
  99                 }
 100
 101                 el := e[l:lend]
 102                 for i, v := range el {
 103                         el[i] *= v
 104                 }
 105
 106                 // Choose between QL and QR iteration.
 107                 if math.Abs(d[lend]) < math.Abs(d[l]) {
 108                         lend = lsv
 109                         l = lendsv
 110                 }
 111                 if lend >= l {
 112                         // QL Iteration.
 113                         // Look for small sub-diagonal element.
 114                         for {
 115                                 if l != lend {
 116                                         for m = l; m < lend; m++ {
 117                                                 if math.Abs(e[m]) <= eps2*(math.Abs(d[m]*d[m+1])) {
 118                                                         break
 119                                                 }
 120                                         }
 121                                 } else {
 122                                         m = lend
 123                                 }
 124                                 if m < lend {
 125                                         e[m] = 0
 126                                 }
 127                                 p := d[l]
 128                                 if m == l {
 129                                         // Eigenvalue found.
 130                                         l++
 131                                         if l > lend {
 132                                                 break
 133                                         }
 134                                         continue
 135                                 }
 136                                 // If remaining matrix is 2 by 2, use Dlae2 to compute its eigenvalues.
 137                                 if m == l+1 {
 138                                         d[l], d[l+1] = impl.Dlae2(d[l], math.Sqrt(e[l]), d[l+1])
 139                                         e[l] = 0
 140                                         l += 2
 141                                         if l > lend {
 142                                                 break
 143                                         }
 144                                         continue
 145                                 }
 146                                 if jtot == nmaxit {
 147                                         break
 148                                 }
 149                                 jtot++
 150
 151                                 // Form shift.
 152                                 rte := math.Sqrt(e[l])
 153                                 sigma := (d[l+1] - p) / (2 * rte)
 154                                 r := impl.Dlapy2(sigma, 1)
 155                                 sigma = p - (rte / (sigma + math.Copysign(r, sigma)))
 156
 157                                 c := 1.0
 158                                 s := 0.0
 159                                 gamma := d[m] - sigma
 160                                 p = gamma * gamma
 161
 162                                 // Inner loop.
 163                                 for i := m - 1; i >= l; i-- {
 164                                         bb := e[i]
 165                                         r := p + bb
 166                                         if i != m-1 {
 167                                                 e[i+1] = s * r
 168                                         }
 169                                         oldc := c
 170                                         c = p / r
 171                                         s = bb / r
 172                                         oldgam := gamma
 173                                         alpha := d[i]
 174                                         gamma = c*(alpha-sigma) - s*oldgam
 175                                         d[i+1] = oldgam + (alpha - gamma)
 176                                         if c != 0 {
 177                                                 p = (gamma * gamma) / c
 178                                         } else {
 179                                                 p = oldc * bb
 180                                         }
 181                                 }
 182                                 e[l] = s * p
 183                                 d[l] = sigma + gamma
 184                         }
 185                 } else {
 186                         for {
 187                                 // QR Iteration.
 188                                 // Look for small super-diagonal element.
 189                                 for m = l; m > lend; m-- {
 190                                         if math.Abs(e[m-1]) <= eps2*math.Abs(d[m]*d[m-1]) {
 191                                                 break
 192                                         }
 193                                 }
 194                                 if m > lend {
 195                                         e[m-1] = 0
 196                                 }
 197                                 p := d[l]
 198                                 if m == l {
 199                                         // Eigenvalue found.
 200                                         l--
 201                                         if l < lend {
 202                                                 break
 203                                         }
 204                                         continue
 205                                 }
 206
 207                                 // If remaining matrix is 2 by 2, use Dlae2 to compute its eigenvalues.
 208                                 if m == l-1 {
 209                                         d[l], d[l-1] = impl.Dlae2(d[l], math.Sqrt(e[l-1]), d[l-1])
 210                                         e[l-1] = 0
 211                                         l -= 2
 212                                         if l < lend {
 213                                                 break
 214                                         }
 215                                         continue
 216                                 }
 217                                 if jtot == nmaxit {
 218                                         break
 219                                 }
 220                                 jtot++
 221
 222                                 // Form shift.
 223                                 rte := math.Sqrt(e[l-1])
 224                                 sigma := (d[l-1] - p) / (2 * rte)
 225                                 r := impl.Dlapy2(sigma, 1)
 226                                 sigma = p - (rte / (sigma + math.Copysign(r, sigma)))
 227
 228                                 c := 1.0
 229                                 s := 0.0
 230                                 gamma := d[m] - sigma
 231                                 p = gamma * gamma
 232
 233                                 // Inner loop.
 234                                 for i := m; i < l; i++ {
 235                                         bb := e[i]
 236                                         r := p + bb
 237                                         if i != m {
 238                                                 e[i-1] = s * r
 239                                         }
 240                                         oldc := c
 241                                         c = p / r
 242                                         s = bb / r
 243                                         oldgam := gamma
 244                                         alpha := d[i+1]
 245                                         gamma = c*(alpha-sigma) - s*oldgam
 246                                         d[i] = oldgam + alpha - gamma
 247                                         if c != 0 {
 248                                                 p = (gamma * gamma) / c
 249                                         } else {
 250                                                 p = oldc * bb
 251                                         }
 252                                 }
 253                                 e[l-1] = s * p
 254                                 d[l] = sigma + gamma
 255                         }
 256                 }
 257
 258                 // Undo scaling if necessary
 259                 switch iscale {
 260                 case down:
 261                         impl.Dlascl(lapack.General, 0, 0, ssfmax, anorm, lendsv-lsv+1, 1, d[lsv:], n)
 262                 case up:
 263                         impl.Dlascl(lapack.General, 0, 0, ssfmin, anorm, lendsv-lsv+1, 1, d[lsv:], n)
 264                 }
 265
 266                 // Check for no convergence to an eigenvalue after a total of n*maxit iterations.
 267                 if jtot >= nmaxit {
 268                         break
 269                 }
 270         }
 271         for _, v := range e[:n-1] {
 272                 if v != 0 {
 273                         return false
 274                 }
 275         }
 276         impl.Dlasrt(lapack.SortIncreasing, n, d)
 277         return true
 278 }