vendor/gonum.org/v1/gonum/lapack/gonum/dlasv2.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 // Dlasv2 computes the singular value decomposition of a 2×2 matrix.
  10 //  [ csl snl] [f g] [csr -snr] = [ssmax     0]
  11 //  [-snl csl] [0 h] [snr  csr] = [    0 ssmin]
  12 // ssmax is the larger absolute singular value, and ssmin is the smaller absolute
  13 // singular value. [cls, snl] and [csr, snr] are the left and right singular vectors.
  14 //
  15 // Dlasv2 is an internal routine. It is exported for testing purposes.
  16 func (impl Implementation) Dlasv2(f, g, h float64) (ssmin, ssmax, snr, csr, snl, csl float64) {
  17         ft := f
  18         fa := math.Abs(ft)
  19         ht := h
  20         ha := math.Abs(h)
  21         // pmax points to the largest element of the matrix in terms of absolute value.
  22         // 1 if F, 2 if G, 3 if H.
  23         pmax := 1
  24         swap := ha > fa
  25         if swap {
  26                 pmax = 3
  27                 ft, ht = ht, ft
  28                 fa, ha = ha, fa
  29         }
  30         gt := g
  31         ga := math.Abs(gt)
  32         var clt, crt, slt, srt float64
  33         if ga == 0 {
  34                 ssmin = ha
  35                 ssmax = fa
  36                 clt = 1
  37                 crt = 1
  38                 slt = 0
  39                 srt = 0
  40         } else {
  41                 gasmall := true
  42                 if ga > fa {
  43                         pmax = 2
  44                         if (fa / ga) < dlamchE {
  45                                 gasmall = false
  46                                 ssmax = ga
  47                                 if ha > 1 {
  48                                         ssmin = fa / (ga / ha)
  49                                 } else {
  50                                         ssmin = (fa / ga) * ha
  51                                 }
  52                                 clt = 1
  53                                 slt = ht / gt
  54                                 srt = 1
  55                                 crt = ft / gt
  56                         }
  57                 }
  58                 if gasmall {
  59                         d := fa - ha
  60                         l := d / fa
  61                         if d == fa { // deal with inf
  62                                 l = 1
  63                         }
  64                         m := gt / ft
  65                         t := 2 - l
  66                         s := math.Hypot(t, m)
  67                         var r float64
  68                         if l == 0 {
  69                                 r = math.Abs(m)
  70                         } else {
  71                                 r = math.Hypot(l, m)
  72                         }
  73                         a := 0.5 * (s + r)
  74                         ssmin = ha / a
  75                         ssmax = fa * a
  76                         if m == 0 {
  77                                 if l == 0 {
  78                                         t = math.Copysign(2, ft) * math.Copysign(1, gt)
  79                                 } else {
  80                                         t = gt/math.Copysign(d, ft) + m/t
  81                                 }
  82                         } else {
  83                                 t = (m/(s+t) + m/(r+l)) * (1 + a)
  84                         }
  85                         l = math.Hypot(t, 2)
  86                         crt = 2 / l
  87                         srt = t / l
  88                         clt = (crt + srt*m) / a
  89                         slt = (ht / ft) * srt / a
  90                 }
  91         }
  92         if swap {
  93                 csl = srt
  94                 snl = crt
  95                 csr = slt
  96                 snr = clt
  97         } else {
  98                 csl = clt
  99                 snl = slt
 100                 csr = crt
 101                 snr = srt
 102         }
 103         var tsign float64
 104         switch pmax {
 105         case 1:
 106                 tsign = math.Copysign(1, csr) * math.Copysign(1, csl) * math.Copysign(1, f)
 107         case 2:
 108                 tsign = math.Copysign(1, snr) * math.Copysign(1, csl) * math.Copysign(1, g)
 109         case 3:
 110                 tsign = math.Copysign(1, snr) * math.Copysign(1, snl) * math.Copysign(1, h)
 111         }
 112         ssmax = math.Copysign(ssmax, tsign)
 113         ssmin = math.Copysign(ssmin, tsign*math.Copysign(1, f)*math.Copysign(1, h))
 114         return ssmin, ssmax, snr, csr, snl, csl
 115 }