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
5package gonum
6
7import "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.
16func (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}
116