vendor/gonum.org/v1/gonum/lapack/gonum/dlartg.go

   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 // 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 }