1// Copyright ©2016 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// Dlanv2 computes the Schur factorization of a real 2×2 matrix:
10//  [ a b ] = [ cs -sn ] * [ aa bb ] * [ cs sn ]
11//  [ c d ]   [ sn  cs ]   [ cc dd ] * [-sn cs ]
12// If cc is zero, aa and dd are real eigenvalues of the matrix. Otherwise it
13// holds that aa = dd and bb*cc < 0, and aa ± sqrt(bb*cc) are complex conjugate
14// eigenvalues. The real and imaginary parts of the eigenvalues are returned in
15// (rt1r,rt1i) and (rt2r,rt2i).
16func (impl Implementation) Dlanv2(a, b, c, d float64) (aa, bb, cc, dd float64, rt1r, rt1i, rt2r, rt2i float64, cs, sn float64) {
17	switch {
18	case c == 0: // Matrix is already upper triangular.
19		aa = a
20		bb = b
21		cc = 0
22		dd = d
23		cs = 1
24		sn = 0
25	case b == 0: // Matrix is lower triangular, swap rows and columns.
26		aa = d
27		bb = -c
28		cc = 0
29		dd = a
30		cs = 0
31		sn = 1
32	case a == d && math.Signbit(b) != math.Signbit(c): // Matrix is already in the standard Schur form.
33		aa = a
34		bb = b
35		cc = c
36		dd = d
37		cs = 1
38		sn = 0
39	default:
40		temp := a - d
41		p := temp / 2
42		bcmax := math.Max(math.Abs(b), math.Abs(c))
43		bcmis := math.Min(math.Abs(b), math.Abs(c))
44		if b*c < 0 {
45			bcmis *= -1
46		}
47		scale := math.Max(math.Abs(p), bcmax)
48		z := p/scale*p + bcmax/scale*bcmis
49		eps := dlamchP
50
51		if z >= 4*eps {
52			// Real eigenvalues. Compute aa and dd.
53			if p > 0 {
54				z = p + math.Sqrt(scale)*math.Sqrt(z)
55			} else {
56				z = p - math.Sqrt(scale)*math.Sqrt(z)
57			}
58			aa = d + z
59			dd = d - bcmax/z*bcmis
60			// Compute bb and the rotation matrix.
61			tau := impl.Dlapy2(c, z)
62			cs = z / tau
63			sn = c / tau
64			bb = b - c
65			cc = 0
66		} else {
67			// Complex eigenvalues, or real (almost) equal eigenvalues.
68			// Make diagonal elements equal.
69			sigma := b + c
70			tau := impl.Dlapy2(sigma, temp)
71			cs = math.Sqrt((1 + math.Abs(sigma)/tau) / 2)
72			sn = -p / (tau * cs)
73			if sigma < 0 {
74				sn *= -1
75			}
76			// Compute [ aa bb ] = [ a b ] [ cs -sn ]
77			//         [ cc dd ]   [ c d ] [ sn  cs ]
78			aa = a*cs + b*sn
79			bb = -a*sn + b*cs
80			cc = c*cs + d*sn
81			dd = -c*sn + d*cs
82			// Compute [ a b ] = [ cs sn ] [ aa bb ]
83			//         [ c d ]   [-sn cs ] [ cc dd ]
84			a = aa*cs + cc*sn
85			b = bb*cs + dd*sn
86			c = -aa*sn + cc*cs
87			d = -bb*sn + dd*cs
88
89			temp = (a + d) / 2
90			aa = temp
91			bb = b
92			cc = c
93			dd = temp
94
95			if cc != 0 {
96				if bb != 0 {
97					if math.Signbit(bb) == math.Signbit(cc) {
98						// Real eigenvalues, reduce to
99						// upper triangular form.
100						sab := math.Sqrt(math.Abs(bb))
101						sac := math.Sqrt(math.Abs(cc))
102						p = sab * sac
103						if cc < 0 {
104							p *= -1
105						}
106						tau = 1 / math.Sqrt(math.Abs(bb+cc))
107						aa = temp + p
108						bb = bb - cc
109						cc = 0
110						dd = temp - p
111						cs1 := sab * tau
112						sn1 := sac * tau
113						cs, sn = cs*cs1-sn*sn1, cs*sn1+sn+cs1
114					}
115				} else {
116					bb = -cc
117					cc = 0
118					cs, sn = -sn, cs
119				}
120			}
121		}
122	}
123
124	// Store eigenvalues in (rt1r,rt1i) and (rt2r,rt2i).
125	rt1r = aa
126	rt2r = dd
127	if cc != 0 {
128		rt1i = math.Sqrt(math.Abs(bb)) * math.Sqrt(math.Abs(cc))
129		rt2i = -rt1i
130	}
131	return
132}
133