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

   1 // Copyright ©2015 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 // Dlasq2 computes all the eigenvalues of the symmetric positive
  14 // definite tridiagonal matrix associated with the qd array Z. Eigevalues
  15 // are computed to high relative accuracy avoiding denormalization, underflow
  16 // and overflow.
  17 //
  18 // To see the relation of Z to the tridiagonal matrix, let L be a
  19 // unit lower bidiagonal matrix with sub-diagonals Z(2,4,6,,..) and
  20 // let U be an upper bidiagonal matrix with 1's above and diagonal
  21 // Z(1,3,5,,..). The tridiagonal is L*U or, if you prefer, the
  22 // symmetric tridiagonal to which it is similar.
  23 //
  24 // info returns a status error. The return codes mean as follows:
  25 //  0: The algorithm completed successfully.
  26 //  1: A split was marked by a positive value in e.
  27 //  2: Current block of Z not diagonalized after 100*n iterations (in inner
  28 //     while loop). On exit Z holds a qd array with the same eigenvalues as
  29 //     the given Z.
  30 //  3: Termination criterion of outer while loop not met (program created more
  31 //     than N unreduced blocks).
  32 //
  33 // z must have length at least 4*n, and must not contain any negative elements.
  34 // Dlasq2 will panic otherwise.
  35 //
  36 // Dlasq2 is an internal routine. It is exported for testing purposes.
  37 func (impl Implementation) Dlasq2(n int, z []float64) (info int) {
  38         // TODO(btracey): make info an error.
  39         if len(z) < 4*n {
  40                 panic(badZ)
  41         }
  42         const cbias = 1.5
  43
  44         eps := dlamchP
  45         safmin := dlamchS
  46         tol := eps * 100
  47         tol2 := tol * tol
  48         if n < 0 {
  49                 panic(nLT0)
  50         }
  51         if n == 0 {
  52                 return info
  53         }
  54         if n == 1 {
  55                 if z[0] < 0 {
  56                         panic(negZ)
  57                 }
  58                 return info
  59         }
  60         if n == 2 {
  61                 if z[1] < 0 || z[2] < 0 {
  62                         panic("lapack: bad z value")
  63                 } else if z[2] > z[0] {
  64                         z[0], z[2] = z[2], z[0]
  65                 }
  66                 z[4] = z[0] + z[1] + z[2]
  67                 if z[1] > z[2]*tol2 {
  68                         t := 0.5 * (z[0] - z[2] + z[1])
  69                         s := z[2] * (z[1] / t)
  70                         if s <= t {
  71                                 s = z[2] * (z[1] / (t * (1 + math.Sqrt(1+s/t))))
  72                         } else {
  73                                 s = z[2] * (z[1] / (t + math.Sqrt(t)*math.Sqrt(t+s)))
  74                         }
  75                         t = z[0] + s + z[1]
  76                         z[2] *= z[0] / t
  77                         z[0] = t
  78                 }
  79                 z[1] = z[2]
  80                 z[5] = z[1] + z[0]
  81                 return info
  82         }
  83         // Check for negative data and compute sums of q's and e's.
  84         z[2*n-1] = 0
  85         emin := z[1]
  86         var d, e, qmax float64
  87         var i1, n1 int
  88         for k := 0; k < 2*(n-1); k += 2 {
  89                 if z[k] < 0 || z[k+1] < 0 {
  90                         panic("lapack: bad z value")
  91                 }
  92                 d += z[k]
  93                 e += z[k+1]
  94                 qmax = math.Max(qmax, z[k])
  95                 emin = math.Min(emin, z[k+1])
  96         }
  97         if z[2*(n-1)] < 0 {
  98                 panic("lapack: bad z value")
  99         }
 100         d += z[2*(n-1)]
 101         qmax = math.Max(qmax, z[2*(n-1)])
 102         // Check for diagonality.
 103         if e == 0 {
 104                 for k := 1; k < n; k++ {
 105                         z[k] = z[2*k]
 106                 }
 107                 impl.Dlasrt(lapack.SortDecreasing, n, z)
 108                 z[2*(n-1)] = d
 109                 return info
 110         }
 111         trace := d + e
 112         // Check for zero data.
 113         if trace == 0 {
 114                 z[2*(n-1)] = 0
 115                 return info
 116         }
 117         // Rearrange data for locality: Z=(q1,qq1,e1,ee1,q2,qq2,e2,ee2,...).
 118         for k := 2 * n; k >= 2; k -= 2 {
 119                 z[2*k-1] = 0
 120                 z[2*k-2] = z[k-1]
 121                 z[2*k-3] = 0
 122                 z[2*k-4] = z[k-2]
 123         }
 124         i0 := 0
 125         n0 := n - 1
 126
 127         // Reverse the qd-array, if warranted.
 128         // z[4*i0-3] --> z[4*(i0+1)-3-1] --> z[4*i0]
 129         if cbias*z[4*i0] < z[4*n0] {
 130                 ipn4Out := 4 * (i0 + n0 + 2)
 131                 for i4loop := 4 * (i0 + 1); i4loop <= 2*(i0+n0+1); i4loop += 4 {
 132                         i4 := i4loop - 1
 133                         ipn4 := ipn4Out - 1
 134                         z[i4-3], z[ipn4-i4-4] = z[ipn4-i4-4], z[i4-3]
 135                         z[i4-1], z[ipn4-i4-6] = z[ipn4-i4-6], z[i4-1]
 136                 }
 137         }
 138
 139         // Initial split checking via dqd and Li's test.
 140         pp := 0
 141         for k := 0; k < 2; k++ {
 142                 d = z[4*n0+pp]
 143                 for i4loop := 4*n0 + pp; i4loop >= 4*(i0+1)+pp; i4loop -= 4 {
 144                         i4 := i4loop - 1
 145                         if z[i4-1] <= tol2*d {
 146                                 z[i4-1] = math.Copysign(0, -1)
 147                                 d = z[i4-3]
 148                         } else {
 149                                 d = z[i4-3] * (d / (d + z[i4-1]))
 150                         }
 151                 }
 152                 // dqd maps Z to ZZ plus Li's test.
 153                 emin = z[4*(i0+1)+pp]
 154                 d = z[4*i0+pp]
 155                 for i4loop := 4*(i0+1) + pp; i4loop <= 4*n0+pp; i4loop += 4 {
 156                         i4 := i4loop - 1
 157                         z[i4-2*pp-2] = d + z[i4-1]
 158                         if z[i4-1] <= tol2*d {
 159                                 z[i4-1] = math.Copysign(0, -1)
 160                                 z[i4-2*pp-2] = d
 161                                 z[i4-2*pp] = 0
 162                                 d = z[i4+1]
 163                         } else if safmin*z[i4+1] < z[i4-2*pp-2] && safmin*z[i4-2*pp-2] < z[i4+1] {
 164                                 tmp := z[i4+1] / z[i4-2*pp-2]
 165                                 z[i4-2*pp] = z[i4-1] * tmp
 166                                 d *= tmp
 167                         } else {
 168                                 z[i4-2*pp] = z[i4+1] * (z[i4-1] / z[i4-2*pp-2])
 169                                 d = z[i4+1] * (d / z[i4-2*pp-2])
 170                         }
 171                         emin = math.Min(emin, z[i4-2*pp])
 172                 }
 173                 z[4*(n0+1)-pp-3] = d
 174
 175                 // Now find qmax.
 176                 qmax = z[4*(i0+1)-pp-3]
 177                 for i4loop := 4*(i0+1) - pp + 2; i4loop <= 4*(n0+1)+pp-2; i4loop += 4 {
 178                         i4 := i4loop - 1
 179                         qmax = math.Max(qmax, z[i4])
 180                 }
 181                 // Prepare for the next iteration on K.
 182                 pp = 1 - pp
 183         }
 184
 185         // Initialise variables to pass to DLASQ3.
 186         var ttype int
 187         var dmin1, dmin2, dn, dn1, dn2, g, tau float64
 188         var tempq float64
 189         iter := 2
 190         var nFail int
 191         nDiv := 2 * (n0 - i0)
 192         var i4 int
 193 outer:
 194         for iwhila := 1; iwhila <= n+1; iwhila++ {
 195                 // Test for completion.
 196                 if n0 < 0 {
 197                         // Move q's to the front.
 198                         for k := 1; k < n; k++ {
 199                                 z[k] = z[4*k]
 200                         }
 201                         // Sort and compute sum of eigenvalues.
 202                         impl.Dlasrt(lapack.SortDecreasing, n, z)
 203                         e = 0
 204                         for k := n - 1; k >= 0; k-- {
 205                                 e += z[k]
 206                         }
 207                         // Store trace, sum(eigenvalues) and information on performance.
 208                         z[2*n] = trace
 209                         z[2*n+1] = e
 210                         z[2*n+2] = float64(iter)
 211                         z[2*n+3] = float64(nDiv) / float64(n*n)
 212                         z[2*n+4] = 100 * float64(nFail) / float64(iter)
 213                         return info
 214                 }
 215
 216                 // While array unfinished do
 217                 // e[n0] holds the value of sigma when submatrix in i0:n0
 218                 // splits from the rest of the array, but is negated.
 219                 var desig float64
 220                 var sigma float64
 221                 if n0 != n-1 {
 222                         sigma = -z[4*(n0+1)-2]
 223                 }
 224                 if sigma < 0 {
 225                         info = 1
 226                         return info
 227                 }
 228                 // Find last unreduced submatrix's top index i0, find qmax and
 229                 // emin. Find Gershgorin-type bound if Q's much greater than E's.
 230                 var emax float64
 231                 if n0 > i0 {
 232                         emin = math.Abs(z[4*(n0+1)-6])
 233                 } else {
 234                         emin = 0
 235                 }
 236                 qmin := z[4*(n0+1)-4]
 237                 qmax = qmin
 238                 zSmall := false
 239                 for i4loop := 4 * (n0 + 1); i4loop >= 8; i4loop -= 4 {
 240                         i4 = i4loop - 1
 241                         if z[i4-5] <= 0 {
 242                                 zSmall = true
 243                                 break
 244                         }
 245                         if qmin >= 4*emax {
 246                                 qmin = math.Min(qmin, z[i4-3])
 247                                 emax = math.Max(emax, z[i4-5])
 248                         }
 249                         qmax = math.Max(qmax, z[i4-7]+z[i4-5])
 250                         emin = math.Min(emin, z[i4-5])
 251                 }
 252                 if !zSmall {
 253                         i4 = 3
 254                 }
 255                 i0 = (i4+1)/4 - 1
 256                 pp = 0
 257                 if n0-i0 > 1 {
 258                         dee := z[4*i0]
 259                         deemin := dee
 260                         kmin := i0
 261                         for i4loop := 4*(i0+1) + 1; i4loop <= 4*(n0+1)-3; i4loop += 4 {
 262                                 i4 := i4loop - 1
 263                                 dee = z[i4] * (dee / (dee + z[i4-2]))
 264                                 if dee <= deemin {
 265                                         deemin = dee
 266                                         kmin = (i4+4)/4 - 1
 267                                 }
 268                         }
 269                         if (kmin-i0)*2 < n0-kmin && deemin <= 0.5*z[4*n0] {
 270                                 ipn4Out := 4 * (i0 + n0 + 2)
 271                                 pp = 2
 272                                 for i4loop := 4 * (i0 + 1); i4loop <= 2*(i0+n0+1); i4loop += 4 {
 273                                         i4 := i4loop - 1
 274                                         ipn4 := ipn4Out - 1
 275                                         z[i4-3], z[ipn4-i4-4] = z[ipn4-i4-4], z[i4-3]
 276                                         z[i4-2], z[ipn4-i4-3] = z[ipn4-i4-3], z[i4-2]
 277                                         z[i4-1], z[ipn4-i4-6] = z[ipn4-i4-6], z[i4-1]
 278                                         z[i4], z[ipn4-i4-5] = z[ipn4-i4-5], z[i4]
 279                                 }
 280                         }
 281                 }
 282                 // Put -(initial shift) into DMIN.
 283                 dmin := -math.Max(0, qmin-2*math.Sqrt(qmin)*math.Sqrt(emax))
 284
 285                 // Now i0:n0 is unreduced.
 286                 // PP = 0 for ping, PP = 1 for pong.
 287                 // PP = 2 indicates that flipping was applied to the Z array and
 288                 //              and that the tests for deflation upon entry in Dlasq3
 289                 //              should not be performed.
 290                 nbig := 100 * (n0 - i0 + 1)
 291                 for iwhilb := 0; iwhilb < nbig; iwhilb++ {
 292                         if i0 > n0 {
 293                                 continue outer
 294                         }
 295
 296                         // While submatrix unfinished take a good dqds step.
 297                         i0, n0, pp, dmin, sigma, desig, qmax, nFail, iter, nDiv, ttype, dmin1, dmin2, dn, dn1, dn2, g, tau =
 298                                 impl.Dlasq3(i0, n0, z, pp, dmin, sigma, desig, qmax, nFail, iter, nDiv, ttype, dmin1, dmin2, dn, dn1, dn2, g, tau)
 299
 300                         pp = 1 - pp
 301                         // When emin is very small check for splits.
 302                         if pp == 0 && n0-i0 >= 3 {
 303                                 if z[4*(n0+1)-1] <= tol2*qmax || z[4*(n0+1)-2] <= tol2*sigma {
 304                                         splt := i0 - 1
 305                                         qmax = z[4*i0]
 306                                         emin = z[4*(i0+1)-2]
 307                                         oldemn := z[4*(i0+1)-1]
 308                                         for i4loop := 4 * (i0 + 1); i4loop <= 4*(n0-2); i4loop += 4 {
 309                                                 i4 := i4loop - 1
 310                                                 if z[i4] <= tol2*z[i4-3] || z[i4-1] <= tol2*sigma {
 311                                                         z[i4-1] = -sigma
 312                                                         splt = i4 / 4
 313                                                         qmax = 0
 314                                                         emin = z[i4+3]
 315                                                         oldemn = z[i4+4]
 316                                                 } else {
 317                                                         qmax = math.Max(qmax, z[i4+1])
 318                                                         emin = math.Min(emin, z[i4-1])
 319                                                         oldemn = math.Min(oldemn, z[i4])
 320                                                 }
 321                                         }
 322                                         z[4*(n0+1)-2] = emin
 323                                         z[4*(n0+1)-1] = oldemn
 324                                         i0 = splt + 1
 325                                 }
 326                         }
 327                 }
 328                 // Maximum number of iterations exceeded, restore the shift
 329                 // sigma and place the new d's and e's in a qd array.
 330                 // This might need to be done for several blocks.
 331                 info = 2
 332                 i1 = i0
 333                 n1 = n0
 334                 for {
 335                         tempq = z[4*i0]
 336                         z[4*i0] += sigma
 337                         for k := i0 + 1; k <= n0; k++ {
 338                                 tempe := z[4*(k+1)-6]
 339                                 z[4*(k+1)-6] *= tempq / z[4*(k+1)-8]
 340                                 tempq = z[4*k]
 341                                 z[4*k] += sigma + tempe - z[4*(k+1)-6]
 342                         }
 343                         // Prepare to do this on the previous block if there is one.
 344                         if i1 <= 0 {
 345                                 break
 346                         }
 347                         n1 = i1 - 1
 348                         for i1 >= 1 && z[4*(i1+1)-6] >= 0 {
 349                                 i1 -= 1
 350                         }
 351                         sigma = -z[4*(n1+1)-2]
 352                 }
 353                 for k := 0; k < n; k++ {
 354                         z[2*k] = z[4*k]
 355                         // Only the block 1..N0 is unfinished.  The rest of the e's
 356                         // must be essentially zero, although sometimes other data
 357                         // has been stored in them.
 358                         if k < n0 {
 359                                 z[2*(k+1)-1] = z[4*(k+1)-1]
 360                         } else {
 361                                 z[2*(k+1)] = 0
 362                         }
 363                 }
 364                 return info
 365         }
 366         info = 3
 367         return info
 368 }