Im Beispiel wird die Inverse der Hilbert-Matrix 6. Ordnung berechnet:

1/1, 1/2, 1/3, 1/4, 1/5, 1/6 
1/2, 1/3, 1/4, 1/5, 1/6, 1/7
1/3, 1/4, 1/5, 1/6, 1/7, 1/8 
1/4, 1/5, 1/6, 1/7, 1/8, 1/9
1/5, 1/6, 1/7, 1/8, 1/9, 1/10
1/6, 1/7, 1/8, 1/9, 1/10, 1/11

Exakte Inverse der Hilbert-Matrix:

36 | -630 | 3360 | -7560 | 7560 | -2772
-630 | 14700 | -88200 | 211680 | -220500 | 83160
3360 | -88200 | 564480 | -1411200 | 1512000| -582120
-7560 | 211680 | -1411200| 3628800 | -3969000 | 1552320
7560 | -220500 | 1512000 | -3969000 | 4410000 | -1746360
-2772 | 83160 | -582120 | 1552320 | -1746360| 698544

Die Funktion 'Inverse_einer_quadratischen_Matrix' verwendet den Datentyp DECIMAL und berechnet diese Werte mit einer Genauigkeit von mehr als 10 Nachkommastellen.
Die gleichen Funktion, aber unter Verwendung des Datentyps DOUBLE, erzielt nur eine Genauigkeit von etwa vier Nachkommastellen.
Die Invertierungsfunktion kommt ohne Typumwandlungen per CDec aus. Das ist möglich, weil sichergestellt ist, dass alle Berechnungen auf der Grundlage von DECIMAL-Werten in der Eingabe-Matrix durchgeführt werden.

Im Beispiel wird die Multiplikation der Hilbert-Matrix mit ihrer Inversen durchgeführt. Im Idealfall sollte dabei die Einheitsmatrix entstehen.
Wer es ausprobiert, wird feststellen, dass die Funktion noch bei Invertierung der Hilbert-Matrix 14. Ordnung eine Lösung errechnet, aus der die Einheitsmatrix mit 8 Nachkommastellen Genauigkeit reproduziert werden kann.
(Es genügt, in der Funktion 'Demo_Matrix_Inverse' die Variable Ordnung = 14 zu setzen. Ab einer Ordnung von 17 stuft die Funktion die Hilbert-Matrix als 'singulär' ein und bricht die Berechnung der Inverse ab.)

Hinweis:
Wer singuläre Matrizen bearbeiten möchte, wird bei der Suche im Internet oder in der Literatur unter den Schlagworten 'generalisierte Inverse' oder 'Moore-Penrose-Inverse' fündig.

Für die Umsetzung von zweidimensionalen Arrays, deren Elemente von Typ DOUBLE sind, auf den Datentyp DECIMAL ist die Funktion 'DoubleArray2DecimalArray' beigefügt.

Zur Demonstration der genaueren Kontrolle des Inhalts der Datenfelder eines VARIANT-Arrays, dient die Funktion 'ConfirmDecimalMatrix'.

Function Demo_Matrix_Inverse()
  ' Beispiel:
  ' Berechnung der Inversen einer Hilbert-Matrix 6. Ordnung
 
  Dim Matrix() As Variant, Inverse() As Variant, Test() As Variant
  Dim Ordnung As Long, i As Long, k As Long
 
  Ordnung = 17
 
  Erzeugung_einer_HilbertMatrix Ordnung, Matrix()
 
  If Not Inverse_einer_Quadratischen_Matrix(Matrix(), Inverse()) Then
    MsgBox "Inverse kann nicht berechnet werden"
  Else
 
    ' Test des Ergebnisses:
    ' theoretisch müsste in der Matrix 'Test' die Einheitsmatrix
    ' stehen (Hauptdiagonale = 1, sonst = 0)
    Matrix_Multiplikation Matrix(), Inverse(), Test()
 
    ' Ausgabe der Ergebnismatrix im Direktfenster
    Debug.Print vbCrLf + "Einheitsmatrix ???"
    For i = 1 To Ordnung
     For k = 1 To Ordnung
       Debug.Print Round(Test(i, k), 8); " ";
     Next k
     Debug.Print ""
    Next i
  End If    
End Function

Function Erzeugung_einer_HilbertMatrix(ByVal Ordnung As Long, _
  Matrix() As Variant) As Boolean
 
  ' Die Funktion erzeugt eine Hilbert-Matrix beliebiger Ordnung
 
  Dim i As Long, k As Long
 
  ' Eingabe prüfen
  If Ordnung < 1 Then Exit Function
 
  ReDim Matrix(1 To Ordnung, 1 To Ordnung)
 
  ' Hilbert-Matrix füllen
  For i = 1 To Ordnung
    For k = 1 To Ordnung
      ' durch Verwendung der Umwandlungsfunktion CDEC
      ' ist explizit sicherzustellen, dass die Brüche mit
      ' max. Genauigkeit in der Matrix gespeichert werden
      Matrix(i, k) = CDec(1) / CDec(k + (i - 1))
    Next k
  Next i
 
  Erzeugung_einer_HilbertMatrix = True
End Function

Function Inverse_einer_Quadratischen_Matrix(Matrix() As Variant, _
  Inverse() As Variant) As Boolean
 
  ' Berechnung der Inversen einer nicht-singulären quadratischen Matrix
  ' gegeben als 2-dimensionales Array mit Feldern des Typs Decimal
  ' Austauschverfahren / Pivotisieren
 
  ' Die Eingabe-Matrix bleibt unverändert
  ' erwartete Deklaration: Matrix(1 To N, 1 To N)
 
  ' Rückgabe: false - falls 'Matrix' singulär, nicht quadratisch
  '                   oder falscher Datentyp in mind. einem Feld
  '           true  - in 'Inverse' steht die Inverse der Matrix
 
  Dim N As Long, M As Long                   ' Array-Dimensionen
  Dim vekx() As Long, veky() As Long         ' Hilfsvektoren
  Dim i As Long, j As Long, k As Long        ' Laufvariable
  Dim ii As Long, ij As Long                 ' innere Laufvariable
  Dim pi As Long, pj As Long                 ' Indices: Pivot-Element
 
  Dim epsilon As Variant                     ' Schranke für Singularität
  Dim eins As Variant                        ' Decimal 1
  Dim PivotElement As Variant                ' Pivot-Element
  Dim faktor As Variant                      ' Skalierungsfaktor
 
  On Error GoTo fehler
 
  ReDim Inverse(0 To 0, 0 To 0)              ' Ausgabe-Matrix löschen
 
  ' Eingabe-Matrix prüfen
 
  ' Größe der Matrix
  N = UBound(Matrix, 1)
 
  ' quadratische Matrix?
  If UBound(Matrix, 2) <> N Then Exit Function
 
  ' Array plausibel deklariert?
  If LBound(Matrix, 1) > 1 Or LBound(Matrix, 2) > 1 Then Exit Function
 
  ' Datentyp DECIMAL ?
  For i = 1 To N
    For k = 1 To N
      If VarType(Matrix(i, k)) <> vbDecimal Then Exit Function
    Next k
  Next i
 
  ' Hilfsvektoren zum Speichern der Zeilen-/Spalten-Vertauschungen
  ReDim vekx(N) As Long, veky(N) As Long
 
  ' Zur Prüfung, ob die Matrix singulär ist
  epsilon = CDec("0,0000000000000000000001")
 
  ' Hilfswert bilden (Decimal 1)
  eins = CDec(1)
 
  ' Übertragung der Eingabe-Matrix
  Inverse() = Matrix()
 
  ' Austauschregister vorbesetzen
  For i = 1 To N
    vekx(i) = 0: veky(i) = 0
  Next i
 
  For i = 1 To N
    ' Suche nach dem Pivotelement
    PivotElement = CDec(0)
 
    For ii = 1 To N
      If vekx(ii) = 0 Then
        For ij = 1 To N
          If veky(ij) = 0 Then
            If Abs(Inverse(ii, ij)) > Abs(PivotElement) Then
              PivotElement = Inverse(ii, ij)
              pi = ii: pj = ij
            End If
          End If
        Next ij
      End If
    Next ii
 
    ' Matrix singulär ?
    If Abs(PivotElement) < epsilon Then
      ' Abbruch wegen singulärer Matrix / Inverse löschen
      ReDim Inverse(0 To 0, 0 To 0)
      Exit Function
    End If
 
    ' Pivot-Indices
    vekx(pi) = pj: veky(pj) = pi
 
    ' Austauschschritt
    For j = 1 To N
      If j <> pi Then
        faktor = Inverse(j, pj) / PivotElement
        For k = 1 To N
          Inverse(j, k) = Inverse(j, k) - Inverse(pi, k) * faktor
        Next k
        Inverse(j, pj) = -faktor
      End If
    Next j
 
    For k = 1 To N
      Inverse(pi, k) = Inverse(pi, k) / PivotElement
    Next k
 
    ' Explizite Typumwandlung erforderlich!
    Inverse(pi, pj) = eins / PivotElement
  Next i
 
  ' Zeilen- und Spaltenvertauschungen aufheben
  For i = 1 To N - 1   
    For M = i To N
      If vekx(M) = i Then Exit For
    Next M
    j = M
    If j <> i Then
      For k = 1 To N
        Swap Inverse(i, k), Inverse(j, k)
      Next k
      vekx(j) = vekx(i)
      vekx(i) = i
    End If
 
    For M = i To N
      If veky(M) = i Then Exit For
    Next M
    j = M
    If j <> i Then
      For k = 1 To N
        Swap Inverse(k, i), Inverse(k, j)
      Next k
      veky(j) = veky(i)
      veky(i) = i
    End If   
  Next i
 
  ' Operation erfolgreich durchgeführt
  Inverse_einer_Quadratischen_Matrix = True
  Exit Function
 
fehler:
  ' ggf. Überlauf abfangen
End Function

Function Swap(a As Variant, b As Variant, _
  Optional Check_Decimal As Boolean = True) As Boolean
 
  ' Hilfsfunktion
  ' der Swap  a <-> b  wird stets durchgeführt
  ' falls DatenTyp unterschiedlich oder nicht-numerisch ist,
  ' wird 'false' gemeldet
  ' optional: falls Datentyp nicht 'Decimal' ist,
  ' wird false gemeldet
 
  Dim c As Variant
 
  c = a: a = b: b = c
 
  If Not VarType(a) = VarType(b) Then Exit Function
  If Not IsNumeric(a) Then Exit Function
 
  If Check_Decimal Then
    If VarType(a) <> vbDecimal Then Exit Function
  End If
 
  Swap = True
End Function

Function Matrix_Multiplikation(a() As Variant, b() As Variant, _
  c() As Variant) As Boolean
 
  ' Multiplikation der Matrix a(x,y) mit b(y,z) in c(x,z)
  ' erwartete Untergrenze der Dimensionen: 1
 
  Dim x As Long, y As Long, xt As Long, z As Long
  Dim i As Long, k As Long, l As Long, ok As Boolean
 
  On Error GoTo fehler
 
  ' Eingabe prüfen
  x = UBound(a, 1): y = UBound(a, 2)
  xt = UBound(b, 1): z = UBound(b, 2)
  If LBound(a, 1) > 1 Or LBound(a, 2) > 1 Then Exit Function
 
  ' Voraussetzung für Multiplikation erfüllt ?
  If y <> xt Then Exit Function
 
  ' Ergebnis-Matrix dimensionieren
  ReDim c(1 To x, 1 To z) As Variant
 
  For i = 1 To x
    For k = 1 To z
      c(i, k) = CDec(0)
      For l = 1 To y
        ' Falls alle Felder der Eingabe-Matrizen
        ' vom Datentyp DECIMAL sind ...
        c(i, k) = c(i, k) + a(i, l) * b(l, k)
        ' sonst ....
        ' c(i, k) = c(i, k) + CDec(a(i, l)) * CDec(b(l, k))
      Next l
    Next k
  Next i
 
  ' Multiplikation erfolgreich durchgeführt
  Matrix_Multiplikation = True
 
  Exit Function
 
fehler:
  ' Fehler wegen
  ' ungeeigneter/nicht-dimensionierter Eingabe-Matrizen
  ' oder wegen Überlauf des Datentyps abfangen
End Function

Function DoubleMatrix2DecimalMatrix(a() As Double, b() As Variant, _
  Optional ByVal ConfirmNumeric As Boolean = False) As Boolean
 
  ' eine zweidimensionale Matrix (a) mit Elementen vom Typ DOUBLE
  ' wird in eine Matrix (b) mit Elementen des Typs VARIANT/DECIMAL
  ' übertragen
 
  Dim i As Long, k As Long  'Laufvariable
  On Error GoTo fehler
 
  ' Ausgabematrix dimensionieren
  ReDim b(LBound(a, 1) To UBound(a, 1), LBound(a, 2) To UBound(a, 2))
 
  For i = LBound(a, 1) To UBound(a, 1)
    For k = LBound(a, 1) To UBound(a, 2)
      ' explizite Typumwandlung jedes Feldes ist erforderlich!
      b(i, k) = CDec(a(i, k))
    Next k
  Next i
 
  ' Übertragung erfolgreich durchgeführt
  DoubleMatrix2DecimalMatrix = True
  Exit Function
 
fehler:
  ' abfangen,
  ' falls Matrix a nicht korrekt dimensioniert
  ' oder Matrix b statisch deklariert ist
End Function

Function ConfirmDecimalMatrix(a() As Variant, b() As Variant, _
  Optional ByVal ConfirmNumeric As Boolean = False) As Boolean
 
  ' eine zweidimensionale Matrix (a) mit Elementen vom Typ VARIANT
  ' wird in eine Matrix (b) mit Elementen des Typs VARIANT/DECIMAL
  ' übertragen
 
  ' falls ConfirmNumeric = true
  ' nur numerische Datentypen werden akzeptiert, sonst Abbruch
 
  ' falls ConfirmNumeric = false
  ' auch STRING mit numerischem Inhalt, BYTE, BOOLEAN oder DATE
  ' wird akzeptiert und in DECIMAL gewandelt
 
  ' Hinweis:
  ' Die Abfrage 'IsNumeric' filtert nur DATE, akzeptiert aber
  ' ansonsten alles, was sich irgendwie numerisch interpretieren
  ' läßt. Sie wird deshalb hier nicht eingesetzt
 
  Dim i As Long, k As Long  ' Laufvariable
  Dim vt As Integer         ' VarType
 
  On Error GoTo fehler
 
  ' Ausgabematrix dimensionieren
  ReDim b(LBound(a, 1) To UBound(a, 1), LBound(a, 2) To UBound(a, 2))
 
  For i = LBound(a, 1) To UBound(a, 1)
    For k = LBound(a, 1) To UBound(a, 2)
      If ConfirmNumeric Then
        ' Nur numerische Datentypen sollen zugelassen sein
        ' sonst erfolgt Abbruch
        vt = VarType(a(i, k))
        If vt <> vbDecimal And vt <> vbCurrency _
          And vt <> vbDouble And vt <> vbSingle _
          And vt <> vbInteger And vt <> vbLong Then
          ' Ausgabe löschen
          ReDim b(0 To 0, 0 To 0)
          Exit Function
        End If
      End If
      ' explizite Typumwandlung jedes Feldes ist erforderlich!
      b(i, k) = CDec(a(i, k))
    Next k
  Next i
 
  ' Übertragung erfolgreich durchgeführt
  ConfirmDecimalMatrix = True
 
   Exit Function
 
fehler:
  ' abfangen,
  ' falls Matrix a nicht korrekt dimensioniert
  ' oder nicht-wandelbare Elemente enthält
  ' oder Matrix b statisch deklariert ist
End Function