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.
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).
16 func (impl Implementation) Dlanv2(a, b, c, d float64) (aa, bb, cc, dd float64, rt1r, rt1i, rt2r, rt2i float64, cs, sn float64) {
18 case c == 0: // Matrix is already upper triangular.
25 case b == 0: // Matrix is lower triangular, swap rows and columns.
32 case a == d && math.Signbit(b) != math.Signbit(c): // Matrix is already in the standard Schur form.
42 bcmax := math.Max(math.Abs(b), math.Abs(c))
43 bcmis := math.Min(math.Abs(b), math.Abs(c))
47 scale := math.Max(math.Abs(p), bcmax)
48 z := p/scale*p + bcmax/scale*bcmis
52 // Real eigenvalues. Compute aa and dd.
54 z = p + math.Sqrt(scale)*math.Sqrt(z)
56 z = p - math.Sqrt(scale)*math.Sqrt(z)
59 dd = d - bcmax/z*bcmis
60 // Compute bb and the rotation matrix.
61 tau := impl.Dlapy2(c, z)
67 // Complex eigenvalues, or real (almost) equal eigenvalues.
68 // Make diagonal elements equal.
70 tau := impl.Dlapy2(sigma, temp)
71 cs = math.Sqrt((1 + math.Abs(sigma)/tau) / 2)
76 // Compute [ aa bb ] = [ a b ] [ cs -sn ]
77 // [ cc dd ] [ c d ] [ sn cs ]
82 // Compute [ a b ] = [ cs sn ] [ aa bb ]
83 // [ c d ] [-sn cs ] [ cc dd ]
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))
106 tau = 1 / math.Sqrt(math.Abs(bb+cc))
113 cs, sn = cs*cs1-sn*sn1, cs*sn1+sn+cs1
124 // Store eigenvalues in (rt1r,rt1i) and (rt2r,rt2i).
128 rt1i = math.Sqrt(math.Abs(bb)) * math.Sqrt(math.Abs(cc))