VB.Net矩阵求特征值

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

推荐文章：那些年，做的几个应用