Public Function Math_Matrix_EigenValue(ByVal K1(,) As Single, ByVal n As Integer, ByVal LoopNumber As Integer, ByVal Errro As Int16, ByRef Ret(,) As Double) As Boolean 'ret里是n*2的数组,第一列是实数部分,第2列为虚数部分
Dim i As Integer = K1.Length / n
If i * n <> K1.Length Then
Return False
End If
Dim j As Integer
Dim k As Integer
Dim t As Integer
Dim m As Integer
Dim A(0, 0) As Single
ReDim Ret(n - 1, 1) 'u v
Dim erro As Double = Math.Pow(0.1, Errro)
Dim b As Single
Dim c As Single
Dim d As Single
Dim g As Single
Dim xy As Single
Dim p As Single
Dim q As Single
Dim r As Single
Dim x As Single
Dim s As Single
Dim e As Single
Dim f As Single
Dim z As Single
Dim y As Single
Dim loop1 As Integer = LoopNumber
Math_Matrix_Hessenberg(K1, n, A) '将方阵K1转化成上Hessenberg矩阵A
m = n
While m <> 0
t = m - 1
While t > 0
If Math.Abs(A(t, t - 1)) > erro * (Math.Abs(A(t - 1, t - 1)) + Math.Abs(A(t, t))) Then
t -= 1
Else
Exit While
End If
End While
If t = m - 1 Then
Ret(m - 1, 0) = A(m - 1, m - 1)
Ret(m - 1, 1) = 0
m -= 1
loop1 = LoopNumber
ElseIf t = m - 2 Then
b = -(A(m - 1, m - 1) + A(m - 2, m - 2))
c = A(m - 1, m - 1) * A(m - 2, m - 2) - A(m - 1, m - 2) * A(m - 2, m - 1)
d = b * b - 4 * c
y = Math.Pow(Math.Abs(d), 0.5)
If d > 0 Then
xy = 1
If b < 0 Then
xy = -1
End If
Ret(m - 1, 0) = -(b + xy * y) / 2
Ret(m - 1, 1) = 0
Ret(m - 2, 0) = c / Ret(m - 1, 0)
Ret(m - 2, 1) = 0
Else
Ret(m - 1, 0) = -b / 2
Ret(m - 2, 0) = Ret(m - 1, 0)
Ret(m - 1, 1) = y / 2
Ret(m - 2, 1) = -Ret(m - 1, 1)
End If
m -= 2
loop1 = LoopNumber
Else
If loop1 < 1 Then
Return False
End If
loop1 -= 1
j = t + 2
While j < m
A(j, j - 2) = 0
j += 1
End While
j = t + 3
While j < m
A(j, j - 3) = 0
j += 1
End While
k = t
While k < m - 1
If k <> t Then
p = A(k, k - 1)
q = A(k + 1, k - 1)
If k <> m - 2 Then
r = A(k + 2, k - 1)
Else
r = 0
End If
Else
b = A(m - 1, m - 1)
c = A(m - 2, m - 2)
x = b + c
y = c * b - A(m - 2, m - 1) * A(m - 1, m - 2)
p = A(t, t) * (A(t, t) - x) + A(t, t + 1) * A(t + 1, t) + y
q = A(t + 1, t) * (A(t, t) + A(t + 1, t + 1) - x)
r = A(t + 1, t) * A(t + 2, t + 1)
End If
If p <> 0 Or q <> 0 Or r <> 0 Then
If p < 0 Then
xy = -1
Else
xy = 1
End If
s = xy * Math.Pow(p * p + q * q + r * r, 0.5)
If k <> t Then
A(k, k - 1) = -s
End If
e = -q / s
f = -r / s
x = -p / s
y = -x - f * r / (p + s)
g = e * r / (p + s)
z = -x - e * q / (p + s)
For j = k To m - 1
b = A(k, j)
c = A(k + 1, j)
p = x * b + e * c
q = e * b + y * c
r = f * b + g * c
If k <> m - 2 Then
b = A(k + 2, j)
p += f * b
q += g * b
r += z * b
A(k + 2, j) = r
End If
A(k + 1, j) = q
A(k, j) = p
Next
j = k + 3
If j >= m - 1 Then
j = m - 1
End If
For i = t To j
b = A(i, k)
c = A(i, k + 1)
p = x * b + e * c
q = e * b + y * c
r = f * b + g * c
If k <> m - 2 Then
b = A(i, k + 2)
p += f * b
q += g * b
r += z * b
A(i, k + 2) = r
End If
A(i, k + 1) = q
A(i, k) = p
Next
End If
k += 1
End While
End If
End While
Return True
End Function
Public Function Math_Matrix_Hessenberg(ByVal A(,) As Single, ByVal n As Integer, ByRef ret(,) As Single) As Integer
Dim i As Integer
Dim j As Integer
Dim k As Integer
Dim temp As Single
Dim MaxNumber As Integer
n -= 1
ReDim ret(n, n)
For k = 1 To n - 1
i = k - 1
MaxNumber = k
temp = Math.Abs(A(k, i))
For j = k + 1 To n
If Math.Abs(A(j, i)) > temp Then
MaxNumber = j
End If
Next
ret(0, 0) = A(MaxNumber, i) '储存最大值
i = MaxNumber
If ret(0, 0) <> 0 Then
If i <> k Then
For j = k - 1 To n
temp = A(i, j)
A(i, j) = A(k, j)
A(k, j) = temp
Next
For j = 0 To n
temp = A(j, i)
A(j, i) = A(j, k)
A(j, k) = temp
Next
End If
For i = k + 1 To n
temp = A(i, k - 1) / ret(0, 0)
A(i, k - 1) = 0
For j = k To n
A(i, j) -= temp * A(k, j)
Next
For j = 0 To n
A(j, k) += temp * A(j, i)
Next
Next
End If
Next
For i = 0 To n
For j = 0 To n
ret(i, j) = A(i, j)
Next
Next
Return n + 1
End Function