求教，lu 分解法的vb 实现

Question

求教，lu 分解法的vb 实现

Answer 1

凑热闹的：
http://www.icl.miyazaki-u.ac.jp/k_aoyama/fpro/lu.html
http://www.kzm.info.gifu-u.ac.jp/~k03kimura/numerical_calculation.html

Answer 2

建议用matlab做吧，一条指令就搞定了
[L,U,P]=lu(A)

Answer 3

找到了VB的算法。
包括LU分解，倒推过程，及解方程组的过程。'LU分解
'A中保存的是L+U，为了节省内存，因为A中的前面一个元素用过后就可以不用了。
Sub LUDCMP(A(), N, INDX(), D)
    NMAX = 100
    TINY = 1E-20
    Dim VV(100)
    D = 1#
    For i = 1 To N
        AAMAX = 0#
        For j = 1 To N
            If Abs(A(i, j)) > AAMAX Then AAMAX = Abs(A(i, j))
        Next j
        If AAMAX = 0# Then Print "Singular matrix."
        VV(i) = 1# / AAMAX
    Next i
    For j = 1 To N
        If j > 1 Then
            For i = 1 To j - 1
                Sum = A(i, j)
                If i > 1 Then
                    For K = 1 To i - 1
                        Sum = Sum - A(i, K) * A(K, j)
                    Next K
                    A(i, j) = Sum
                End If
            Next i
        End If
        AAMAX = 0#
        For i = j To N
            Sum = A(i, j)
            If j > 1 Then
                For K = 1 To j - 1
                    Sum = Sum - A(i, K) * A(K, j)
                Next K
                A(i, j) = Sum
            End If
            DUM = VV(i) * Abs(Sum)
            If DUM >= AAMAX Then
                IMAX = i
                AAMAX = DUM
            End If
        Next i
        If j <> IMAX Then
            For K = 1 To N
                DUM = A(IMAX, K)
                A(IMAX, K) = A(j, K)
                A(j, K) = DUM
            Next K
            D = -D
            VV(IMAX) = VV(j)
        End If
        INDX(j) = IMAX
        If j <> N Then
            If A(j, j) = 0# Then A(j, j) = TINY
            DUM = 1# / A(j, j)
            For i = j + 1 To N
                A(i, j) = A(i, j) * DUM
            Next i
        End If
    Next j
    If A(N, N) = 0# Then A(N, N) = TINY
End Sub
'倒推过程，参看数值分析
Sub LUBKSB(A(), N, INDX(), B())
    II = 0
    For i = 1 To N
        LL = INDX(i)
        Sum = B(LL)
        B(LL) = B(i)
        If II <> 0 Then
            For j = II To i - 1
                Sum = Sum - A(i, j) * B(j)
            Next j
        ElseIf Sum <> 0# Then
            II = i
        End If
        B(i) = Sum
    Next i
    For i = N To 1 Step -1
        Sum = B(i)
        If i < N Then
            For j = i + 1 To N
                Sum = Sum - A(i, j) * B(j)
            Next j
        End If
        B(i) = Sum / A(i, i)
    Next i
End Sub
Private Sub Command1_Click()
    'program D1R2
    'Driver program for routine LUBKSB,LUDCMP
    N = 3
    Dim A(3, 3), B(3), A1(3, 3), INDX(3), X(3)
    '输入已知的方程组的系数矩阵
    A(1, 1) = 1#: A(1, 2) = 2#: A(1, 3) = 3#
    A(2, 1) = 3#: A(2, 2) = 1#: A(2, 3) = 5#
    A(3, 1) = 2#: A(3, 2) = 5#: A(3, 3) = 2#
    '输入已知的方程组的右端向量
    B(1) = 1#
    B(2) = 2#
    B(3) = 3#
    Print
    Print Tab(5); "已知的方程组的右端向量"
    Print Tab(14); Format$(B(1), "##.##")
    Print Tab(14); Format$(B(2), "##.##")
    Print Tab(14); Format$(B(3), "##.##")    Call LUDCMP(A1(), N, INDX(), P)
   'A1中的元素就是有L+U（相加）对应的元素
    
    For L = 1 To N
        X(L) = B(L)
    Next L
    Call LUBKSB(A1(), N, INDX(), X())
    Print
    Print Tab(5); "计算出的方程组的解"
    Print Tab(14); Format$(X(1), "#.####E+00")
    Print Tab(14); Format$(X(2), "#.####E+00")
    Print Tab(14); Format$(X(3), "#.####E+00")
    '将计算出的B乘以系数矩阵，以验证计算结果正确
    For L = 1 To N
        B(L) = 0#
        For j = 1 To N
            B(L) = B(L) + A(L, j) * X(j)
        Next j
    Next L
    Print
    Print Tab(5); "计算出的解乘以系数矩阵的结果"
    Print Tab(14); Format$(B(1), "##.##")
    Print Tab(14); Format$(B(2), "##.##")
    Print Tab(14); Format$(B(3), "##.##")
End Sub

Answer 4

哎呀，上面的
    Print Tab(5); "已知的方程组的右端向量"
    Print Tab(14); Format$(B(1), "##.##")
    Print Tab(14); Format$(B(2), "##.##")
    Print Tab(14); Format$(B(3), "##.##")
后面少了一段，请加上
   For i = 1 To N
        For j = 1 To N
            A1(i, j) = A(i, j)
        Next j
    Next i