mdladjust.bas

采用VB编写

    Row = UBound(A, 1) - LBound(A, 1) + 1
    Col = UBound(A, 2) - LBound(A, 2) + 1
    If Row <> Col Then
        MsgBox "方程组的系数矩阵有误！"
        Exit Sub
    End If
    '准备约化过程的变量和数组
    n = UBound(b) - LBound(b) + 1
    If n <> Row Then
        MsgBox "方程组的系数矩阵与常数项大小不符！"
        Exit Sub
    End If
    ReDim L(2 To Row) As Double
    Dim sumAX As Double, iPos%, temp#
    '约化过程
    For iStep = 1 To n - 1
        '列选主元
        iPos = 0
        For iRow = iStep + 1 To n
            If Abs(A(iRow, iStep)) > Abs(A(iStep, iStep)) Then
                iPos = iRow
            End If
        Next iRow
        If iPos > iStep Then    '需要换主元
            For iCol = iStep To n
                temp = A(iStep, iCol)
                A(iStep, iCol) = A(iPos, iCol)
                A(iPos, iCol) = temp
            Next iCol
            temp = b(iStep)
            b(iStep) = b(iPos)
            b(iPos) = temp
        End If
        '约化过程
        For iRow = iStep + 1 To n
            L(iRow) = A(iRow, iStep) / A(iStep, iStep)
            For iCol = iStep To n
                A(iRow, iCol) = A(iRow, iCol) - L(iRow) * A(iStep, iCol)
            Next iCol
            b(iRow) = b(iRow) - L(iRow) * b(iStep)
        Next iRow
        ShowMatrix A
    Next iStep
    '回代过程
    X(n) = b(n) / A(n, n)
    For iRow = n - 1 To 1 Step -1
        sumAX = 0
        For iCol = n To iRow + 1 Step -1
            sumAX = sumAX + A(iRow, iCol) * X(iCol)
        Next iCol
        X(iRow) = (b(iRow) - sumAX) / A(iRow, iRow)
    Next iRow
End Sub

'Guass－Seidel迭代法求解线性方程组
Private Function Seidel(A, b, X, eps#) As Boolean
    Dim i%, j%
    Dim P#, Q#, s#, t#
    Dim Row%, Col%, n%
    Row = UBound(A, 1) - LBound(A, 1) + 1
    Col = UBound(A, 2) - LBound(A, 2) + 1
    n = UBound(b) - LBound(b) + 1
    If n <> Row Then
        MsgBox "方程组的系数矩阵与常数项大小不符！"
        Exit Function
    End If
    
    For i = 1 To n
        P = 0#
        X(i) = 0#
        For j = 1 To n
            If i <> j Then P = P + Abs(A(i, j))
        Next j
        If P >= Abs(A(i, i)) Then
            Seidel = False
            Exit Function
        End If
    Next i
    
    P = eps + 1#
    While P >= eps
        P = 0#
        For i = 1 To n
            t = X(i)
            s = 0#
            For j = 1 To n
                If j <> i Then s = s + A(i, j) * X(j)
            Next j
            
            X(i) = (b(i) - s) / (A(i, i))
            Q = Abs(X(i) - t) '/ (1# + Abs(x(i)))
            If Q > P Then P = Q
        Next i
    Wend
    
    Seidel = True
End Function

Public Sub ShowMatrix(tt)
    Dim i%, j%, n%, m%
    m = UBound(tt, 1) - LBound(tt, 1) + 1
    n = UBound(tt, 2) - LBound(tt, 2) + 1
    For i = 1 To m
        For j = 1 To n
            Debug.Print tt(i, j),
        Next j
        Debug.Print
    Next i
End Sub

'通用的间接平差解算过程：输入系数矩阵A、权矩阵P、常数向量L和解向量X，求出X，并通过参数传出去
Public Sub InAdjust(A, P, L, X)
    Dim a1%, a2%, p1%, p2%, L1%, x1%    '输入矩阵或向量的大小
    Dim At() As Double, AtP() As Double, Naa#(), W() As Double  '几个中间矩阵
    
    '计算并检查输入矩阵或向量的大小
    On Error Resume Next
    a1 = UBound(A, 1) - LBound(A, 1) + 1
    If Err Then
        MsgBox "系数矩阵A大小错误！"
        Exit Sub
    End If
    On Error Resume Next
    a2 = UBound(A, 2) - LBound(A, 2) + 1
    If Err Then
        MsgBox "系数矩阵A大小错误！"
        Exit Sub
    End If
    On Error Resume Next
    L1 = UBound(L) - LBound(L) + 1
    If Err Then
        MsgBox "常数向量L大小错误！"
        Exit Sub
    End If
    On Error Resume Next
    x1 = UBound(X) - LBound(X) + 1
    If Err Then
        MsgBox "解向量X大小错误！"
        Exit Sub
    End If
    On Error Resume Next
    p1 = UBound(P, 1) - LBound(P, 1) + 1
    If Err Then
        MsgBox "权矩阵P大小错误！"
        Exit Sub
    End If
    On Error Resume Next
    p2 = UBound(P, 2) - LBound(P, 2) + 1
    If Err Then
        MsgBox "权矩阵P大小错误！"
        Exit Sub
    End If
    If p1 <> p2 Then
        MsgBox "权矩阵P不是方阵！"
        Exit Sub
    End If
    If p1 <> a1 Or p2 <> a1 Then
        MsgBox "权矩阵P与系数矩阵A大小不符！"
        Exit Sub
    End If
    If a2 <> x1 Then
        MsgBox "系数矩阵A大小与解向量X大小不符！"
        Exit Sub
    End If
    If a1 <> L1 Then
        MsgBox "系数矩阵A大小与常数向量L大小不符！"
        Exit Sub
    End If
    
    '定义中间矩阵的大小
    ReDim At(1 To a2, 1 To a1), AtP(1 To a2, 1 To a1)
    ReDim Naa(1 To a2, 1 To a2), W(1 To a2)
    
    '组成法方程并计算
    Debug.Print "The A matrix is:"
    ShowMatrix A
    MatrixTrans A, At               '求A的转置矩阵
    Debug.Print "The At matrix is:"
    ShowMatrix At
    
    Debug.Print "The P matrix is:"
    ShowMatrix P
    Matrix_Multy AtP, At, P         '求AtP
    Debug.Print "and The AtP matrix is:"
    ShowMatrix AtP
    
    Matrix_Multy Naa, AtP, A        '法方程系数矩阵
    Debug.Print "the Naa matrix is:"
    ShowMatrix Naa
    
    Debug.Print "the L matrix is:"
    For x1 = LBound(L) To UBound(L)
        Debug.Print L(x1)
    Next x1
    Matrix_Multy W, AtP, L          '法方程常数向量
    Debug.Print "the W matrix is:"
    For x1 = LBound(W) To UBound(W)
        Debug.Print W(x1)
    Next x1
    
    MajorInColGuass Naa, W, X
    Debug.Print "the X matrix is:"
    For x1 = LBound(X) To UBound(X)
        Debug.Print X(x1)
    Next x1
    'Seidel Naa, W, x, 0.000001
End Sub


'通用的条件平差解算过程：输入系数矩阵A、权矩阵P、常数向量L和解向量X，求出X，并通过参数传出去
Public Sub CondiAdjust(b, P, W, V)
    Dim b1%, b2%, p1%, p2%, w1%, v1%    '输入矩阵或向量的大小
    Dim Q#(), Bt#(), QBt#(), Nbb#(), K#(), i% '几个中间矩阵
    
    '计算并检查输入矩阵或向量的大小
    On Error Resume Next
    b1 = UBound(b, 1) - LBound(b, 1) + 1
    If Err Then
        MsgBox "系数矩阵B大小错误！"
        Exit Sub
    End If
    On Error Resume Next
    b2 = UBound(b, 2) - LBound(b, 2) + 1
    If Err Then
        MsgBox "系数矩阵B大小错误！"
        Exit Sub
    End If
    On Error Resume Next
    w1 = UBound(W) - LBound(W) + 1
    If Err Then
        MsgBox "常数向量W大小错误！"
        Exit Sub
    End If
    On Error Resume Next
    v1 = UBound(V) - LBound(V) + 1
    If Err Then
        MsgBox "改正数向量V大小错误！"
        Exit Sub
    End If
    On Error Resume Next
    p1 = UBound(P, 1) - LBound(P, 1) + 1
    If Err Then
        MsgBox "权矩阵P大小错误！"
        Exit Sub
    End If
    On Error Resume Next
    p2 = UBound(P, 2) - LBound(P, 2) + 1
    If Err Then
        MsgBox "权矩阵P大小错误！"
        Exit Sub
    End If
    If p1 <> p2 Then
        MsgBox "权矩阵P不是方阵！"
        Exit Sub
    End If
    If p1 <> b2 Then
        MsgBox "权矩阵P与系数矩阵A大小不符！"
        Exit Sub
    End If
    If b2 <> v1 Then
        MsgBox "系数矩阵B大小与解向量V大小不符！"
        Exit Sub
    End If
    If b1 <> w1 Then
        MsgBox "系数矩阵B大小与常数向量W大小不符！"
        Exit Sub
    End If
    
    '定义中间矩阵的大小
    ReDim Bt(1 To b2, 1 To b1), QBt(1 To b2, 1 To b1)
    ReDim Nbb(1 To b1, 1 To b1), K(1 To b1), Q(1 To p1, 1 To p2)
    
    '组成法方程并计算
    For i = 1 To p1             '求Q矩阵
        Q(i, i) = 1 / P(i, i)
    Next i
    MatrixTrans b, Bt
    Matrix_Multy QBt, Q, Bt
    Matrix_Multy Nbb, b, QBt        '法方程系数矩阵
    ShowMatrix Nbb
    MajorInColGuass Nbb, W, K       '解法方程
    'Seidel Nbb, W, K, 0.0000001
    Matrix_Multy V, QBt, K          '求改正数
End Sub