github.com/gonum/lapack@v0.0.0-20181123203213-e4cdc5a0bff9/native/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 native
     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  }