jiesuanjuzhen.txt

解算高阶矩阵

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
  
  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

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

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
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
Function lemiv(m As Integer, n As Integer, dbla() As Double, dblb() As Double, dblx() 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
   
  If (Not minv(m, n, dbla(), dblap(), dblu(), dblv(), ka, eps)) Then
     lemiv = False
     Exit Function
  End If
  
  For i = 1 To n
    dblx(i) = 0#
    For j = 1 To m
      dblx(i) = dblx(i) + dblap(i, j) * dblb(j)
    Next j
  Next i
  
  lemiv = True
  
End Function
Sub mmul(m As Integer, n As Integer, l As Integer, mtxa() As Double, mtxb() As Double, mtxc() As Double)
  Dim i As Integer, j As Integer, k As Integer
   For i = 1 To m
     For j = 1 To l
       mtxc(i, j) = 0#
       For k = 1 To n
         mtxc(i, j) = mtxc(i, j) + mtxa(i, k) * mtxb(k, j)
       Next k
     Next j
  Next i
  

End Sub