github.com/gonum/lapack@v0.0.0-20181123203213-e4cdc5a0bff9/native/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 native
     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  }