matrixmodule.bas

清华大学2002年出版的《科学与工程数值计算算法VB》配套源码

        Next i

        ' 求解失败
        If (d + 1# = 1#) Then
            MCinv = False
            Exit Function
        End If

        If (nIs(k) <> k) Then
          For j = 1 To n
              t = mtxAR(k, j)
              mtxAR(k, j) = mtxAR(nIs(k), j)
              mtxAR(nIs(k), j) = t

              t = mtxAI(k, j)
              mtxAI(k, j) = mtxAI(nIs(k), j)
              mtxAI(nIs(k), j) = t
            Next j
        End If
        If (nJs(k) <> k) Then
            For i = 1 To n
                t = mtxAR(i, k)
                mtxAR(i, k) = mtxAR(i, nJs(k))
                mtxAR(i, nJs(k)) = t
                t = mtxAI(i, k)
                mtxAI(i, k) = mtxAI(i, nJs(k))
                mtxAI(i, nJs(k)) = t
            Next i
        End If

        mtxAR(k, k) = mtxAR(k, k) / d
        mtxAI(k, k) = -mtxAI(k, k) / d
        For j = 1 To n
            If (j <> k) Then
                p = mtxAR(k, j) * mtxAR(k, k)
                q = mtxAI(k, j) * mtxAI(k, k)
                s = (mtxAR(k, j) + mtxAI(k, j)) * (mtxAR(k, k) + mtxAI(k, k))
                mtxAR(k, j) = p - q
                mtxAI(k, j) = s - p - q
            End If
        Next j
        For i = 1 To n
            If (i <> k) Then
              For j = 1 To n
                    If (j <> k) Then
                        p = mtxAR(k, j) * mtxAR(i, k)
                        q = mtxAI(k, j) * mtxAI(i, k)
                        s = (mtxAR(k, j) + mtxAI(k, j)) * (mtxAR(i, k) + mtxAI(i, k))
                        t = p - q
                        b = s - p - q
                        mtxAR(i, j) = mtxAR(i, j) - t
                        mtxAI(i, j) = mtxAI(i, j) - b
                    End If
                Next j
            End If
        Next i

        For i = 1 To n
            If (i <> k) Then
                p = mtxAR(i, k) * mtxAR(k, k)
                q = mtxAI(i, k) * mtxAI(k, k)
                s = (mtxAR(i, k) + mtxAI(i, k)) * (mtxAR(k, k) + mtxAI(k, k))
                mtxAR(i, k) = q - p
                mtxAI(i, k) = p + q - s
            End If
        Next i
    Next k

    ' 调整恢复行列次序
    For k = n To 1 Step -1
        If (nJs(k) <> k) Then
          For j = 1 To n
              t = mtxAR(k, j)
              mtxAR(k, j) = mtxAR(nJs(k), j)
              mtxAR(nJs(k), j) = t
              t = mtxAI(k, j)
              mtxAI(k, j) = mtxAI(nJs(k), j)
              mtxAI(nJs(k), j) = t
            Next j
        End If
        If (nIs(k) <> k) Then
            For i = 1 To n
                t = mtxAR(i, k)
                mtxAR(i, k) = mtxAR(i, nIs(k))
                mtxAR(i, nIs(k)) = t
                t = mtxAI(i, k)
                mtxAI(i, k) = mtxAI(i, nIs(k))
                mtxAI(i, nIs(k)) = t
            Next i
        End If
    Next k
    
    ' 求解成功
    MCinv = True

End Function

'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：MatrixModule.bas
'  函数名：MSsgj
'  功能：  计算对称正定矩阵的逆
'  参数：  n     - Integer型变量，矩阵阶数。
'          mtxA  - Double型二维数组，体积为n x n。存放实对称正定矩阵A；返回时存放逆矩阵A-。
'  返回值：Boolean型，成功为True，失败为False。
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function MSsgj(n As Integer, mtxA() As Double) As Boolean
    ' 局部变量
    Dim i As Integer, j, k, m
    Dim w As Double, g As Double
    ReDim b(n) As Double

    ' 循环重新编号求解
    For k = 1 To n
        w = a(1, 1)
         ' 求解失败
        If (Abs(w) + 1# = 1#) Then
            MSsgj = False
            Exit Function
        End If

        m = n - k - 1
        For i = 2 To n
            g = a(i, 1)
            b(i) = g / w
            If (i <= m) Then b(i) = -b(i)
            For j = 2 To i
              a(i - 1, j - 1) = a(i, j) + g * b(j)
            Next j
        Next i
        a(n, n) = 1# / w
        For i = 2 To n
          a(n, i - 1) = b(i)
        Next i
    Next k

    For i = 1 To n - 1
        For j = i + 1 To n
            a(i, j) = a(j, i)
        Next j
    Next i

' 求解成功
MSsgj = True

End Function

'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：MatrixModule.bas
'  函数名：MTinv
'  功能：  用特兰持（Trench）方法求托伯利兹（Toeplitz）矩阵的逆矩阵
'  参数：  n      - Integer型变量，T型矩阵阶数。
'          dblT   - Double型一维数组，长度为n。存放n阶T型矩阵中的上三角元素t0,t1,…tn-1。
'          dblTT  - Double型一维数组，长度为n。其中后n-1个元素tt(1),…,tt(n-1)依次存放n阶T型矩阵中的元素。
'          dblB    - Double型二维数组，体积为n x n。返回n阶T型矩阵的逆矩阵。
'  返回值：Boolean型，成功为True，失败为False。
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function MTinv(n As Integer, dblT() As Double, dblTT() As Double, dblB() As Double) As Boolean
    ' 局部变量
    Dim i As Integer, j As Integer, k As Integer
    Dim a As Double, s As Double
    ReDim c(n) As Double, r(n) As Double, p(n) As Double

    ' 矩阵非T型矩阵
    If (Abs(dblT(1)) + 1# = 1#) Then
        MTinv = False
        Exit Function
    End If

    ' 取初值
    a = dblT(1)
    c(1) = dblTT(2) / dblT(1)
    r(1) = dblT(2) / dblT(1)

    ' 循环计算
    For k = 1 To n - 2
        s = 0#
        For j = 2 To k + 1
          s = s + c(k + 1 - j + 1) * dblTT(j)
        Next j

        s = (s - dblTT(k + 2)) / a

        For i = 1 To k
          p(i) = c(i) + s * r(k - i + 1)
        Next i

        c(k + 1) = -s
        s = 0#
        For j = 2 To k + 1
          s = s + r(k + 1 - j + 1) * dblT(j)
        Next j

        s = (s - dblT(k + 2)) / a
        For i = 1 To k
            r(i) = r(i) + s * c(k - i + 1)
            c(k - i + 1) = p(k - i + 1)
        Next i
        
        r(k + 1) = -s
        a = 0#
        For j = 2 To k + 2
          a = a + dblT(j) * c(j - 1)
        Next j

        a = dblT(1) - a
        If (Abs(a) + 1# = 1#) Then
            MTinv = False
            Exit Function
        End If
    Next k

    dblB(1, 1) = 1# / a
    For i = 1 To n - 1
        dblB(1, i + 1) = -r(i) / a
        dblB(i + 1, 1) = -c(i) / a
    Next i

    ' 计算逆矩阵中的各元素
    For i = 1 To n - 1
        For j = 1 To n - 1
            dblB(i + 1, j + 1) = dblB(i, j) - c(i) * dblB(1, j + 1)
            dblB(i + 1, j + 1) = dblB(i + 1, j + 1) + c(n - j) * dblB(1, n - i + 1)
        Next j
    Next i

    ' 求解成功
    MTinv = True
    
End Function

'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：MatrixModule.bas
'  函数名：MDetGauss
'  功能：  用全选主元高斯消去法求行列式的值
'  参数：  n       - Integer型变量，方阵的阶数。
'          mtxA    - Double型二维数组，体积为n x n，存放方阵。
'  返回值：Double型，行列式的值。
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function MDetGauss(n As Integer, mtxA() As Double) As Double
    ' 局部变量
    Dim i As Integer, j As Integer, k As Integer, nIs As Integer, nJs As Integer
    Dim f As Double, det As Double, q As Double, d As Double

    f = 1#
    det = 1#

    ' 选主元
    For k = 1 To n - 1
        q = 0#
        For i = k To n
            For j = k To n
                d = Abs(mtxA(i, j))
                If (d > q) Then
                    q = d
                    nIs = i
                    nJs = j
                End If
            Next j
        Next i

        ' 求解失败
        If (q + 1# = 1#) Then
            MDetGauss = 0
            Exit Function
        End If

        If (nIs <> k) Then
            f = -f
            For j = k To n
                d = mtxA(k, j)
                mtxA(k, j) = mtxA(nIs, j)
                mtxA(nIs, j) = d
            Next j
        End If

        ' 调整
        If (nJs <> k) Then
            f = -f
            For i = k To n
                d = mtxA(i, nJs)
                mtxA(i, nJs) = mtxA(i, k)
                mtxA(i, k) = d
            Next i
        End If

        ' 计算行列式的值
        det = det * mtxA(k, k)
        
        ' 调整方阵为上三角矩阵
        For i = k + 1 To n
            d = mtxA(i, k) / mtxA(k, k)
            For j = k + 1 To n
                mtxA(i, j) = mtxA(i, j) - d * mtxA(k, j)
            Next j
        Next i
    Next k

    ' 计算行列式的值
    det = f * det * mtxA(n, n)

    ' 求解成功
    MDetGauss = det

End Function

'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：MatrixModule.bas
'  函数名：MRank
'  功能：  用全选主元高斯消去法求矩阵的秩
'  参数：  m       - Integer型变量，矩阵的行数。
'          n       - Integer型变量，矩阵的列数。
'          mtxA    - Double型二维数组，体积为m x n，存放待求秩的矩阵。
'  返回值：Integer型，矩阵的秩。
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function MRank(m As Integer, n As Integer, mtxA() As Double) As Integer
    ' 局部变量
    Dim i As Integer, j As Integer, k As Integer, l As Integer, nIs As Integer, nJs As Integer, nn As Integer
    Dim q As Double, d As Double

    ' 基数
    nn = m
    If (m >= n) Then nn = n

    ' 消元求解
    k = 0
    For l = 1 To nn
        q = 0#
        For i = 2 To m
            For j = 2 To n
                d = Abs(mtxA(i, j))
                If (d > q) Then
                    q = d
                    nIs = i
                    nJs = j
                End If
            Next j
        Next i