matrixmodule.bas

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

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

        k = k + 1
        If (nIs <> l) Then
            For j = l To n
                d = mtxA(l, j)
                mtxA(l, j) = mtxA(nIs, j)
                mtxA(nIs, j) = d
            Next j
        End If

        If (nJs <> l) Then
            For i = l To m
                d = mtxA(i, nJs)
                mtxA(i, nJs) = mtxA(i, l)
                mtxA(i, l) = d
             Next i
        End If

        For i = l + 1 To n
            d = mtxA(i, l) / mtxA(l, l)
            For j = l + 1 To n
                mtxA(i, j) = mtxA(i, j) - d * mtxA(l, j)
            Next j
        Next i
    Next l

    ' 求解成功
    MRank = k

End Function

'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：MatrixModule.bas
'  函数名：MDetChol
'  功能：  用乔里斯基分解法求对称正定矩阵行列式的值
'  参数：  n       - Integer型变量，对称正定矩阵的阶数。
'          mtxA    - Double型二维数组，体积为n x n，存放对称正定矩阵，返回时，其下三角部分存放分解后的下三角矩阵，其余元素为0。
'          dblDet  - Double型变量，返回对称正定矩阵行列式的值。
'  返回值：Boolean型，成功为True，失败为False。
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function MDetChol(n As Integer, mtxA() As Double, det As Double) As Boolean
    ' 局部变量
    Dim i As Integer, j As Integer, k As Integer
    Dim d As Double

    ' 矩阵校验失败
    If ((mtxA(1, 1) + 1# = 1#) Or (mtxA(1, 1) < 0#)) Then
        MDetChol = False
        Exit Function
    End If

    ' 开始，赋初值
    mtxA(1, 1) = Sqr(mtxA(1, 1))
    d = mtxA(1, 1)

    For i = 2 To n
        mtxA(i, 1) = mtxA(i, 1) / mtxA(1, 1)
    Next i

    ' 循环计算
    For j = 2 To n
        For k = 1 To j - 1
            mtxA(j, j) = mtxA(j, j) - mtxA(j, k) * mtxA(j, k)
        Next k

        If ((mtxA(j, j) + 1# = 1#) Or (mtxA(j, j) < 0#)) Then
            MDetChol = False
            Exit Function
        End If

        mtxA(j, j) = Sqr(mtxA(j, j))
        d = d * mtxA(j, j)

        For i = j + 1 To n
            For k = 1 To j - 1
                mtxA(i, j) = mtxA(i, j) - mtxA(i, k) * mtxA(j, k)
            Next k

            mtxA(i, j) = mtxA(i, j) / mtxA(j, j)
        Next i
    Next j

    ' 计算行列式值
    det = d * d
    
    ' 下三角矩阵
    For i = 1 To n - 1
        For j = i + 1 To n
            mtxA(i, j) = 0#
        Next j
    Next i

    ' 求解成功
    MDetChol = True

End Function

'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：MatrixModule.bas
'  函数名：MLU
'  功能：  矩阵的三角分解
'  参数：  n       - Integer型变量，矩阵的阶数。
'          mtxA    - Double型二维数组，体积为n x n，存放n阶矩阵，返回时存放Q矩阵。
'          mtxL    - Double型二维数组，体积为n x n，返回时存放下三角矩阵L。
'          mtxU    - Double型二维数组，体积为n x n，返回时存放上三角矩阵L。
'  返回值：Boolean型，成功为True，失败为False。
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function MLU(n As Integer, mtxA() As Double, mtxL() As Double, mtxU() As Double) As Boolean
    ' 局部变量
    Dim i As Integer, j As Integer, k As Integer

    For k = 1 To n - 1
        ' 分解失败
        If (Abs(mtxA(k, k)) + 1# = 1#) Then
            MLU = False
            Exit Function
        End If

        For i = k + 1 To n
            mtxA(i, k) = mtxA(i, k) / mtxA(k, k)
        Next i

        For i = k + 1 To n
            For j = k + 1 To n
                mtxA(i, j) = mtxA(i, j) - mtxA(i, k) * mtxA(k, j)
            Next j
        Next i
    Next k

    For i = 1 To n
        For j = 1 To i
          mtxL(i, j) = mtxA(i, j)
          mtxU(i, j) = 0#
        Next j

        mtxL(i, i) = 1#
        mtxU(i, i) = mtxA(i, i)
        For j = i + 1 To n
          mtxL(i, j) = 0#
          mtxU(i, j) = mtxA(i, j)
        Next j
    Next i

    ' 分解成功
    MLU = True

End Function

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：MatrixModule.bas
'  函数名：MMqr
'  功能：  用豪斯荷尔德变换法对矩阵进行QR分解
'  参数：   m    - Integer型变量。矩阵的行数， m>=n
'           n    - Integer型变量。矩阵的列数，n<=m
'          dblA  - Double型二维数组，体积为n x n。存放待分解矩阵；返回时，存放分解式中的R矩阵.
'          dblQ  - Double型二维数组，体积为m x m。返回时，存放分解式中的Q矩阵
'  返回值： Boolean型。False，失败；True, 成功
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function MMqr(m As Integer, n As Integer, dblA() As Double, dblQ() As Double) As Boolean
    Dim i As Integer, j As Integer, k As Integer, nn As Integer, jj As Integer
    Dim u As Double, alpha As Double, w As Double, t As Double

    If (m < n) Then
        MMqr = False
        Exit Function
    End If

    For i = 1 To m
        For j = 1 To m
            dblQ(i, j) = 0#
            If (i = j) Then
                dblQ(i, j) = 1#
            End If
        Next j
    Next i

    nn = n
    If (m = n) Then
        nn = m - 1
    End If

    For k = 1 To nn
        u = 0#
        For i = k To m
            w = Abs(dblA(i, k))
            If (w > u) Then
                u = w
            End If
        Next i

        alpha = 0#
        For i = k To m
            t = dblA(i, k) / u
            alpha = alpha + t * t
        Next i

        If (dblA(k, k) > 0#) Then
            u = -u
        End If

        alpha = u * Sqr(alpha)
        If (Abs(alpha) + 1# = 1#) Then
            MMqr = False
            Exit Function
        End If

        u = Sqr(2# * alpha * (alpha - dblA(k, k)))
        If ((u + 1#) <> 1#) Then
            dblA(k, k) = (dblA(k, k) - alpha) / u
            For i = k + 1 To m
                dblA(i, k) = dblA(i, k) / u
            Next i

            For j = 1 To m
                t = 0#
                For jj = k To m
                  t = t + dblA(jj, k) * dblQ(jj, j)
                Next jj
                For i = k To m
                    dblQ(i, j) = dblQ(i, j) - 2# * t * dblA(i, k)
                Next i
            Next j

            For j = k + 1 To n
                t = 0#
                For jj = k To m
                  t = t + dblA(jj, k) * dblA(jj, j)
                Next jj
                For i = k To m
                    dblA(i, j) = dblA(i, j) - 2# * t * dblA(i, k)
                Next i
            Next j

            dblA(k, k) = alpha
            For i = k + 1 To m
              dblA(i, k) = 0#
            Next i
        End If
    Next k

    For i = 1 To m - 1
        For j = i + 1 To m
            t = dblQ(i, j)
            dblQ(i, j) = dblQ(j, i)
            dblQ(j, i) = t
        Next j
    Next i

    MMqr = True

End Function

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：MatrixModule.bas
'  函数名：MUav
'  功能：  用豪斯荷尔德变换及变形QR算法对矩阵进行奇异值分解
'  参数：   m    - Integer型变量。系数矩阵的行数， m>=n
'           n    - Integer型变量。系数矩阵的列数，n<=m
'          dblA  - Double型二维数组，体积为m x n。存放待分解矩阵；
'                  返回时，其对角线存放矩阵的奇异值(以非递增次序排列)，其余元素为0。
'          dblU  - Double型二维数组，体积为m x m。返回时，存放奇异值分解式中的左奇异向量U。
'          dblV  - Double型二维数组，体积为n x n。返回时，存放奇异值分解式中的右奇异向量VT。
'           ka  - Integer型变量。ka=max(m,n)+1
'          eps  - Double型变量。奇异值分解函数中的控制精度参数。
'  返回值： Boolean型。False，失败无解；True, 成功
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function MUav(m As Integer, n As Integer, dblA() As Double, dblU() As Double, dblV() As Double, ka As Integer, eps As Double) As Boolean
    ' 局部变量
    Dim i As Integer, j As Integer, k As Integer, l As Integer, it As Integer
    Dim ll As Integer, kk As Integer, mm As Integer, nn As Integer, m1 As Integer, ks As Integer
    Dim d As Double, dd As Double, t As Double, sm As Double, sm1 As Double, em1 As Double, sk As Double, ek As Double
    Dim b As Double, c As Double, shh As Double, fg(2) As Double, cs(2) As Double
    ReDim s(ka) As Double, e(ka) As Double, w(ka) As Double

    it = 60
    k = n
    
    If (m - 1 < n) Then
        k = m - 1
    End If

    l = m
    If (n - 2 < m) Then
        l = n - 2
    End If
    If (l < 0) Then
        l = 0
    End If

    ll = k
    
    If (l > k) Then
        ll = l
    End If

    If (ll >= 1) Then
        For kk = 1 To ll
            If (kk <= k) Then
                d = 0#
                For i = kk To m
                    d = d + dblA(i, kk) * dblA(i, kk)
                Next i

                s(kk) = Sqr(d)
                If s(kk) <> 0# Then
                    If (dblA(kk, kk) <> 0#) Then
                        s(kk) = Abs(s(kk))
                        If (dblA(kk, kk) < 0#) Then
                            s(kk) = -s(kk)
                        End If
                    End If
                    For i = kk To m
                        dblA(i, kk) = dblA(i, kk) / s(kk)
                    Next i
                    dblA(kk, kk) = 1# + dblA(kk, kk)
                End If
                s(kk) = -s(kk)
            End If

            If (n >= kk + 1) Then
                For j = kk + 1 To n
                    If ((kk <= k) And (s(kk) <> 0#)) Then
                        d = 0#
                        For i = kk To m
                            d = d + dblA(i, kk) * dblA(i, j)
                        Next i
                        d = -d / dblA(kk, kk)
                        For i = kk To m
                            dblA(i, j) = dblA(i, j) + d * dblA(i, kk)
                        Next i
                    End If
                    e(j) = dblA(kk, j)
                Next j