mmath.bas

用于矩阵加减乘除运算、高斯消元法解方程等

Attribute VB_Name = "mMath"
Option Explicit

Private Declare Sub CopyMemory Lib "kernel32" Alias "RtlMoveMemory" (Destination As Any, Source As Any, ByVal Length As Long)


Public Function MatrixTo1DimArray(ByRef SrcMat() As Double, ByRef DArr() As Double)
''矩阵必是2维的，转成数组是1维的
  Dim i As Integer
  If UBound(SrcMat, 1) - LBound(SrcMat, 1) + 1 = 1 Then
     ReDim DArr(LBound(SrcMat, 2) To UBound(SrcMat, 2)) As Double
     For i = LBound(SrcMat, 2) To UBound(SrcMat, 2)
        DArr(i) = SrcMat(LBound(SrcMat, 1), i)
     Next i
     MatrixTo1DimArray = True
  ElseIf UBound(SrcMat, 2) - LBound(SrcMat, 2) + 1 = 1 Then
     ReDim DArr(LBound(SrcMat, 1) To UBound(SrcMat, 1)) As Double
     For i = LBound(SrcMat, 1) To UBound(SrcMat, 1)
        DArr(i) = SrcMat(i, LBound(SrcMat, 1))
     Next i
     MatrixTo1DimArray = True
  Else
     MatrixTo1DimArray = False
  End If
End Function

Public Function CreateMatrixFrom1DimArray(ByRef SrcArr() As Double, ByRef MRsl() As Double, ByVal bAsRowMatrix As Boolean) As Boolean
  Dim al As Integer
  Dim au As Integer
  Dim i As Integer
  al = LBound(SrcArr)
  au = UBound(SrcArr)
  If bAsRowMatrix Then
     ReDim MRsl(0 To 0, al To au) As Double
     For i = al To au
       MRsl(0, i) = SrcArr(i)
     Next i
  Else
     ReDim MRsl(al To au, 0 To 0) As Double
     For i = al To au
       MRsl(i, 0) = SrcArr(i)
     Next i
  End If
End Function

Public Function CreateMatrixFromMatrix(ByRef MSrc() As Double, ByVal iRows As Integer, ByVal iCols As Integer, ByVal bAllowAutoFill0 As Boolean, ByRef MR() As Double) As Boolean
  '根据另一矩阵构造新的矩阵(从第一元素始取)
  'iRows 构成几行
  'iCols 构成几列
  '
  Dim al1 As Integer, al2 As Integer, au1 As Integer, au2 As Integer
  Dim bl1 As Integer, bl2 As Integer, bu1 As Integer, bu2 As Integer
  Dim i As Integer
  Dim j As Integer
  '
  al1 = LBound(MSrc, 1)
  au1 = UBound(MSrc, 1)
  al2 = LBound(MSrc, 2)
  au2 = UBound(MSrc, 2)
  If iRows > au1 - al1 + 1 Or iCols > au2 - al2 + 1 Then
     If Not bAllowAutoFill0 Then     '允许在不够取时自动填0
        CreateMatrixFromMatrix = False
        Exit Function
     End If
  End If
  
   '
     ReDim MR(al1 To iRows - al1 - 1, al2 To iCols - al2 - 1) As Double
     For i = al1 To iRows - al1 - 1
       For j = al2 To iCols - al2 - 1
          If j > au2 Or i > au1 Then
             MR(i, j) = 0#
          Else
             MR(i, j) = MSrc(i, j)
          End If
       Next j
     Next i
     CreateMatrixFromMatrix = True
  
End Function

Public Sub CreateSlopeMatrix(ByRef aSlope() As Double, ByRef MR() As Double)
   '根据一维数组aSlope()值构建对角线矩阵，如平差中的权系数阵
   '             a  0  0  0
   '             0  b  0  0
   '             0  0  c  0
   '             0  0  0  d
   Dim al As Integer
   Dim au As Integer
   Dim i As Integer
   Dim j As Integer
   al = LBound(aSlope)
   au = UBound(aSlope)
   ReDim MR(al To au, al To au) As Double
   For i = al To au
     For j = al To au
        If i = j Then
           MR(i, j) = aSlope(j)
        Else
           MR(i, j) = 0#
        End If
     Next j
  Next i
End Sub


Public Sub CreateUnitMatrix(ByVal iRows As Integer, ByRef MR() As Double)
  Dim i As Integer
  Dim j As Integer
  ReDim MR(0 To iRows - 1, 0 To iRows - 1) As Double
  For i = 0 To iRows - 1
     For j = 0 To iRows - 1
       If i = j Then
          MR(i, j) = 1#
       Else
          MR(i, j) = 0#
       End If
     Next j
  Next i
End Sub

Public Function CreateMatrixByPara(ByRef RslMatrix() As Double, ByVal iNeedRows As Integer, ByVal iNeedCols As Integer, ByVal bAllowMoreAtom As Boolean, ParamArray pSA()) As Boolean
  '根据参数表构造矩阵
  'iNeedRows 构成几行
  'iNeedCols 构成几列
  '
   Dim pNum As Integer
   Dim colNum As Integer
   Dim pCount As Integer
   Dim i As Integer
   Dim j As Integer
   '
   pNum = UBound(pSA) - LBound(pSA) + 1
   colNum = Sqr(pNum)
   If iNeedRows = 0 Then '自动根据列数测行数
      If iNeedCols = 0 Then
         colNum = Fix(Sqr(pNum))
         iNeedRows = colNum
         iNeedCols = colNum
      Else
         iNeedRows = Fix(pNum / iNeedCols)
      End If
   Else
      If iNeedCols = 0 Then
         iNeedCols = Fix(pNum / iNeedRows)
      End If
   End If
   If iNeedRows = 0 Or iNeedCols = 0 Then
      CreateMatrixByPara = False
      Exit Function
   End
     If iNeedRows * iNeedCols > pNum Then
        CreateMatrixByPara = False
        Exit Function
     End If
     If iNeedRows * iNeedCols < pNum And (Not bAllowMoreAtom) Then
        CreateMatrixByPara = False
        Exit Function
     End If
   End If
   '
     ReDim RslMatrix(0 To iNeedRows - 1, 0 To iNeedCols - 1) As Double
     pCount = LBound(pSA)
     For i = 0 To iNeedRows - 1
       For j = 0 To iNeedCols - 1
         RslMatrix(i, j) = CDbl(pSA(pCount))
         pCount = pCount + 1
       Next j
     Next i
     ShowMatrix RslMatrix
     CreateMatrixByPara = True
End Function



Public Sub MatrixTransPos(ByRef MA() As Double, ByRef MR() As Double)
  '矩阵转置，MA转置后返回MR
  Dim al1 As Integer, al2 As Integer, au1 As Integer, au2 As Integer
  Dim bl1 As Integer, bl2 As Integer, bu1 As Integer, bu2 As Integer
  Dim i As Integer
  Dim j As Integer
  '
  al1 = LBound(MA, 1)
  au1 = UBound(MA, 1)
  al2 = LBound(MA, 2)
  au2 = UBound(MA, 2)
  ReDim MR(al2 To au2, al1 To au1) As Double
  For i = al1 To au1
       For j = al2 To au2
          MR(j, i) = MA(i, j)
       Next j
  Next i
    
End Sub

Public Function MatrixMulti(ByRef MA() As Double, ByRef MB() As Double, ByRef MR() As Double) As Boolean
   'Mr=MA*MB
   'MA、MB可以是一维数组
Dim al1 As Integer, al2 As Integer, au1 As Integer, au2 As Integer
Dim bl1 As Integer, bl2 As Integer, bu1 As Integer, bu2 As Integer
Dim i As Integer
Dim j As Integer
Dim k As Integer
'
al1 = LBound(MA, 1)
au1 = UBound(MA, 1)
'
al2 = LBound(MA, 2)
au2 = UBound(MA, 2)
   
bl1 = LBound(MB, 1)
bu1 = UBound(MB, 1)
bl2 = LBound(MB, 2)
bu2 = UBound(MB, 2)
'
If au2 - al2 <> bu1 - bl1 Then
   MatrixMulti = False
   Exit Function
Else
   ReDim MR(al1 To au1, bl2 To bu2) As Double
   For i = al1 To au1
       For j = bl2 To bu2
           MR(i, j) = 0#
       Next j
   Next i
   For i = al1 To au1
       For j = bl2 To bu2
           For k = al2 To au2
               MR(i, j) = MR(i, j) + MA(i, k) * MB(k - al2 + bl2, j)
           Next k
       Next j
   Next i
   '
   MatrixMulti = True
End If
End Function

Public Function EQSolution(ByRef a() As Double, ByRef rsl() As Double, ByVal KeepCoefficientMatrix As Boolean) As Integer
'列主消元法解方程组
'A() 为方程源系数阵,调用前定义为 Dim A(MMM,MMM+1) as double 格式,MMM为方程个数
'    如果A()的第二维数超过一维+1,则多余第二维值忽略
'a0 X + a1 Y + a2 Z + a3 = 0
'b0 X + b1 Y + b2 Z + b3 = 0
'c0 X + c1 Y + c2 Z + c3 = 0 euquation
'A(*)={ {a0,a1,a2,a3}
'       {b0,b1,b2,b3}
'       {c0,c1,c2,c3}
'     }
'返回Rsl()一维数组为方程组的解, 下标从lbound(rsl) to Ubound(Rsl)
'When KeepCoefficientMatrix =false, A() will be changed after EQSolution
'Returns (Not 0) if errors or return 0 if OK!
'0----OK!
'1----方程组维数不对 Rows+1 > Cols
'3----方程无解
Dim i As Integer, j As Integer, k As Integer
Dim ExChg As Double
Dim Dims1 As Integer, Dims2 As Integer
Dim l1 As Integer, l2 As Integer
Dim u1 As Integer, u2 As Integer
Dim MaxElement As Double, MainElement As Double
Dim MaxNo As Integer
l1 = LBound(a, 1)
u1 = UBound(a, 1)
l2 = LBound(a, 2)
u2 = UBound(a, 2)

Dim ss As String
Dim ii As Integer, jj As Integer

Dims1 = u1 - l1 + 1: Dims2 = u2 - l2 + 1
If Dims1 + 1 > Dims2 Then
  EQSolution = 1  '方程组无效
  Exit Function
End If
'
ReDim rsl(l1 To u1) As Double
u2 = u1 - l1 + l2 + 1
If KeepCoefficientMatrix Then
  Dim BackA() As Double
  ReDim BackA(l1 To u1, l2 To u2) As Double
  Call CopyMemory(ByVal VarPtr(BackA(l1, l2)), ByVal VarPtr(a(l1, l2)), Dims1 * Dims2 * Len(BackA(l1, l2)))
End If
'方程消元变换
For i = l1 To u1
  '求最大主元
  'MsgBox "消元 i=" & i
  'GoSub eqSub
  MaxElement = a(i, i - l1 + l2)
  MaxNo = i
  For j = i To u1
    If Abs(a(j, i - l1 + l2)) > Abs(MaxElement) Then MaxElement = Abs(a(j, i - l1 + l2)): MaxNo = j
  Next j
  If MaxNo <> i Then
    For j = l2 To u2
      ExChg = a(MaxNo, j): a(MaxNo, j) = a(i, j): a(i, j) = ExChg
    Next j
  End If
  If a(i, i - l1 + l2) = 0 Then
    EQSolution = 3          '主元系数为0,方程无解
    'MsgBox "奇异方程！无解！", vbCritical, "解方程"
    GoTo SUEND
  End If
  '消主元
  MainElement = a(i, i - l1 + l2)
  For j = i - l1 + l2 To u2
    a(i, j) = a(i, j) / MainElement
  Next j
  For j = i + 1 To u1
    If a(j, i - l1 + l2) <> 0 Then
      MainElement = a(j, i - l1 + l2)
      For k = i - l1 + l2 To u2
        a(j, k) = a(j, k) / MainElement - a(i, k)
      Next k
    End If
  Next j
Next i

  'MsgBox "消元 后="
  'GoSub eqSub


'求解
rsl(u1) = 0 - a(u1, u2)
For i = u1 - 1 To l1 Step -1
  rsl(i) = 0 - a(i, u2)
  For j = i + l2 - l1 + 1 To u2 - 1
    rsl(i) = rsl(i) - a(i, j) * rsl(j + l1 - l2)
  Next j
Next i

'  MsgBox "解算后="
'  GoSub RslSub

EQSolution = 0
SUEND:
If KeepCoefficientMatrix Then
  Call CopyMemory(ByVal VarPtr(a(l1, l2)), ByVal VarPtr(BackA(l1, l2)), Dims1 * Dims2 * Len(BackA(l1, l2)))
End If
Exit Function

eqSub:
ss = ""
For ii = l1 To u1
   ss = ss & vbCrLf
   For jj = l2 To u2
      ss = ss & Format(Format(a(ii, jj), "0.00000"), "@@@@@@@@@@@@@") & ","
   Next jj
Next ii
MsgBox ss
Return

RslSub:
Return
ss = ""
For i = LBound(rsl) To UBound(rsl)
   ss = ss & vbCrLf & rsl(i)
Next i
MsgBox ss
Return
End Function

Public Sub ShowArray(ByRef aa() As Double, Optional ByVal sTitle As String = "")
  Dim i As Integer
  Dim ss As String
   For i = LBound(aa) To UBound(aa)
     ss = ss & vbCrLf & " " & StrExtendToLength(i, 5, " ", AlignmentLeft, True, "0") & " :     " & aa(i)
   Next i
   If sTitle = "" Then
      MsgBox ss, vbOKOnly, "数组显示"
   Else
      MsgBox ss, vbOKOnly, "数组显示 ： " & sTitle
   End If
   
End Sub


 Public Sub ShowMatrix(ByRef MA() As Double, Optional ByVal sTitle As String = "")
  Dim i As Integer
  Dim j As Integer
  Dim ss As String
  Dim spLen As Integer
  Dim iTestDim2 As Integer
  
  On Error Resume Next
  Err.Clear

  iTestDim2 = UBound(MA, 2)
  If Err Then
     On Error GoTo 0
     ShowArray MA, sTitle & "(一维数组)"
  Else
     spLen = 12
     ss = "行/列    "
     For j = LBound(MA, 2) To UBound(MA, 2)
        ss = ss & StrExtendToLength(j, spLen, " ", AlignmentLeft, True, "0") & "    "
     Next j
     ss = ss & vbCrLf
     For j = LBound(MA, 2) To UBound(MA, 2)
        ss = ss & String(spLen + 5, "-")
     Next j
      
     For i = LBound(MA, 1) To UBound(MA, 1)
       ss = ss & vbCrLf & " " & StrExtendToLength(i, 5, " ", AlignmentLeft, True, "0") & "  |  "
       For j = LBound(MA, 2) To UBound(MA, 2)
         ss = ss & StrExtendToLength(MA(i, j), spLen, " ", AlignmentLeft, True, "0.000") & "    "
       Next j
     Next i
     If sTitle = "" Then
        MsgBox ss, vbOKOnly, "矩阵显示"
     Else
        MsgBox ss, vbOKOnly, "矩阵显示 ： " & sTitle
     End If
  End If
  On Error GoTo 0
End Sub


Public Function ParaMatrixInvert(ByRef rsl() As Double, ParamArray pSA()) As Integer
   Dim pNum As Integer
   Dim colNum As Integer
   Dim pCount As Integer
   Dim i As Integer
   Dim j As Integer
   Dim bl As Boolean
   Dim MIrsl As Integer
   Dim SA() As Double
   '
   pNum = UBound(pSA) - LBound(pSA) + 1
   colNum = Sqr(pNum)
   If colNum * colNum <> pNum Then
      ParaMatrixInvert = False
      Exit Function
   Else
     ReDim SA(0 To colNum - 1, 0 To colNum - 1)
     pCount = LBound(pSA)
     For i = 0 To colNum - 1
       For j = 0 To colNum - 1
         SA(i, j) = CDbl(pSA(pCount))
         pCount = pCount + 1
       Next j
     Next i
     ShowMatrix SA
     MIrsl = MatrixInvert(SA, rsl)
     bl = (MIrsl = 0)
     If bl Then
        ShowMatrix rsl
     Else
        MsgBox "求逆失败！"
     End If
     ParaMatrixInvert = bl
   End If
End Function


'----------------------------------------------------------------------
Public Function MatrixInvert(ByRef SA() As Double, ByRef rsl() As Double) As Integer
'矩阵求逆
'SA() 源矩阵, Dim A(MMM,MMM) as double
'a0  a1  a2
'b0  b1  b2
'c0  c1  c2
'SA(*)={ {a0,a1,a2}
'        {b0,b1,b2}
'        {c0,c1,c2}
'      }
'返回Rsl()与A()同维逆阵
'Returns (Not 0) if errors or return 0 if OK!
'0----OK!
'1----维数不对
'3----不可逆
Dim Dims1, Dims2 As Integer
Dim l1, l2, u1, u2 As Integer
Dim rl1, ru1, rl2, ru2 As Integer
Dim MaxElement, MainElement As Double
Dim MaxNo As Integer
Dim i As Integer, j As Integer, k As Integer
Dim ExChg As Double
'
l1 = LBound(SA, 1)
u1 = UBound(SA, 1)
l2 = LBound(SA, 2)
u2 = UBound(SA, 2)
ReDim rsl(l1 To u1, l2 To u2) As Double
rl1 = LBound(rsl, 1)
ru1 = UBound(rsl, 1)
rl2 = LBound(rsl, 2)
ru2 = UBound(rsl, 2)