1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
|
// 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
}
|
