gonum.org/v1/gonum@v0.14.0/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  //
    11  //	[ cs sn] * [f] = [r]
    12  //	[-sn cs]   [g] = [0]
    13  //
    14  // where cs*cs + sn*sn = 1.
    15  //
    16  // This is a more accurate version of BLAS Drotg that uses scaling to avoid
    17  // overflow or underflow, with the other differences that
    18  //   - cs >= 0
    19  //   - if g = 0, then cs = 1 and sn = 0
    20  //   - if f = 0 and g != 0, then cs = 0 and sn = sign(1,g)
    21  //
    22  // Dlartg is an internal routine. It is exported for testing purposes.
    23  func (impl Implementation) Dlartg(f, g float64) (cs, sn, r float64) {
    24  	// Implementation based on Supplemental Material to:
    25  	//
    26  	// Edward Anderson
    27  	// Algorithm 978: Safe Scaling in the Level 1 BLAS
    28  	// ACM Trans. Math. Softw. 44, 1, Article 12 (2017)
    29  	// DOI: https://doi.org/10.1145/3061665
    30  	//
    31  	// For further details see:
    32  	//
    33  	// W. Pereira, A. Lotfi, J. Langou
    34  	// Numerical analysis of Givens rotation
    35  	// DOI: https://doi.org/10.48550/arXiv.2211.04010
    36  
    37  	if g == 0 {
    38  		return 1, 0, f
    39  	}
    40  
    41  	g1 := math.Abs(g)
    42  
    43  	if f == 0 {
    44  		return 0, math.Copysign(1, g), g1
    45  	}
    46  
    47  	const safmin = dlamchS
    48  	const safmax = 1 / safmin
    49  	rtmin := math.Sqrt(safmin)
    50  	rtmax := math.Sqrt(safmax / 2)
    51  
    52  	f1 := math.Abs(f)
    53  
    54  	if rtmin < f1 && f1 < rtmax && rtmin < g1 && g1 < rtmax {
    55  		d := math.Sqrt(f*f + g*g)
    56  		cs = f1 / d
    57  		r = math.Copysign(d, f)
    58  		sn = g / r
    59  
    60  		return cs, sn, r
    61  	}
    62  
    63  	u := math.Min(math.Max(safmin, math.Max(f1, g1)), safmax)
    64  	fs := f / u
    65  	gs := g / u
    66  	d := math.Sqrt(fs*fs + gs*gs)
    67  	cs = math.Abs(fs) / d
    68  	r = math.Copysign(d, f)
    69  	sn = gs / r
    70  	r *= u
    71  
    72  	return cs, sn, r
    73  }