jacobi.txt

VB常用数值算法集2

Sub JACOBI(A(), N, D(), V(), NROT)
    Dim B(100), Z(100)
    For IP = 1 To N
        For IQ = 1 To N
            V(IP, IQ) = 0#
        Next IQ
        V(IP, IP) = 1#
    Next IP
    For IP = 1 To N
        B(IP) = A(IP, IP)
        D(IP) = B(IP)
        Z(IP) = 0#
    Next IP
    NROT = 0
    For I = 1 To 50
        SM = 0#
        For IP = 1 To N - 1
            For IQ = IP + 1 To N
                SM = SM + Abs(A(IP, IQ))
            Next IQ
        Next IP
        If SM = 0# Then Exit Sub
        If I < 4 Then
            TRESH = 0.2 * SM / N ^ 2
        Else
            TRESH = 0#
        End If
        For IP = 1 To N - 1
            For IQ = IP + 1 To N
                G = 100# * Abs(A(IP, IQ))
                SSS = Abs(D(IP)) + G
                DDD = Abs(D(IQ)) + G
                If I > 4 And SSS = Abs(D(IP)) And DDD = Abs(D(IQ)) Then
                    A(IP, IQ) = 0#
                ElseIf Abs(A(IP, IQ)) > TRESH Then
                    H = D(IQ) - D(IP)
                    If Abs(H) + G = Abs(H) Then
                        T = A(IP, IQ) / H
                    Else
                        THETA = 0.5 * H / A(IP, IQ)
                        T = 1# / (Abs(THETA) + Sqr(1# + THETA ^ 2))
                        If THETA < 0# Then T = -T
                    End If
                    C = 1# / Sqr(1# + T ^ 2)
                    S = T * C
                    TAU = S / (1# + C)
                    H = T * A(IP, IQ)
                    Z(IP) = Z(IP) - H
                    Z(IQ) = Z(IQ) + H
                    D(IP) = D(IP) - H
                    D(IQ) = D(IQ) + H
                    A(IP, IQ) = 0#
                    For J = 1 To IP - 1
                        G = A(J, IP)
                        H = A(J, IQ)
                        A(J, IP) = G - S * (H + G * TAU)
                        A(J, IQ) = H + S * (G - H * TAU)
                    Next J
                    For J = IP + 1 To IQ - 1
                        G = A(IP, J)
                        H = A(J, IQ)
                        A(IP, J) = G - S * (H + G * TAU)
                        A(J, IQ) = H + S * (G - H * TAU)
                    Next J
                    For J = IQ + 1 To N
                        G = A(IP, J)
                        H = A(IQ, J)
                        A(IP, J) = G - S * (H + G * TAU)
                        A(IQ, J) = H + S * (G - H * TAU)
                    Next J
                    For J = 1 To N
                        G = V(J, IP)
                        H = V(J, IQ)
                        V(J, IP) = G - S * (H + G * TAU)
                        V(J, IQ) = H + S * (G - H * TAU)
                    Next J
                    NROT = NROT + 1
                End If
            Next IQ
        Next IP
        For IP = 1 To N
            B(IP) = B(IP) + Z(IP)
            D(IP) = B(IP)
            Z(IP) = 0#
        Next IP
    Next I
    Print " 50 iterations should never happen"
End Sub