github.com/gonum/lapack@v0.0.0-20181123203213-e4cdc5a0bff9/native/dlartg.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  // Dlartg generates a plane rotation so that
    10  //  [ cs sn] * [f] = [r]
    11  //  [-sn cs]   [g] = [0]
    12  // This is a more accurate version of BLAS drotg, with the other differences that
    13  // if g = 0, then cs = 1 and sn = 0, and if f = 0 and g != 0, then cs = 0 and sn = 1.
    14  // If abs(f) > abs(g), cs will be positive.
    15  //
    16  // Dlartg is an internal routine. It is exported for testing purposes.
    17  func (impl Implementation) Dlartg(f, g float64) (cs, sn, r float64) {
    18  	safmn2 := math.Pow(dlamchB, math.Trunc(math.Log(dlamchS/dlamchE)/math.Log(dlamchB)/2))
    19  	safmx2 := 1 / safmn2
    20  	if g == 0 {
    21  		cs = 1
    22  		sn = 0
    23  		r = f
    24  		return cs, sn, r
    25  	}
    26  	if f == 0 {
    27  		cs = 0
    28  		sn = 1
    29  		r = g
    30  		return cs, sn, r
    31  	}
    32  	f1 := f
    33  	g1 := g
    34  	scale := math.Max(math.Abs(f1), math.Abs(g1))
    35  	if scale >= safmx2 {
    36  		var count int
    37  		for {
    38  			count++
    39  			f1 *= safmn2
    40  			g1 *= safmn2
    41  			scale = math.Max(math.Abs(f1), math.Abs(g1))
    42  			if scale < safmx2 {
    43  				break
    44  			}
    45  		}
    46  		r = math.Sqrt(f1*f1 + g1*g1)
    47  		cs = f1 / r
    48  		sn = g1 / r
    49  		for i := 0; i < count; i++ {
    50  			r *= safmx2
    51  		}
    52  	} else if scale <= safmn2 {
    53  		var count int
    54  		for {
    55  			count++
    56  			f1 *= safmx2
    57  			g1 *= safmx2
    58  			scale = math.Max(math.Abs(f1), math.Abs(g1))
    59  			if scale >= safmn2 {
    60  				break
    61  			}
    62  		}
    63  		r = math.Sqrt(f1*f1 + g1*g1)
    64  		cs = f1 / r
    65  		sn = g1 / r
    66  		for i := 0; i < count; i++ {
    67  			r *= safmn2
    68  		}
    69  	} else {
    70  		r = math.Sqrt(f1*f1 + g1*g1)
    71  		cs = f1 / r
    72  		sn = g1 / r
    73  	}
    74  	if math.Abs(f) > math.Abs(g) && cs < 0 {
    75  		cs *= -1
    76  		sn *= -1
    77  		r *= -1
    78  	}
    79  	return cs, sn, r
    80  }