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