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) '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

（编辑：李大同）

【声明】本站内容均来自网络，其相关言论仅代表作者个人观点，不代表本站立场。若无意侵犯到您的权利，请及时与联系站长删除相关内容!