blob: f9b5f6f8d7b47a91f54243a3137b25c9fbb63c9d [file] [log] [blame]
// Copyright ©2016 The Gonum Authors. All rights reserved.
// Use of this source code is governed by a BSD-style
// license that can be found in the LICENSE file.
package gonum
import "math"
// Dlanv2 computes the Schur factorization of a real 2×2 matrix:
// [ a b ] = [ cs -sn ] * [ aa bb ] * [ cs sn ]
// [ c d ] [ sn cs ] [ cc dd ] * [-sn cs ]
// If cc is zero, aa and dd are real eigenvalues of the matrix. Otherwise it
// holds that aa = dd and bb*cc < 0, and aa ± sqrt(bb*cc) are complex conjugate
// eigenvalues. The real and imaginary parts of the eigenvalues are returned in
// (rt1r,rt1i) and (rt2r,rt2i).
func (impl Implementation) Dlanv2(a, b, c, d float64) (aa, bb, cc, dd float64, rt1r, rt1i, rt2r, rt2i float64, cs, sn float64) {
switch {
case c == 0: // Matrix is already upper triangular.
aa = a
bb = b
cc = 0
dd = d
cs = 1
sn = 0
case b == 0: // Matrix is lower triangular, swap rows and columns.
aa = d
bb = -c
cc = 0
dd = a
cs = 0
sn = 1
case a == d && math.Signbit(b) != math.Signbit(c): // Matrix is already in the standard Schur form.
aa = a
bb = b
cc = c
dd = d
cs = 1
sn = 0
default:
temp := a - d
p := temp / 2
bcmax := math.Max(math.Abs(b), math.Abs(c))
bcmis := math.Min(math.Abs(b), math.Abs(c))
if b*c < 0 {
bcmis *= -1
}
scale := math.Max(math.Abs(p), bcmax)
z := p/scale*p + bcmax/scale*bcmis
eps := dlamchP
if z >= 4*eps {
// Real eigenvalues. Compute aa and dd.
if p > 0 {
z = p + math.Sqrt(scale)*math.Sqrt(z)
} else {
z = p - math.Sqrt(scale)*math.Sqrt(z)
}
aa = d + z
dd = d - bcmax/z*bcmis
// Compute bb and the rotation matrix.
tau := impl.Dlapy2(c, z)
cs = z / tau
sn = c / tau
bb = b - c
cc = 0
} else {
// Complex eigenvalues, or real (almost) equal eigenvalues.
// Make diagonal elements equal.
safmn2 := math.Pow(dlamchB, math.Log(dlamchS/dlamchE)/math.Log(dlamchB)/2)
safmx2 := 1 / safmn2
sigma := b + c
loop:
for iter := 0; iter < 20; iter++ {
scale = math.Max(math.Abs(temp), math.Abs(sigma))
switch {
case scale >= safmx2:
sigma *= safmn2
temp *= safmn2
case scale <= safmn2:
sigma *= safmx2
temp *= safmx2
default:
break loop
}
}
p = temp / 2
tau := impl.Dlapy2(sigma, temp)
cs = math.Sqrt((1 + math.Abs(sigma)/tau) / 2)
sn = -p / (tau * cs)
if sigma < 0 {
sn *= -1
}
// Compute [ aa bb ] = [ a b ] [ cs -sn ]
// [ cc dd ] [ c d ] [ sn cs ]
aa = a*cs + b*sn
bb = -a*sn + b*cs
cc = c*cs + d*sn
dd = -c*sn + d*cs
// Compute [ a b ] = [ cs sn ] [ aa bb ]
// [ c d ] [-sn cs ] [ cc dd ]
a = aa*cs + cc*sn
b = bb*cs + dd*sn
c = -aa*sn + cc*cs
d = -bb*sn + dd*cs
temp = (a + d) / 2
aa = temp
bb = b
cc = c
dd = temp
if cc != 0 {
if bb != 0 {
if math.Signbit(bb) == math.Signbit(cc) {
// Real eigenvalues, reduce to
// upper triangular form.
sab := math.Sqrt(math.Abs(bb))
sac := math.Sqrt(math.Abs(cc))
p = sab * sac
if cc < 0 {
p *= -1
}
tau = 1 / math.Sqrt(math.Abs(bb+cc))
aa = temp + p
bb = bb - cc
cc = 0
dd = temp - p
cs1 := sab * tau
sn1 := sac * tau
cs, sn = cs*cs1-sn*sn1, cs*sn1+sn*cs1
}
} else {
bb = -cc
cc = 0
cs, sn = -sn, cs
}
}
}
}
// Store eigenvalues in (rt1r,rt1i) and (rt2r,rt2i).
rt1r = aa
rt2r = dd
if cc != 0 {
rt1i = math.Sqrt(math.Abs(bb)) * math.Sqrt(math.Abs(cc))
rt2i = -rt1i
}
return
}