vendor/gonum.org/v1/gonum/lapack/gonum/dlaev2.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 "math"
   8
   9 // Dlaev2 computes the Eigen decomposition of a symmetric 2×2 matrix.
  10 // The matrix is given by
  11 //  [a b]
  12 //  [b c]
  13 // Dlaev2 returns rt1 and rt2, the eigenvalues of the matrix where |RT1| > |RT2|,
  14 // and [cs1, sn1] which is the unit right eigenvalue for RT1.
  15 //  [ cs1 sn1] [a b] [cs1 -sn1] = [rt1   0]
  16 //  [-sn1 cs1] [b c] [sn1  cs1]   [  0 rt2]
  17 //
  18 // Dlaev2 is an internal routine. It is exported for testing purposes.
  19 func (impl Implementation) Dlaev2(a, b, c float64) (rt1, rt2, cs1, sn1 float64) {
  20         sm := a + c
  21         df := a - c
  22         adf := math.Abs(df)
  23         tb := b + b
  24         ab := math.Abs(tb)
  25         acmx := c
  26         acmn := a
  27         if math.Abs(a) > math.Abs(c) {
  28                 acmx = a
  29                 acmn = c
  30         }
  31         var rt float64
  32         if adf > ab {
  33                 rt = adf * math.Sqrt(1+(ab/adf)*(ab/adf))
  34         } else if adf < ab {
  35                 rt = ab * math.Sqrt(1+(adf/ab)*(adf/ab))
  36         } else {
  37                 rt = ab * math.Sqrt(2)
  38         }
  39         var sgn1 float64
  40         if sm < 0 {
  41                 rt1 = 0.5 * (sm - rt)
  42                 sgn1 = -1
  43                 rt2 = (acmx/rt1)*acmn - (b/rt1)*b
  44         } else if sm > 0 {
  45                 rt1 = 0.5 * (sm + rt)
  46                 sgn1 = 1
  47                 rt2 = (acmx/rt1)*acmn - (b/rt1)*b
  48         } else {
  49                 rt1 = 0.5 * rt
  50                 rt2 = -0.5 * rt
  51                 sgn1 = 1
  52         }
  53         var cs, sgn2 float64
  54         if df >= 0 {
  55                 cs = df + rt
  56                 sgn2 = 1
  57         } else {
  58                 cs = df - rt
  59                 sgn2 = -1
  60         }
  61         acs := math.Abs(cs)
  62         if acs > ab {
  63                 ct := -tb / cs
  64                 sn1 = 1 / math.Sqrt(1+ct*ct)
  65                 cs1 = ct * sn1
  66         } else {
  67                 if ab == 0 {
  68                         cs1 = 1
  69                         sn1 = 0
  70                 } else {
  71                         tn := -cs / tb
  72                         cs1 = 1 / math.Sqrt(1+tn*tn)
  73                         sn1 = tn * cs1
  74                 }
  75         }
  76         if sgn1 == sgn2 {
  77                 tn := cs1
  78                 cs1 = -sn1
  79                 sn1 = tn
  80         }
  81         return rt1, rt2, cs1, sn1
  82 }