d8r1.frm

常用的数值算法的VB程序

VERSION 5.00
Begin VB.Form Form1 
   Caption         =   "Form1"
   ClientHeight    =   7590
   ClientLeft      =   2055
   ClientTop       =   540
   ClientWidth     =   5010
   LinkTopic       =   "Form1"
   ScaleHeight     =   7590
   ScaleWidth      =   5010
   Begin VB.CommandButton Command1 
      Caption         =   "Command1"
      Height          =   375
      Left            =   3480
      TabIndex        =   0
      Top             =   6960
      Width           =   1215
   End
End
Attribute VB_Name = "Form1"
Attribute VB_GlobalNameSpace = False
Attribute VB_Creatable = False
Attribute VB_PredeclaredId = True
Attribute VB_Exposed = False
Private Sub Command1_Click()
    'PROGRAM D8R1
    'Driver for routine JACOBI
    NP = 3
    Dim D(3), V(3, 3), R(3)
    Dim A(3, 3), E(3, 3)
    A(1, 1) = 1#: A(1, 2) = 2#: A(1, 3) = 3#
    A(2, 1) = 2#: A(2, 2) = 2#: A(2, 3) = 3#
    A(3, 1) = 3#: A(3, 2) = 3#: A(3, 3) = 3#
    For II = 1 To 3
        For JJ = 1 To 3
            E(II, JJ) = A(II, JJ)
        Next JJ
    Next II
    Call JACOBI(E(), NP, D(), V(), NROT)
    Print Tab(5); "Number of JACOBI rotations:  "; NROT
    Print Tab(5); " Eigenvalues:"
    For J = 1 To NP
        Print Tab(7); Format$(D(J), "#0.######")
    Next J
    Print
    Print Tab(5); " Eigenvectors:"
    For J = 1 To NP
        Print Tab(5); "Number  "; J
        For K = 1 To NP
            Print Tab(5 + (K - 1) * 15); Format$(V(K, J), "##.#####0");
        Next K
        Print
    Next J
    Print
    'Eigenvector test
    Print Tab(5); " Eigenvector Test"
    For J = 1 To NP
        For L = 1 To NP
            R(L) = 0#
            For K = 1 To NP
               If K > L Then
                    KK = L
                    LL = K
                Else
                    KK = K
                    LL = L
                End If
                R(L) = R(L) + A(LL, KK) * V(K, J)
            Next K
        Next L
        Print
        Print Tab(5); " Vector Number", J
        Print
        Print Tab(5); " Vector        Mtrx*Vec.        Ratio"
        For L = 1 To NP
            RATIO = R(L) / V(L, J)
            Print Tab(5); Format$(V(L, J), "#.#####0"),
            Print Tab(20); Format$(R(L), "#.#####0"),
            Print Tab(35); Format$(RATIO, "#.#####0")
        Next L
    Next J
End Sub
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