github.com/jingcheng-WU/gonum@v0.9.1-0.20210323123734-f1a2a11a8f7b/lapack/gonum/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 gonum
     6  
     7  import "math"
     8  
     9  // Dlartg generates a plane rotation so that
    10  //  [ cs sn] * [f] = [r]
    11  //  [-sn cs]   [g] = [0]
    12  // where cs*cs + sn*sn = 1.
    13  //
    14  // This is a more accurate version of BLAS Drotg, with the other differences
    15  // that
    16  //  - if g = 0, then cs = 1 and sn = 0
    17  //  - if f = 0 and g != 0, then cs = 0 and sn = 1
    18  //  - r takes the sign of f and so cs is always non-negative
    19  //
    20  // Dlartg is an internal routine. It is exported for testing purposes.
    21  func (impl Implementation) Dlartg(f, g float64) (cs, sn, r float64) {
    22  	// Implementation based on Supplemental Material to:
    23  	// Edward Anderson. 2017. Algorithm 978: Safe Scaling in the Level 1 BLAS.
    24  	// ACM Trans. Math. Softw. 44, 1, Article 12 (July 2017), 28 pages.
    25  	// DOI: https://doi.org/10.1145/3061665
    26  	const safmin = dlamchS
    27  	const safmax = 1 / safmin
    28  	f1 := math.Abs(f)
    29  	g1 := math.Abs(g)
    30  	switch {
    31  	case g == 0:
    32  		cs = 1
    33  		sn = 0
    34  		r = f
    35  	case f == 0:
    36  		cs = 0
    37  		sn = math.Copysign(1, g)
    38  		r = g1
    39  	case drtmin < f1 && f1 < drtmax && drtmin < g1 && g1 < drtmax:
    40  		d := math.Sqrt(f*f + g*g)
    41  		p := 1 / d
    42  		cs = f1 * p
    43  		sn = g * math.Copysign(p, f)
    44  		r = math.Copysign(d, f)
    45  	default:
    46  		maxfg := math.Max(f1, g1)
    47  		u := math.Min(math.Max(safmin, maxfg), safmax)
    48  		uu := 1 / u
    49  		fs := f * uu
    50  		gs := g * uu
    51  		d := math.Sqrt(fs*fs + gs*gs)
    52  		p := 1 / d
    53  		cs = math.Abs(fs) * p
    54  		sn = gs * math.Copysign(p, f)
    55  		r = math.Copysign(d, f) * u
    56  	}
    57  	return cs, sn, r
    58  }