vendor/gonum.org/v1/gonum/lapack/gonum/dsteqr.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"
  11         "gonum.org/v1/gonum/blas/blas64"
  12         "gonum.org/v1/gonum/lapack"
  13 )
  14
  15 // Dsteqr computes the eigenvalues and optionally the eigenvectors of a symmetric
  16 // tridiagonal matrix using the implicit QL or QR method. The eigenvectors of a
  17 // full or band symmetric matrix can also be found if Dsytrd, Dsptrd, or Dsbtrd
  18 // have been used to reduce this matrix to tridiagonal form.
  19 //
  20 // d, on entry, contains the diagonal elements of the tridiagonal matrix. On exit,
  21 // d contains the eigenvalues in ascending order. d must have length n and
  22 // Dsteqr will panic otherwise.
  23 //
  24 // e, on entry, contains the off-diagonal elements of the tridiagonal matrix on
  25 // entry, and is overwritten during the call to Dsteqr. e must have length n-1 and
  26 // Dsteqr will panic otherwise.
  27 //
  28 // z, on entry, contains the n×n orthogonal matrix used in the reduction to
  29 // tridiagonal form if compz == lapack.OriginalEV. On exit, if
  30 // compz == lapack.OriginalEV, z contains the orthonormal eigenvectors of the
  31 // original symmetric matrix, and if compz == lapack.TridiagEV, z contains the
  32 // orthonormal eigenvectors of the symmetric tridiagonal matrix. z is not used
  33 // if compz == lapack.None.
  34 //
  35 // work must have length at least max(1, 2*n-2) if the eigenvectors are computed,
  36 // and Dsteqr will panic otherwise.
  37 //
  38 // Dsteqr is an internal routine. It is exported for testing purposes.
  39 func (impl Implementation) Dsteqr(compz lapack.EVComp, n int, d, e, z []float64, ldz int, work []float64) (ok bool) {
  40         if n < 0 {
  41                 panic(nLT0)
  42         }
  43         if len(d) < n {
  44                 panic(badD)
  45         }
  46         if len(e) < n-1 {
  47                 panic(badE)
  48         }
  49         if compz != lapack.None && compz != lapack.TridiagEV && compz != lapack.OriginalEV {
  50                 panic(badEVComp)
  51         }
  52         if compz != lapack.None {
  53                 if len(work) < max(1, 2*n-2) {
  54                         panic(badWork)
  55                 }
  56                 checkMatrix(n, n, z, ldz)
  57         }
  58
  59         var icompz int
  60         if compz == lapack.OriginalEV {
  61                 icompz = 1
  62         } else if compz == lapack.TridiagEV {
  63                 icompz = 2
  64         }
  65
  66         if n == 0 {
  67                 return true
  68         }
  69         if n == 1 {
  70                 if icompz == 2 {
  71                         z[0] = 1
  72                 }
  73                 return true
  74         }
  75
  76         bi := blas64.Implementation()
  77
  78         eps := dlamchE
  79         eps2 := eps * eps
  80         safmin := dlamchS
  81         safmax := 1 / safmin
  82         ssfmax := math.Sqrt(safmax) / 3
  83         ssfmin := math.Sqrt(safmin) / eps2
  84
  85         // Compute the eigenvalues and eigenvectors of the tridiagonal matrix.
  86         if icompz == 2 {
  87                 impl.Dlaset(blas.All, n, n, 0, 1, z, ldz)
  88         }
  89         const maxit = 30
  90         nmaxit := n * maxit
  91
  92         jtot := 0
  93
  94         // Determine where the matrix splits and choose QL or QR iteration for each
  95         // block, according to whether top or bottom diagonal element is smaller.
  96         l1 := 0
  97         nm1 := n - 1
  98
  99         type scaletype int
 100         const (
 101                 down scaletype = iota + 1
 102                 up
 103         )
 104         var iscale scaletype
 105
 106         for {
 107                 if l1 > n-1 {
 108                         // Order eigenvalues and eigenvectors.
 109                         if icompz == 0 {
 110                                 impl.Dlasrt(lapack.SortIncreasing, n, d)
 111                         } else {
 112                                 // TODO(btracey): Consider replacing this sort with a call to sort.Sort.
 113                                 for ii := 1; ii < n; ii++ {
 114                                         i := ii - 1
 115                                         k := i
 116                                         p := d[i]
 117                                         for j := ii; j < n; j++ {
 118                                                 if d[j] < p {
 119                                                         k = j
 120                                                         p = d[j]
 121                                                 }
 122                                         }
 123                                         if k != i {
 124                                                 d[k] = d[i]
 125                                                 d[i] = p
 126                                                 bi.Dswap(n, z[i:], ldz, z[k:], ldz)
 127                                         }
 128                                 }
 129                         }
 130                         return true
 131                 }
 132                 if l1 > 0 {
 133                         e[l1-1] = 0
 134                 }
 135                 var m int
 136                 if l1 <= nm1 {
 137                         for m = l1; m < nm1; m++ {
 138                                 test := math.Abs(e[m])
 139                                 if test == 0 {
 140                                         break
 141                                 }
 142                                 if test <= (math.Sqrt(math.Abs(d[m]))*math.Sqrt(math.Abs(d[m+1])))*eps {
 143                                         e[m] = 0
 144                                         break
 145                                 }
 146                         }
 147                 }
 148                 l := l1
 149                 lsv := l
 150                 lend := m
 151                 lendsv := lend
 152                 l1 = m + 1
 153                 if lend == l {
 154                         continue
 155                 }
 156
 157                 // Scale submatrix in rows and columns L to Lend
 158                 anorm := impl.Dlanst(lapack.MaxAbs, lend-l+1, d[l:], e[l:])
 159                 switch {
 160                 case anorm == 0:
 161                         continue
 162                 case anorm > ssfmax:
 163                         iscale = down
 164                         // Pretend that d and e are matrices with 1 column.
 165                         impl.Dlascl(lapack.General, 0, 0, anorm, ssfmax, lend-l+1, 1, d[l:], 1)
 166                         impl.Dlascl(lapack.General, 0, 0, anorm, ssfmax, lend-l, 1, e[l:], 1)
 167                 case anorm < ssfmin:
 168                         iscale = up
 169                         impl.Dlascl(lapack.General, 0, 0, anorm, ssfmin, lend-l+1, 1, d[l:], 1)
 170                         impl.Dlascl(lapack.General, 0, 0, anorm, ssfmin, lend-l, 1, e[l:], 1)
 171                 }
 172
 173                 // Choose between QL and QR.
 174                 if math.Abs(d[lend]) < math.Abs(d[l]) {
 175                         lend = lsv
 176                         l = lendsv
 177                 }
 178                 if lend > l {
 179                         // QL Iteration. Look for small subdiagonal element.
 180                         for {
 181                                 if l != lend {
 182                                         for m = l; m < lend; m++ {
 183                                                 v := math.Abs(e[m])
 184                                                 if v*v <= (eps2*math.Abs(d[m]))*math.Abs(d[m+1])+safmin {
 185                                                         break
 186                                                 }
 187                                         }
 188                                 } else {
 189                                         m = lend
 190                                 }
 191                                 if m < lend {
 192                                         e[m] = 0
 193                                 }
 194                                 p := d[l]
 195                                 if m == l {
 196                                         // Eigenvalue found.
 197                                         l++
 198                                         if l > lend {
 199                                                 break
 200                                         }
 201                                         continue
 202                                 }
 203
 204                                 // If remaining matrix is 2×2, use Dlae2 to compute its eigensystem.
 205                                 if m == l+1 {
 206                                         if icompz > 0 {
 207                                                 d[l], d[l+1], work[l], work[n-1+l] = impl.Dlaev2(d[l], e[l], d[l+1])
 208                                                 impl.Dlasr(blas.Right, lapack.Variable, lapack.Backward,
 209                                                         n, 2, work[l:], work[n-1+l:], z[l:], ldz)
 210                                         } else {
 211                                                 d[l], d[l+1] = impl.Dlae2(d[l], e[l], d[l+1])
 212                                         }
 213                                         e[l] = 0
 214                                         l += 2
 215                                         if l > lend {
 216                                                 break
 217                                         }
 218                                         continue
 219                                 }
 220
 221                                 if jtot == nmaxit {
 222                                         break
 223                                 }
 224                                 jtot++
 225
 226                                 // Form shift
 227                                 g := (d[l+1] - p) / (2 * e[l])
 228                                 r := impl.Dlapy2(g, 1)
 229                                 g = d[m] - p + e[l]/(g+math.Copysign(r, g))
 230                                 s := 1.0
 231                                 c := 1.0
 232                                 p = 0.0
 233
 234                                 // Inner loop
 235                                 for i := m - 1; i >= l; i-- {
 236                                         f := s * e[i]
 237                                         b := c * e[i]
 238                                         c, s, r = impl.Dlartg(g, f)
 239                                         if i != m-1 {
 240                                                 e[i+1] = r
 241                                         }
 242                                         g = d[i+1] - p
 243                                         r = (d[i]-g)*s + 2*c*b
 244                                         p = s * r
 245                                         d[i+1] = g + p
 246                                         g = c*r - b
 247
 248                                         // If eigenvectors are desired, then save rotations.
 249                                         if icompz > 0 {
 250                                                 work[i] = c
 251                                                 work[n-1+i] = -s
 252                                         }
 253                                 }
 254                                 // If eigenvectors are desired, then apply saved rotations.
 255                                 if icompz > 0 {
 256                                         mm := m - l + 1
 257                                         impl.Dlasr(blas.Right, lapack.Variable, lapack.Backward,
 258                                                 n, mm, work[l:], work[n-1+l:], z[l:], ldz)
 259                                 }
 260                                 d[l] -= p
 261                                 e[l] = g
 262                         }
 263                 } else {
 264                         // QR Iteration.
 265                         // Look for small superdiagonal element.
 266                         for {
 267                                 if l != lend {
 268                                         for m = l; m > lend; m-- {
 269                                                 v := math.Abs(e[m-1])
 270                                                 if v*v <= (eps2*math.Abs(d[m])*math.Abs(d[m-1]) + safmin) {
 271                                                         break
 272                                                 }
 273                                         }
 274                                 } else {
 275                                         m = lend
 276                                 }
 277                                 if m > lend {
 278                                         e[m-1] = 0
 279                                 }
 280                                 p := d[l]
 281                                 if m == l {
 282                                         // Eigenvalue found
 283                                         l--
 284                                         if l < lend {
 285                                                 break
 286                                         }
 287                                         continue
 288                                 }
 289
 290                                 // If remaining matrix is 2×2, use Dlae2 to compute its eigenvalues.
 291                                 if m == l-1 {
 292                                         if icompz > 0 {
 293                                                 d[l-1], d[l], work[m], work[n-1+m] = impl.Dlaev2(d[l-1], e[l-1], d[l])
 294                                                 impl.Dlasr(blas.Right, lapack.Variable, lapack.Forward,
 295                                                         n, 2, work[m:], work[n-1+m:], z[l-1:], ldz)
 296                                         } else {
 297                                                 d[l-1], d[l] = impl.Dlae2(d[l-1], e[l-1], d[l])
 298                                         }
 299                                         e[l-1] = 0
 300                                         l -= 2
 301                                         if l < lend {
 302                                                 break
 303                                         }
 304                                         continue
 305                                 }
 306                                 if jtot == nmaxit {
 307                                         break
 308                                 }
 309                                 jtot++
 310
 311                                 // Form shift.
 312                                 g := (d[l-1] - p) / (2 * e[l-1])
 313                                 r := impl.Dlapy2(g, 1)
 314                                 g = d[m] - p + (e[l-1])/(g+math.Copysign(r, g))
 315                                 s := 1.0
 316                                 c := 1.0
 317                                 p = 0.0
 318
 319                                 // Inner loop.
 320                                 for i := m; i < l; i++ {
 321                                         f := s * e[i]
 322                                         b := c * e[i]
 323                                         c, s, r = impl.Dlartg(g, f)
 324                                         if i != m {
 325                                                 e[i-1] = r
 326                                         }
 327                                         g = d[i] - p
 328                                         r = (d[i+1]-g)*s + 2*c*b
 329                                         p = s * r
 330                                         d[i] = g + p
 331                                         g = c*r - b
 332
 333                                         // If eigenvectors are desired, then save rotations.
 334                                         if icompz > 0 {
 335                                                 work[i] = c
 336                                                 work[n-1+i] = s
 337                                         }
 338                                 }
 339
 340                                 // If eigenvectors are desired, then apply saved rotations.
 341                                 if icompz > 0 {
 342                                         mm := l - m + 1
 343                                         impl.Dlasr(blas.Right, lapack.Variable, lapack.Forward,
 344                                                 n, mm, work[m:], work[n-1+m:], z[m:], ldz)
 345                                 }
 346                                 d[l] -= p
 347                                 e[l-1] = g
 348                         }
 349                 }
 350
 351                 // Undo scaling if necessary.
 352                 switch iscale {
 353                 case down:
 354                         // Pretend that d and e are matrices with 1 column.
 355                         impl.Dlascl(lapack.General, 0, 0, ssfmax, anorm, lendsv-lsv+1, 1, d[lsv:], 1)
 356                         impl.Dlascl(lapack.General, 0, 0, ssfmax, anorm, lendsv-lsv, 1, e[lsv:], 1)
 357                 case up:
 358                         impl.Dlascl(lapack.General, 0, 0, ssfmin, anorm, lendsv-lsv+1, 1, d[lsv:], 1)
 359                         impl.Dlascl(lapack.General, 0, 0, ssfmin, anorm, lendsv-lsv, 1, e[lsv:], 1)
 360                 }
 361
 362                 // Check for no convergence to an eigenvalue after a total of n*maxit iterations.
 363                 if jtot >= nmaxit {
 364                         break
 365                 }
 366         }
 367         for i := 0; i < n-1; i++ {
 368                 if e[i] != 0 {
 369                         return false
 370                 }
 371         }
 372         return true
 373 }