d1r9.frm

矩阵特征值的求解过程之二

VERSION 5.00
Begin VB.Form Form1 
   Caption         =   "Form1"
   ClientHeight    =   6255
   ClientLeft      =   60
   ClientTop       =   345
   ClientWidth     =   6240
   LinkTopic       =   "Form1"
   ScaleHeight     =   6255
   ScaleWidth      =   6240
   StartUpPosition =   3  'Windows Default
   Begin VB.CommandButton Command1 
      Caption         =   "Command1"
      Height          =   300
      Left            =   4560
      TabIndex        =   0
      Top             =   5760
      Width           =   1335
   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 D1R9
    'Driver program for routine SPARSE
    N = 20
    Dim B(20), X(20), BCMP(20)
    '输入已知的方程组的右端向量B
    For I = 1 To N
        X(I) = 0#
        B(I) = 1#
    Next I
    B(1) = 3#
    B(N) = -1#
    '利用下列子程序ASUB输入已知的方程组的系数矩阵
    Call SPARSE(B(), N, X(), RSQ)
    Print
    Print Tab(5); "Sum - squared residual:", RSQ
    Print
    Print Tab(5); "计算出的方程组的解"
    For I = 1 To 4
        For J = I * 5 - 4 To I * 5
            Print Tab((J - (I - 1) * 5) * 9); Format$(X(J), "#.####E+00");
        Next J
        Print
    Next I
    '将计算出的方程组的解乘以系数矩阵，以验证计算结果正确
    Call ASUB(X(), BCMP())
    Print
    Print Tab(5); "计算出的方程组的解乘以系数矩阵的结果"
    Print
    Print Tab(7); "解乘以系数矩阵"; Tab(25); "方程组的右端向量"
    For I = 1 To N
        Print Tab(13); Format$(BCMP(I), "##.##"); Tab(32); Format$(B(I), "##.##")
    Next I
End Sub
Sub ASUB(XIN(), XOUT())
    N = 20
    XOUT(1) = XIN(1) + 2# * XIN(2)
    XOUT(N) = -2# * XIN(N - 1) + XIN(N)
    For I = 2 To N - 1
        XOUT(I) = -2# * XIN(I - 1) + XIN(I) + 2# * XIN(I + 1)
    Next I
End Sub
Sub ATSUB(XIN(), XOUT())
    N = 20
    XOUT(1) = XIN(1) - 2# * XIN(2)
    XOUT(N) = 2# * XIN(N - 1) + XIN(N)
    For I = 2 To N - 1
        XOUT(I) = 2# * XIN(I - 1) + XIN(I) - 2# * XIN(I + 1)
    Next I
End Sub
Sub SPARSE(B(), N, X(), RSQ)
    NMAX = 500
    EPS = 0.000001
    Dim G(500), H(500), XI(500), XJ(500)
    EPS2 = N * EPS ^ 2
    IRST = 0
1   IRST = IRST + 1
    Call ASUB(X(), XI())
    RP = 0#
    BSQ = 0#
    For J = 1 To N
        BSQ = BSQ + B(J) ^ 2
        XI(J) = XI(J) - B(J)
        RP = RP + XI(J) ^ 2
    Next J
    Call ATSUB(XI(), G())
    For J = 1 To N
        G(J) = -G(J)
        H(J) = G(J)
    Next J
    For ITER = 1 To 10 * N
        Call ASUB(H(), XI())
        ANUM = 0#
        ADEN = 0#
        For J = 1 To N
            ANUM = ANUM + G(J) * H(J)
            ADEN = ADEN + XI(J) ^ 2
        Next J
        If ADEN = 0# Then Print "very singular matrix"
        ANUM = ANUM / ADEN
        For J = 1 To N
            XI(J) = X(J)
            X(J) = X(J) + ANUM * H(J)
        Next J
        Call ASUB(X(), XJ())
        RSQ = 0#
        For J = 1 To N
            XJ(J) = XJ(J) - B(J)
            RSQ = RSQ + XJ(J) ^ 2
        Next J
        If RSQ = RP Or RSQ <= BSQ * EPS2 Then Exit Sub
        If RSQ > RP Then
            For J = 1 To N
                X(J) = XI(J)
            Next J
            If IRST >= 3 Then
                Exit Sub
            End If
            GoTo 1
        End If
        RP = RSQ
        Call ATSUB(XJ(), XI())
        GG = 0#
        DGG = 0#
        For J = 1 To N
            GG = GG + G(J) ^ 2
            DGG = DGG + (XI(J) + G(J)) * XI(J)
        Next J
        If GG = 0# Then Exit Sub
        GAM = DGG / GG
        For J = 1 To N
            G(J) = -XI(J)
            H(J) = G(J) + GAM * H(J)
        Next J
    Next ITER
    Print "too many iterations"
End Sub