vb 矩阵计算程序（转载）

测量专用 · 发表于 2020-6-1 10:43

'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：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

测量专用 · 发表于 2020-6-1 10:48

'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：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

      ' 求解失败
      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

测量专用 · 发表于 2020-6-1 20:30

'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：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

测量专用 · 发表于 2020-6-1 20:36

'''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：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

测量专用 · 发表于 2020-6-2 09:13

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：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

测量专用 · 发表于 2020-6-2 09:18

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：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
         End If

         If (kk <= k) Then
            For i = kk To m
                  dblU(i, kk) = dblA(i, kk)
            Next i
         End If

         If (kk <= l) Then
            d = 0#
            For i = kk + 1 To n
                  d = d + e(i) * e(i)
            Next i

            e(kk) = Sqr(d)
            If (e(kk) <> 0#) Then
                  If (e(kk + 1) <> 0#) Then
                     e(kk) = Abs(e(kk))
                     If (e(kk + 1) < 0#) Then
                        e(kk) = -e(kk)
                     End If
                  End If
                  For i = kk + 1 To n
                  e(i) = e(i) / e(kk)
                  Next i
                  e(kk + 1) = 1# + e(kk + 1)
            End If

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

测量专用 · 发表于 2020-6-2 09:18

mm = n
If (m + 1 < n) Then mm = m + 1
If (k < n) Then s(k + 1) = dblA(k + 1, k + 1)
If (m < mm) Then s(mm) = 0#
If (l + 1 < mm) Then e(l + 1) = dblA(l + 1, mm)

e(mm) = 0#
nn = m
If (m > n) Then nn = n
If (nn >= k + 1) Then
      For j = k + 1 To nn
         For i = 1 To m
            dblU(i, j) = 0#
         Next i
         dblU(j, j) = 1#
      Next j
End If

If (k >= 1) Then
      For ll = 1 To k
         kk = k - ll + 1
         If (s(kk) <> 0#) Then
            If (nn >= kk + 1) Then
                  For j = kk + 1 To nn
                     d = 0#
                     For i = kk To m
                        d = d + dblU(i, kk) * dblU(i, j) / dblU(kk, kk)
                     Next i
                     d = -d
                     For i = kk To m
                        dblU(i, j) = dblU(i, j) + d * dblU(i, kk)
                     Next i
                  Next j
            End If

            For i = kk To m
                  dblU(i, kk) = -dblU(i, kk)
            Next i

            dblU(kk, kk) = 1# + dblU(kk, kk)
            If (kk - 1 >= 1) Then
                  For i = 1 To kk - 1
                  dblU(i, kk) = 0#
                  Next i
            End If
         Else
            For i = 1 To m
                  dblU(i, kk) = 0#
            Next i
            dblU(kk, kk) = 1#
         End If
      Next ll
End If

For ll = 1 To n
      kk = n - ll + 1
      If ((kk <= l) And (e(kk) <> 0#)) Then
         For j = kk + 1 To n
            d = 0#
            For i = kk + 1 To n
                  d = d + dblV(i, kk) * dblV(i, j) / dblV(kk + 1, kk)
            Next i
            d = -d
            For i = kk + 1 To n
                  dblV(i, j) = dblV(i, j) + d * dblV(i, kk)
            Next i
         Next j
      End If

      For i = 1 To n
         dblV(i, kk) = 0#
      Next i

      dblV(kk, kk) = 1#
Next ll

For i = 1 To m
      For j = 1 To n
         dblA(i, j) = 0#
      Next j
Next i

m1 = mm
it = 60
While (True)
      If (mm = 0) Then
         Call Cal1(dblA, e, s, dblV, m, n)
         MUav = True
         Exit Function
      End If
      If (it = 0) Then
         Call Cal1(dblA, e, s, dblV, m, n)
         MUav = False
         Exit Function
      End If

      kk = mm - 1
      While ((kk <> 0) And (Abs(e(kk)) <> 0#))
         d = Abs(s(kk)) + Abs(s(kk + 1))
         dd = Abs(e(kk))
         If (dd > eps * d) Then
            kk = kk - 1
         Else
            e(kk) = 0#
         End If
      Wend

      If (kk = mm - 1) Then
         kk = kk + 1
         If (s(kk) < 0#) Then
            s(kk) = -s(kk)
            For i = 1 To n
                  dblV(i, kk) = -dblV(i, kk)
            Next i
         End If

         While ((kk <> m1) And (s(kk) < s(kk + 1)))
            d = s(kk)
            s(kk) = s(kk + 1)
            s(kk + 1) = d
            If (kk < n) Then
                  For i = 1 To n
                     d = dblV(i, kk)
                     dblV(i, kk) = dblV(i, kk + 1)
                     dblV(i, kk + 1) = d
                  Next i
            End If
            If (kk < m) Then
                  For i = 1 To m
                     d = dblU(i, kk)
                     dblU(i, kk) = dblU(i, kk + 1)
                     dblU(i, kk + 1) = d
                  Next i
            End If
            kk = kk + 1
         Wend
         it = 60
         mm = mm - 1
      Else
         ks = mm
         While ((ks > kk) And (Abs(s(ks)) <> 0#))
            d = 0#
            If (ks <> mm) Then d = d + Abs(e(ks))
            If (ks <> kk + 1) Then d = d + Abs(e(ks - 1))
            dd = Abs(s(ks))
            If (dd > eps * d) Then
                  ks = ks - 1
            Else
                  s(ks) = 0#
            End If
         Wend
         If (ks = kk) Then
            kk = kk + 1
            d = Abs(s(mm))
            t = Abs(s(mm - 1))
            If (t > d) Then d = t
            t = Abs(e(mm - 1))
            If (t > d) Then d = t
            t = Abs(s(kk))
            If (t > d) Then d = t
            t = Abs(e(kk))
            If (t > d) Then d = t
            sm = s(mm) / d
            sm1 = s(mm - 1) / d
            em1 = e(mm - 1) / d
            sk = s(kk) / d
            ek = e(kk) / d
            b = ((sm1 + sm) * (sm1 - sm) + em1 * em1) / 2#
            c = sm * em1
            c = c * c
            shh = 0#
            If ((b <> 0#) Or (c <> 0#)) Then
                  shh = Sqr(b * b + c)
                  If (b < 0#) Then shh = -shh
                  shh = c / (b + shh)
            End If
            fg(1) = (sk + sm) * (sk - sm) - shh
            fg(2) = sk * ek
            For i = kk To mm - 1
                  Call Cal2(fg, cs)
                  If (i <> kk) Then e(i - 1) = fg(1)
                  fg(1) = cs(1) * s(i) + cs(2) * e(i)
                  e(i) = cs(1) * e(i) - cs(2) * s(i)
                  fg(2) = cs(2) * s(i + 1)
                  s(i + 1) = cs(1) * s(i + 1)
                  If ((cs(1) <> 1#) Or (cs(2) <> 0#)) Then
                     For j = 1 To n
                        d = cs(1) * dblV(j, i) + cs(2) * dblV(j, i + 1)
                        dblV(j, i + 1) = -cs(2) * dblV(j, i) + cs(1) * dblV(j, i + 1)
                        dblV(j, i) = d
                     Next j
                  End If
                  Call Cal2(fg, cs)
                  s(i) = fg(1)
                  fg(1) = cs(1) * e(i) + cs(2) * s(i + 1)
                  s(i + 1) = -cs(2) * e(i) + cs(1) * s(i + 1)
                  fg(2) = cs(2) * e(i + 1)
                  e(i + 1) = cs(1) * e(i + 1)
                  If (i < m) Then
                     If ((cs(1) <> 1#) Or (cs(2) <> 0#)) Then
                        For j = 1 To m
                              d = cs(1) * dblU(j, i) + cs(2) * dblU(j, i + 1)
                              dblU(j, i + 1) = -cs(2) * dblU(j, i) + cs(1) * dblU(j, i + 1)
                              dblU(j, i) = d
                        Next j
                     End If
                  End If
            Next i
            e(mm - 1) = fg(1)
            it = it - 1
         Else
            If (ks = mm) Then
                  kk = kk + 1
                  fg(2) = e(mm - 1)
                  e(mm - 1) = 0#
                  For ll = kk To mm - 1
                     i = mm + kk - ll - 1
                     fg(1) = s(i)
                     Call Cal2(fg, cs)
                     s(i) = fg(1)
                     If (i <> kk) Then
                        fg(2) = -cs(2) * e(i - 1)
                        e(i - 1) = cs(1) * e(i - 1)
                     End If
                     If ((cs(1) <> 1#) Or (cs(2) <> 0#)) Then
                        For j = 1 To n
                           d = cs(1) * dblV(j, i) + cs(2) * dblV(j, mm)
                           dblV(j, mm) = -cs(2) * dblV(j, i) + cs(1) * dblV(j, mm)
                           dblV(j, i) = d
                        Next j
                     End If
                  Next ll
            Else
                  kk = ks + 1
                  fg(2) = e(kk - 1)
                  e(kk - 1) = 0#
                  For i = kk To mm
                     fg(1) = s(i)
                     Call Cal2(fg, cs)
                     s(i) = fg(1)
                     fg(2) = -cs(2) * e(i)
                     e(i) = cs(1) * e(i)
                     If ((cs(1) <> 1#) Or (cs(2) <> 0#)) Then
                        For j = 1 To m
                           d = cs(1) * dblU(j, i) + cs(2) * dblU(j, kk - 1)
                           dblU(j, kk - 1) = -cs(2) * dblU(j, i) + cs(1) * dblU(j, kk - 1)
                           dblU(j, i) = d
                        Next j
                     End If
                  Next i
            End If
         End If
      End If
Wend

End Function

测量专用 · 发表于 2020-6-2 09:19

注意16楼和17楼是一个函数

测量专用 · 发表于 2020-6-2 14:51

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：MatrixModule.bas
'  函数名：Cal1
'  功能：  内部过程，供MUav函数调用
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
Sub Cal1(dblA() As Double, e() As Double, s() As Double, dblV() As Double, m As Integer, n As Integer)

Dim i As Integer, j As Integer, p As Integer, q As Integer
Dim d As Double

If (m >= n) Then
      i = n
Else
      i = m
End If

For j = 1 To i - 1
      dblA(j, j) = s(j)
      dblA(j, j + 1) = e(j)
Next j

dblA(i, i) = s(i)

If (m < n) Then dblA(i, i + 1) = e(i)

For i = 1 To n - 1
      For j = i + 1 To n
         d = dblV(i, j)
         dblV(i, j) = dblV(j, i)
         dblV(j, i) = d
      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

测量专用 · 发表于 2020-6-2 14:52

''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'  模块名：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

		自动登录	找回密码
密码			免费注册

[分享] vb 矩阵计算程序（转载）