vendor/gonum.org/v1/gonum/lapack/gonum/dlaqr1.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 "math"
   8
   9 // Dlaqr1 sets v to a scalar multiple of the first column of the product
  10 //  (H - (sr1 + i*si1)*I)*(H - (sr2 + i*si2)*I)
  11 // where H is a 2×2 or 3×3 matrix, I is the identity matrix of the same size,
  12 // and i is the imaginary unit. Scaling is done to avoid overflows and most
  13 // underflows.
  14 //
  15 // n is the order of H and must be either 2 or 3. It must hold that either sr1 =
  16 // sr2 and si1 = -si2, or si1 = si2 = 0. The length of v must be equal to n. If
  17 // any of these conditions is not met, Dlaqr1 will panic.
  18 //
  19 // Dlaqr1 is an internal routine. It is exported for testing purposes.
  20 func (impl Implementation) Dlaqr1(n int, h []float64, ldh int, sr1, si1, sr2, si2 float64, v []float64) {
  21         if n != 2 && n != 3 {
  22                 panic(badDims)
  23         }
  24         checkMatrix(n, n, h, ldh)
  25         if len(v) != n {
  26                 panic(badSlice)
  27         }
  28         if !((sr1 == sr2 && si1 == -si2) || (si1 == 0 && si2 == 0)) {
  29                 panic(badShifts)
  30         }
  31
  32         if n == 2 {
  33                 s := math.Abs(h[0]-sr2) + math.Abs(si2) + math.Abs(h[ldh])
  34                 if s == 0 {
  35                         v[0] = 0
  36                         v[1] = 0
  37                 } else {
  38                         h21s := h[ldh] / s
  39                         v[0] = h21s*h[1] + (h[0]-sr1)*((h[0]-sr2)/s) - si1*(si2/s)
  40                         v[1] = h21s * (h[0] + h[ldh+1] - sr1 - sr2)
  41                 }
  42                 return
  43         }
  44
  45         s := math.Abs(h[0]-sr2) + math.Abs(si2) + math.Abs(h[ldh]) + math.Abs(h[2*ldh])
  46         if s == 0 {
  47                 v[0] = 0
  48                 v[1] = 0
  49                 v[2] = 0
  50         } else {
  51                 h21s := h[ldh] / s
  52                 h31s := h[2*ldh] / s
  53                 v[0] = (h[0]-sr1)*((h[0]-sr2)/s) - si1*(si2/s) + h[1]*h21s + h[2]*h31s
  54                 v[1] = h21s*(h[0]+h[ldh+1]-sr1-sr2) + h[ldh+2]*h31s
  55                 v[2] = h31s*(h[0]+h[2*ldh+2]-sr1-sr2) + h21s*h[2*ldh+1]
  56         }
  57 }