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) 'uv
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 + 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(,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