mdladjust.bas

在测量平差中

    '准备约化过程的变量和数组
    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
        Debug.Print "after the "; iStep; "'s Guassion, the A is: "
        ShowMatrix A
        Debug.Print "The W matrix is:"
        For iRow = 1 To n
            Debug.Print B(iRow)
        Next iRow
    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, 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)
    
    '组成法方程并计算
    MatrixTrans A, At
    Matrix_Multy AtP, At, P
    Matrix_Multy Naa, AtP, A        '法方程系数矩阵
    Matrix_Multy W, AtP, L          '法方程常数向量
    
    MajorInColGuass Naa, W, x
    '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)
    
    '组成法方程并计算
    Debug.Print "The P matrix is:"
    ShowMatrix P
    For i = 1 To p1             '求Q矩阵
        Q(i, i) = 1 / P(i, i)
    Next i
    Debug.Print "and the Q matrix is:"
    ShowMatrix Q
    
    Debug.Print "The B matrix is:"
    ShowMatrix B
    MatrixTrans B, Bt           '求B的转置矩阵
    Debug.Print "and the Bt matrix is:"
    ShowMatrix Bt
    
    Matrix_Multy QBt, Q, Bt
    Debug.Print "the QBt matrix is:"
    ShowMatrix QBt
    
    Matrix_Multy Nbb, B, QBt        '法方程系数矩阵
    Debug.Print "The Nbb matrix is:"
    ShowMatrix Nbb
    
    Debug.Print "The W matrix is:"
    For i = 1 To b1
        Debug.Print W(i)
    Next i
    MajorInColGuass Nbb, W, K       '解法方程
    Debug.Print "The K matrix is:"
    For i = 1 To b1
        Debug.Print K(i)
    Next i
    'Seidel Nbb, W, K, 0.0000001
    Matrix_Multy V, QBt, K          '求改正数
    Debug.Print "The V matrix is:"
    For i = 1 To b2
        Debug.Print V(i)
    Next i
End Sub