lu

LU分解法求解线性方程组的VB.NET源码.里面有详细的使用说明.

    REM 用LU分解法求解线性方程组的解
    REM A(,)实型二组数组,输入、输出参数，在LUDCMP中，输入时按列存放实方阵A。输出时，对角线以下部分存放单位下
    REM 三解矩阵L，对角组及其以上部分存放上三角矩阵U，在LUBSKSB中，将LUDCMP中输出结果A作为输入

    REM N 整型变量，输入参数，方程系数矩阵A的阶数

    REM INDX（），整形一维数组，在子过程LUDCMP中为输出参数，用于记录置换矩阵，称为置换向量。在子过程LUBKSB中
    REM 为输入参数，输入子过程LUDCMP的输出结果

    REM 参数D，输出参数，依赖于行交换次数为偶（+1）还是奇（-1）

    REM 参数B（），一维实型数组，输入、输出参数。输入实向量b，输出时，方程组的解x存储在数组B（）中。

    Public Sub LUDCMP(ByRef A(,) As Double, ByVal N As Integer, ByRef INDX() As Integer, ByVal D As Integer)
        Dim NMAX As Integer = 100
        Dim TINY As Double = 1.0E-20
        Dim VV(100) As Double
        Dim i, j, k, imax As Integer
        Dim aamax, sum, dum As Double
        D = 1.0#
        For i = 1 To N
            aamax = 0.0#
            For j = 1 To N
                If Math.Abs(A(i, j)) > aamax Then aamax = Math.Abs(A(i, j))
            Next j
            If aamax = 0.0# Then Print("Singular matrix.")
            VV(i) = 1.0# / 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.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) * Math.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.0# Then A(j, j) = TINY
                dum = 1.0# / 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.0# Then A(N, N) = TINY
    End Sub

    Public Sub LUBKSB(ByRef A(,) As Double, ByVal N As Integer, ByRef INDX() As Integer, ByRef B() As Double)
        Dim II As Integer = 0
        Dim sum, LL As Double
        Dim i, j As Integer
        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.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