gonum.org/v1/gonum@v0.14.0/lapack/gonum/dlaev2.go (about)

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