matrixmodule.bas

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

        Next j
    Next i
End Sub


''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：MatrixModule.bas
'  函数名：Cal2
'  功能：  内部过程，供MUav函数调用
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Sub Cal2(fg() As Double, cs() As Double)

    Dim r As Double, d As Double

    If ((Abs(fg(1)) + Abs(fg(2))) = 0#) Then
        cs(1) = 1#
        cs(2) = 0#
        d = 0#
    Else
        d = Sqr(fg(1) * fg(1) + fg(2) * fg(2))
        If (Abs(fg(1)) > Abs(fg(2))) Then
            d = Abs(d)
            If (fg(1) < 0#) Then d = -d
        End If

        If (Abs(fg(2)) >= Abs(fg(1))) Then
            d = Abs(d)
            If (fg(2) < 0#) Then d = -d
        End If
        cs(1) = fg(1) / d
        cs(2) = fg(2) / d
    End If

    r = 1#
    If (Abs(fg(1)) > Abs(fg(2))) Then
        r = cs(2)
    Else
        If (cs(1) <> 0#) Then
            r = 1# / cs(1)
        End If
    End If

    fg(1) = d
    fg(2) = r

End Sub

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：MatrixModule.bas
'  函数名：MInv
'  功能：  求矩阵的广义逆
'  参数：   m    - Integer型变量。系数矩阵的行数， m>=n
'           n    - Integer型变量。系数矩阵的列数，n<=m
'          dblA  - Double型二维数组，体积为m x n。存放待分解矩阵；
'                  返回时，其对角线存放矩阵的奇异值(以非递增次序排列)，其余元素为0。
'          dblAP - Double型二维数组，体积为n x m。返回时存放矩阵的广义逆。
'          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 MInv(m As Integer, n As Integer, dblA() As Double, dblAP() 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

    If Not MUav(m, n, dblA, dblU, dblV, ka, eps) Then
        MInv = False
        Exit Function
    End If

    j = n
    If (m < n) Then j = m
    j = j - 1
    k = 0
    While (k <= j)
        If (dblA(k + 1, k + 1) = 0#) Then GoTo o_lable
        k = k + 1
    Wend
    
o_lable:
    k = k - 1
    For i = 0 To n - 1
        For j = 0 To m - 1
            dblAP(i + 1, j + 1) = 0#
            For l = 0 To k
                dblAP(i + 1, j + 1) = dblAP(i + 1, j + 1) + dblV(l + 1, i + 1) * dblU(j + 1, l + 1) / dblA(l + 1, l + 1)
            Next l
        Next j
    Next i

    MInv = True
    
End Function

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：MatrixModule.bas
'  函数名：MSymTri
'  功能：  用豪斯荷尔德变换约化对称矩阵为对称三对角矩阵
'  参数：   n    - Integer型变量，对称矩阵的阶数。
'          dblA  - Double型二维数组，体积为n x n。存放n阶对称矩阵。
'          dblQ  - Double型二维数组，体积为n x n。返回时存放豪斯荷尔德变换的乘积矩阵Q。
'          dblT - Double型二维数组，体积为n x n。返回时存放对称三对角矩阵T。
'          dblB - Double型一维数组，长度为n。返回时存放对称三对角矩阵T主对角线上的元素。
'          dblC - Double型一维数组，长度为n。返回时前n-1个元素存放对称三对角矩阵T次对角线上的元素。
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Sub MSymTri(n As Integer, dblA() As Double, dblQ() As Double, dblT() As Double, dblB() As Double, dblC() As Double)
    Dim i As Integer, j As Integer, k As Integer
    Dim h As Double, f As Double, g As Double, h2 As Double

    For i = 1 To n
        For j = 1 To n
            dblQ(i, j) = dblA(i, j)
        Next j
    Next i

    For i = n To 2 Step -1
        h = 0#
        If (i > 2) Then
            For k = 1 To i - 1
                h = h + dblQ(i, k) * dblQ(i, k)
            Next k
        End If

        If (h + 1# = 1#) Then
            dblC(i) = 0#
            If (i = 2) Then dblC(i) = dblQ(i, i - 1)
            dblB(i) = 0#
        Else
            dblC(i) = Sqr(h)
            
            If (dblQ(i, i - 1) > 0#) Then dblC(i) = -dblC(i)
            
            h = h - dblQ(i, i - 1) * dblC(i)
            dblQ(i, i - 1) = dblQ(i, i - 1) - dblC(i)
            f = 0#

            For j = 1 To i - 1
                dblQ(j, i) = dblQ(i, j) / h
                g = 0#
                For k = 1 To j
                  g = g + dblQ(j, k) * dblQ(i, k)
                Next k

                If (j + 1 <= i - 1) Then
                    For k = j + 1 To i - 1
                        g = g + dblQ(k, j) * dblQ(i, k)
                    Next k
                End If

                dblC(j) = g / h
                f = f + g * dblQ(j, i)
            Next j

            h2 = f / (h + h)
            For j = 1 To i - 1
                f = dblQ(i, j)
                g = dblC(j) - h2 * f
                dblC(j) = g
                For k = 1 To j
                    dblQ(j, k) = dblQ(j, k) - f * dblC(k) - g * dblQ(i, k)
                Next k
            Next j

            dblB(i) = h
        End If
    Next i

    For i = 1 To n - 1
        dblC(i) = dblC(i + 1)
    Next i

    dblC(n) = 0#
    dblB(1) = 0#
    For i = 1 To n
        If ((dblB(i) <> 0#) And (i - 1 >= 0)) Then
          For j = 1 To i - 1
                g = 0#
                For k = 1 To i - 1
                    g = g + dblQ(i, k) * dblQ(k, j)
                Next k
                For k = 1 To i - 1
                    dblQ(k, j) = dblQ(k, j) - g * dblQ(k, i)
                Next k
          Next j
        End If

        dblB(i) = dblQ(i, i)
        dblQ(i, i) = 1#
        If (i - 1 >= 0) Then
            For j = 1 To i - 1
                dblQ(i, j) = 0#
                dblQ(j, i) = 0#
            Next j
        End If
    Next i
    
    ' 构造对称三对角矩阵
    For i = 1 To n
        For j = 1 To n
            dblT(i, j) = 0
            k = i - j
            If k = 0 Then dblT(i, j) = dblB(j)
            If k = 1 Then dblT(i, j) = dblC(j)
            If k = -1 Then dblT(i, j) = dblC(i)
        Next j
    Next i
    
End Sub

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：MatrixModule.bas
'  函数名：MSymTriEigenv
'  功能：  用变形QR方法计算对称三对角矩阵的全部特征值和特征向量
'  参数：   n    - Integer型变量，对称三对角矩阵的阶数。
'          dblB - Double型一维数组，长度为n，存放对称三对角矩阵T主对角线上的元素。返回时，存放全部特征值。
'          dblC - Double型一维数组，长度为n，前n-1个元素存放对称三对角矩阵T次对角线上的元素。
'          dblQ - Double型二维数组，体积为n x n。
'                 1）如果存放单位矩阵，则返回n阶对称三对角矩阵的特征向量组。
'                 2）如果存放对称矩阵A的豪斯荷尔德变换的乘积矩阵Q（可由函数MSymTri求得），则返回n阶对称矩阵A的特征向量组。
'                    其中dblQ中的第j列为与数组dblB中第j个特征值对应的特征向量。
'          eps  - Double型变量。迭代过程中的控制精度参数。
'          nMaxItNum    - Integer。为求得一个特征值所允许的最大迭代次数。
'  返回值： Boolean型。False，失败无解；True, 成功
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Function MSymTriEigenv(n As Integer, dblB() As Double, dblC() As Double, dblQ() As Double, eps As Double, nMaxItNum As Integer) As Boolean
    ' 局部变量
    Dim i As Integer, j As Integer, k As Integer, m As Integer, it As Integer
    Dim d As Double, f As Double, h As Double, g As Double, p As Double, r As Double, e As Double, s As Double

    dblC(n) = 0#
    d = 0#
    f = 0#
    
    For j = 1 To n
        it = 0
        h = eps * (Abs(dblB(j)) + Abs(dblC(j)))
        If (h > d) Then d = h
        m = j
        While ((m <= n) And (Abs(dblC(m)) > d))
            m = m + 1
        Wend

        If (m <> j) Then
            Do
                If (it = nMaxItNum) Then
                    MSymTriEigenv = False
                    Exit Function
                End If

                it = it + 1
                g = dblB(j)
                p = (dblB(j + 1) - g) / (2# * dblC(j))
                r = Sqr(p * p + 1#)

                If (p >= 0#) Then
                    dblB(j) = dblC(j) / (p + r)
                Else
                    dblB(j) = dblC(j) / (p - r)
                End If

                h = g - dblB(j)
                For i = j + 1 To n
                  dblB(i) = dblB(i) - h
                Next i

                f = f + h
                p = dblB(m)
                e = 1#
                s = 0#
                
                For i = m - 1 To j Step -1
                    g = e * dblC(i)
                    h = e * p
                    If (Abs(p) >= Abs(dblC(i))) Then
                        e = dblC(i) / p
                        r = Sqr(e * e + 1#)
                        dblC(i + 1) = s * p * r
                        s = e / r
                        e = 1# / r
                    Else
                        e = p / dblC(i)
                        r = Sqr(e * e + 1#)
                        dblC(i + 1) = s * dblC(i) * r
                        s = 1# / r
                        e = e / r
                    End If

                    p = e * dblB(i) - s * g
                    dblB(i + 1) = h + s * (e * g + s * dblB(i))
                    For k = 1 To n
                        h = dblQ(k, i + 1)
                        dblQ(k, i + 1) = s * dblQ(k, i) + e * h
                        dblQ(k, i) = e * dblQ(k, i) - s * h
                    Next k
                Next i

                dblC(j) = s * p
                dblB(j) = e * p

            Loop While (Abs(dblC(j)) > d)

        End If

        dblB(j) = dblB(j) + f

    Next j

    For i = 1 To n
        k = i
        p = dblB(i)
        If (i + 1 <= n) Then
            j = i + 1
            While ((j <= n) And (dblB(j) <= p))
                k = j
                p = dblB(j)
                j = j + 1
            Wend
        End If

        If (k <> i) Then
            dblB(k) = dblB(i)
            dblB(i) = p
            
            For j = 1 To n
                p = dblQ(j, i)
                dblQ(j, i) = dblQ(j, k)
                dblQ(j, k) = p
            Next j
        End If
    Next i

    MSymTriEigenv = True
    
End Function

'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''