vendor/gonum.org/v1/gonum/lapack/gonum/dlanv2.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 // Dlanv2 computes the Schur factorization of a real 2×2 matrix:
  10 //  [ a b ] = [ cs -sn ] * [ aa bb ] * [ cs sn ]
  11 //  [ c d ]   [ sn  cs ]   [ cc dd ] * [-sn cs ]
  12 // If cc is zero, aa and dd are real eigenvalues of the matrix. Otherwise it
  13 // holds that aa = dd and bb*cc < 0, and aa ± sqrt(bb*cc) are complex conjugate
  14 // eigenvalues. The real and imaginary parts of the eigenvalues are returned in
  15 // (rt1r,rt1i) and (rt2r,rt2i).
  16 func (impl Implementation) Dlanv2(a, b, c, d float64) (aa, bb, cc, dd float64, rt1r, rt1i, rt2r, rt2i float64, cs, sn float64) {
  17         switch {
  18         case c == 0: // Matrix is already upper triangular.
  19                 aa = a
  20                 bb = b
  21                 cc = 0
  22                 dd = d
  23                 cs = 1
  24                 sn = 0
  25         case b == 0: // Matrix is lower triangular, swap rows and columns.
  26                 aa = d
  27                 bb = -c
  28                 cc = 0
  29                 dd = a
  30                 cs = 0
  31                 sn = 1
  32         case a == d && math.Signbit(b) != math.Signbit(c): // Matrix is already in the standard Schur form.
  33                 aa = a
  34                 bb = b
  35                 cc = c
  36                 dd = d
  37                 cs = 1
  38                 sn = 0
  39         default:
  40                 temp := a - d
  41                 p := temp / 2
  42                 bcmax := math.Max(math.Abs(b), math.Abs(c))
  43                 bcmis := math.Min(math.Abs(b), math.Abs(c))
  44                 if b*c < 0 {
  45                         bcmis *= -1
  46                 }
  47                 scale := math.Max(math.Abs(p), bcmax)
  48                 z := p/scale*p + bcmax/scale*bcmis
  49                 eps := dlamchP
  50
  51                 if z >= 4*eps {
  52                         // Real eigenvalues. Compute aa and dd.
  53                         if p > 0 {
  54                                 z = p + math.Sqrt(scale)*math.Sqrt(z)
  55                         } else {
  56                                 z = p - math.Sqrt(scale)*math.Sqrt(z)
  57                         }
  58                         aa = d + z
  59                         dd = d - bcmax/z*bcmis
  60                         // Compute bb and the rotation matrix.
  61                         tau := impl.Dlapy2(c, z)
  62                         cs = z / tau
  63                         sn = c / tau
  64                         bb = b - c
  65                         cc = 0
  66                 } else {
  67                         // Complex eigenvalues, or real (almost) equal eigenvalues.
  68                         // Make diagonal elements equal.
  69                         sigma := b + c
  70                         tau := impl.Dlapy2(sigma, temp)
  71                         cs = math.Sqrt((1 + math.Abs(sigma)/tau) / 2)
  72                         sn = -p / (tau * cs)
  73                         if sigma < 0 {
  74                                 sn *= -1
  75                         }
  76                         // Compute [ aa bb ] = [ a b ] [ cs -sn ]
  77                         //         [ cc dd ]   [ c d ] [ sn  cs ]
  78                         aa = a*cs + b*sn
  79                         bb = -a*sn + b*cs
  80                         cc = c*cs + d*sn
  81                         dd = -c*sn + d*cs
  82                         // Compute [ a b ] = [ cs sn ] [ aa bb ]
  83                         //         [ c d ]   [-sn cs ] [ cc dd ]
  84                         a = aa*cs + cc*sn
  85                         b = bb*cs + dd*sn
  86                         c = -aa*sn + cc*cs
  87                         d = -bb*sn + dd*cs
  88
  89                         temp = (a + d) / 2
  90                         aa = temp
  91                         bb = b
  92                         cc = c
  93                         dd = temp
  94
  95                         if cc != 0 {
  96                                 if bb != 0 {
  97                                         if math.Signbit(bb) == math.Signbit(cc) {
  98                                                 // Real eigenvalues, reduce to
  99                                                 // upper triangular form.
 100                                                 sab := math.Sqrt(math.Abs(bb))
 101                                                 sac := math.Sqrt(math.Abs(cc))
 102                                                 p = sab * sac
 103                                                 if cc < 0 {
 104                                                         p *= -1
 105                                                 }
 106                                                 tau = 1 / math.Sqrt(math.Abs(bb+cc))
 107                                                 aa = temp + p
 108                                                 bb = bb - cc
 109                                                 cc = 0
 110                                                 dd = temp - p
 111                                                 cs1 := sab * tau
 112                                                 sn1 := sac * tau
 113                                                 cs, sn = cs*cs1-sn*sn1, cs*sn1+sn+cs1
 114                                         }
 115                                 } else {
 116                                         bb = -cc
 117                                         cc = 0
 118                                         cs, sn = -sn, cs
 119                                 }
 120                         }
 121                 }
 122         }
 123
 124         // Store eigenvalues in (rt1r,rt1i) and (rt2r,rt2i).
 125         rt1r = aa
 126         rt2r = dd
 127         if cc != 0 {
 128                 rt1i = math.Sqrt(math.Abs(bb)) * math.Sqrt(math.Abs(cc))
 129                 rt2i = -rt1i
 130         }
 131         return
 132 }