hqr.txt

用VB实现矩阵特征值的求解算法集。包括：对称矩阵的雅可比变换

Sub HQR(A(), N, WR(), WI())
    ANORM = Abs(A(1, 1))
    For I = 2 To N

        For J = I - 1 To N
            ANORM = ANORM + Abs(A(I, J))
        Next J
    Next I
    NN = N
    T = 0#
1   If NN >= 1 Then
        ITS = 0
2       For L = NN To 2 Step -1
            S = Abs(A(L - 1, L - 1)) + Abs(A(L, L))
            If S = 0# Then S = ANORM
            If Abs(A(L, L - 1)) + S = S Then GoTo 3
        Next L
        L = 1
3       X = A(NN, NN)
        If L = NN Then
            WR(NN) = X + T
            WI(NN) = 0#
            NN = NN - 1
        Else
            Y = A(NN - 1, NN - 1)
            W = A(NN, NN - 1) * A(NN - 1, NN)
            If L = NN - 1 Then
                P = 0.5 * (Y - X)
                Q = P ^ 2 + W
                Z = Sqr(Abs(Q))
                X = X + T
                If Q >= 0# Then
                    Z = P + Abs(Z) * Sgn(P)
                    WR(NN) = Z + X
                    WR(NN - 1) = WR(NN)
                    If Z <> 0# Then WR(NN) = X - W / Z
                    WI(NN) = 0#
                    WI(NN - 1) = 0#
                Else
                    WR(NN) = X + P
                    WR(NN - 1) = WR(NN)
                    WI(NN) = Z
                    WI(NN - 1) = -Z
                End If
                NN = NN - 2
            Else
                If ITS = 30 Then Print " too many iterations "
                If ITS = 10 Or ITS = 20 Then
                    T = T + X
                    For I = 1 To NN
                        A(I, I) = A(I, I) - X
                    Next I
                    S = Abs(A(NN, NN - 1)) + Abs(A(NN - 1, NN - 2))
                    X = 0.75 * S
                    Y = X
                    W = -0.4375 * S ^ 2
                End If
                ITS = ITS + 1
                For M = NN - 2 To L Step -1
                    Z = A(M, M)
                    R = X - Z
                    S = Y - Z
                    P = (R * S - W) / A(M + 1, M) + A(M, M + 1)
                    Q = A(M + 1, M + 1) - Z - R - S
                    R = A(M + 2, M + 1)
                    S = Abs(P) + Abs(Q) + Abs(R)
                    P = P / S
                    Q = Q / S
                    R = R / S
                    If M = L Then GoTo 4
                    U = Abs(A(M, M - 1)) * (Abs(Q) + Abs(R))
                    BBB = Abs(A(M + 1, M + 1))
                    AAA = Abs(A(M - 1, M - 1)) + Abs(Z) + BBB
                    V = Abs(P) * AAA
                    If U + V = V Then GoTo 4
                Next M
4               For I = M + 2 To NN
                    A(I, I - 2) = 0#
                    If I <> M + 2 Then A(I, I - 3) = 0#
                Next I
                For K = M To NN - 1
                    If K <> M Then
                        P = A(K, K - 1)
                        Q = A(K + 1, K - 1)
                        R = 0#
                        If K <> NN - 1 Then R = A(K + 2, K - 1)
                        X = Abs(P) + Abs(Q) + Abs(R)
                        If X <> 0# Then
                            P = P / X
                            Q = Q / X
                            R = R / X
                        End If
                    End If
                    S = Sqr(P ^ 2 + Q ^ 2 + R ^ 2) * Sgn(P)
                    If S <> 0# Then
                        If K = M Then
                            If L <> M Then A(K, K - 1) = -A(K, K - 1)
                        Else
                            A(K, K - 1) = -S * X
                        End If
                        P = P + S
                        X = P / S
                        Y = Q / S
                        Z = R / S
                        Q = Q / P
                        R = R / P
                        For J = K To NN
                            P = A(K, J) + Q * A(K + 1, J)
                            If K <> NN - 1 Then
                                P = P + R * A(K + 2, J)
                                A(K + 2, J) = A(K + 2, J) - P * Z
                            End If
                            A(K + 1, J) = A(K + 1, J) - P * Y
                            A(K, J) = A(K, J) - P * X
                        Next J
                        If NN > K + 3 Then
                           III = K + 3
                        Else
                           III = NN
                        End If
                        For I = L To III
                            P = X * A(I, K) + Y * A(I, K + 1)
                            If K <> NN - 1 Then
                                P = P + Z * A(I, K + 2)
                                A(I, K + 2) = A(I, K + 2) - P * R
                            End If
                            A(I, K + 1) = A(I, K + 1) - P * Q
                            A(I, K) = A(I, K) - P
                        Next I
                    End If
                Next K
                GoTo 2
            End If
        End If
        GoTo 1
    End If
End Sub