STRUCTURE MATRICE : SURCHARGE DES OPÉRATEURS POUR LE CALCUL MATRICIEL, INVERSION TRACE DÉTERMINANT ET AUTRES OPÉRATIONS

Public Structure Matrice
    ' structure pour les matrices à 2 dimensions 

    Private A(,) As Long ' on peut changer de TYPE mais il faut alors le changer (partout)(ou presque)
    Private mDim1 As Long, mDim2 As Long

#Region "Constructeurs et propriétés"

    Public Sub New(ByVal dim1 As Long, ByVal dim2 As Long)
        Dim longA(dim1, dim2) As Long
        A = longA
        mDim1 = dim1
        mDim2 = dim2
    End Sub

    ' ceci est la propriété par défaut pour une matrice
    Default Public Property Item(ByVal dim1 As Long, ByVal dim2 As Long) As Long ' TYPE
        Get
            If dim1 <= mDim1 And dim2 <= mDim2 And dim1 >= 0 And dim2 >= 0 Then
                Return A(dim1, dim2)
            End If
        End Get
        Set(ByVal value As Long)
            If dim1 <= mDim1 And dim2 <= mDim2 And dim1 >= 0 And dim2 >= 0 Then
                A(dim1, dim2) = value
            End If
        End Set
    End Property

    Public ReadOnly Property Dim1() As Long
        Get
            Return mDim1
        End Get
    End Property
    Public ReadOnly Property Dim2() As Long
        Get
            Return mDim2
        End Get
    End Property

    Function Transpose() As Matrice
        'TRANSPOSE MATRICIELLE
        Dim result As New Matrice(mDim2, mDim1)

        'Transpose
        For i As Long = 1 To mDim1
            For j As Long = 1 To mDim2
                result(i, j) = A(j, i)
        Next j, i

        Return result
    End Function

    Public ReadOnly Property Trace() As Long ' TYPE
        Get
            'TRACE DE LA MATRICE

            'Paramètres
            Dim tmp As Long

            'Vérifications
            If mDim1 <> mDim2 Then
                Err.Raise(Number:=513, Source:="MATRICE", Description:="MATRICE NON CARRE")
            Else
                'Addition ou soustraction
                For i As Long = 1 To mDim1
                    tmp += Me(i, i)
                Next i
            End If
        End Get
    End Property

    Public ReadOnly Property Inverse() As Matrice
        Get
            'INVERSION MATRICE METHODE DE GAUSS JORDAN

            'Paramètres
            Dim i As Long, j As Long, k As Double
            Dim Dummy As Double
            Dim B As New Matrice(mDim1, mDim1)

            'Vérifications
            If mDim1 <> mDim2 Then
                Err.Raise(Number:=513, Source:="MATRICE", Description:="MATRICE NON CARRE")
                Return B
            Else
                'Algo

                For k = 1 To mDim1
                    B(k, k) = 1
                Next k
                For k = 1 To mDim1
                    Dummy = A(k, k)
                    For j = 1 To mDim1
                        A(k, j) = A(k, j) / Dummy
                        B(k, j) = B(k, j) / Dummy
                    Next j
                    For i = 1 To mDim1
                        If i <> k Then
                            Dummy = A(i, k)
                            For j = 1 To mDim1
                                A(i, j) = A(i, j) - Dummy * A(k, j)
                                B(i, j) = B(i, j) - Dummy * B(k, j)
                            Next j
                        End If
                    Next i
                Next k

                'renvoi
                Return B
            End If
        End Get
    End Property

    Public ReadOnly Property Determinant() As Long
        Get
            'DETERMINANT D'UNE MATRICE PAR TRIANGULATION DE GAUSS

            'Paramètres
            Dim i As Long, j As Long, k As Double
            Dim Factor As Double, Sum As Double = 1

            'Vérifications
            If mDim1 <> mDim2 Then
                Err.Raise(Number:=513, Source:="MATRICE", Description:="MATRICE NON CARRE")
                Return 0
            Else
                'Algo élimination
                For k = 1 To mDim1 - 1
                    For i = k + 1 To mDim1
                        Factor = A(i, k) / A(k, k)
                        For j = k + 1 To mDim1
                            A(i, j) = A(i, j) - Factor * A(k, j)
                        Next j
                    Next i
                Next k

                'Produit diagonale
                For k = 1 To mDim1
                    Sum = Sum * A(k, k)
                Next k
            End If

        End Get
    End Property
#End Region

#Region "Opérateurs surchargés"

    Public Shared Operator *(ByVal A As Matrice, ByVal B As Matrice) As Matrice
        'PRODUIT MATRICE : R=A*B

        'Paramètres
        Dim A1 As Long = A.Dim1, A2 As Long = A.Dim2
        Dim B1 As Long = B.Dim1, B2 As Long = B.Dim2
        Dim result As New Matrice(A1, B2)
        Dim s As Long

        'Vérification
        If A2 <> B1 Then
            Err.Raise(Number:=513, Source:="MATRICE", Description:="MATRICE NON CARRE")
            Return result
        Else
            'Produit
            For i As Long = 1 To A1
                For j As Long = 1 To B2
                    s = 0
                    For k As Long = 1 To B1
                        s = s + A(i, k) * B(k, j)
                    Next k
                    result(i, j) = s
            Next j, i

            'renvoi
            Return result
        End If
    End Operator

    Public Shared Operator *(ByVal A As Matrice, ByVal B As Long) As Matrice
        'PRODUIT MATRICE PAR UN SCALAIRE : A=cR

        'Paramètres
        Dim A1 As Long = A.Dim1, A2 As Long = A.Dim2
        Dim result As New Matrice(A1, A2)

        'Vérification
        If A1 <> A2 Then
            Err.Raise(Number:=513, Source:="MATRICE", Description:="MATRICE NON CARRE")
            Return result
        Else
            'Produit scalaire
            For i As Long = 1 To A1
                For j As Long = 1 To A2
                    result(i, j) = B * A(i, j)
            Next j, i

            'renvoi
            Return result
        End If

    End Operator

    Public Shared Operator /(ByVal A As Matrice, ByVal B As Matrice) As Matrice
        Return A * B.Inverse
    End Operator

    Public Shared Operator ^(ByVal A As Matrice, ByVal Expo As Long) As Matrice
        'PUISSANCE ENTIERE D'UNE MATRICE

        'Paramètres
        Dim A1 As Long = A.Dim1, A2 As Long = A.Dim2
        Dim B As New Matrice(A1, A2)

        'Vérification
        If A1 <> A2 Then
            Err.Raise(Number:=513, Source:="MATRICE", Description:="MATRICE NON CARRE")
            Return B
        Else
            'Matrice unité
            Dim i As Long
            For i = 1 To A1 : B(i, i) = 1 : Next i

            'Cas Triviaux
            If Expo = 0 Then Return B
            'Si puissance négative
            If Expo < 0 Then
                A = A.Inverse
                Expo = Math.Abs(Expo)
            End If
            If Expo = 1 Then Return A

            'Algo exponentiation rapide
            Do
                If Expo And 1 Then B = A * B
                Expo = Expo \ 2
                A = A * A
            Loop While Expo > 1
            Return B * A
        End If

    End Operator

    Public Shared Operator +(ByVal A As Matrice, ByVal B As Matrice) As Matrice
        'ADDITION MATRICE : R=A+B

        'Paramètres
        Dim A1 As Long = A.Dim1, A2 As Long = A.Dim2
        Dim B1 As Long = B.Dim1, B2 As Long = B.Dim2
        Dim result As New Matrice(A1, B2)

        'Vérifications
        If A1 <> B1 Or A2 <> B2 Then
            Err.Raise(Number:=513, Source:="MATRICE", Description:="MATRICE NON CARRE")
            Return result
        Else
            'Addition 
            For i As Long = 1 To A1
                For j As Long = 1 To B2
                    result(i, j) = A(i, j) + B(i, j)
            Next j, i

            'renvoi
            Return result
        End If
    End Operator

    Public Shared Operator -(ByVal A As Matrice, ByVal B As Matrice) As Matrice
        'SOUSTRACTION MATRICE : R=A-B

        'Paramètres
        Dim A1 As Long = A.Dim1, A2 As Long = A.Dim2
        Dim B1 As Long = B.Dim1, B2 As Long = B.Dim2
        Dim result As New Matrice(A1, B2)

        'Vérifications
        If A1 <> B1 Or A2 <> B2 Then
            Err.Raise(Number:=513, Source:="MATRICE", Description:="MATRICE NON CARRE")
            Return result
        Else
            'soustraction
            For i As Long = 1 To A1
                For j As Long = 1 To B2
                    result(i, j) = A(i, j) - B(i, j)
            Next j, i

            'renvoi
            Return result
        End If
    End Operator

    Public Shared Operator -(ByVal A As Matrice) As Matrice
        'NEGATION MATRICE : R=-B

        'Paramètres
        Dim A1 As Long = A.Dim1, A2 As Long = A.Dim2
        Dim result As New Matrice(A1, A2)

        'Négation
        For i As Long = 1 To A1
            For j As Long = 1 To A2
                result(i, j) = -A(i, j)
        Next j, i

        'renvoi
        Return result
    End Operator

#End Region

#Region "Autres méthodes"

    Public Sub PermuteLignes(ByVal L1 As Long, ByVal L2 As Long)
        'PERMUTE LES LIGNES L1 ET L2

        'Paramètres
        Dim v As Double

        'Vérifications
        If mDim1 < L1 Or mDim1 < L2 Then
            Err.Raise(Number:=514, Source:="MATRICE", Description:="MANIPULATION LIGNE : MATRICE PAS ASSEZ GRANDE")
        Else
            'Permutation
            For j As Long = 1 To mDim2
                v = A(L1, j) : A(L1, j) = A(L2, j) : A(L2, j) = v
            Next j
        End If
    End Sub

    Public Sub PermuteColonnes(ByVal C1 As Long, ByVal C2 As Long)
        'PERMUTE LES COLONNES C1 ET C2

        'Paramètres
        Dim v As Double

        'Vérifications
        If mDim2 < C1 Or mDim2 < C2 Then
            Err.Raise(Number:=514, Source:="MATRICE", Description:="MANIPULATION COLONNE : MATRICE PAS ASSEZ GRANDE")
        Else
            'Permutation
            For j As Long = 1 To mDim1
                v = A(j, C1) : A(j, C1) = A(j, C2) : A(j, C2) = v
            Next j
        End If
    End Sub

    Public Function Pivot(ByVal R As Matrice) As Matrice
        'PIVOTage DE MATRICES d'un système Ax=r

        'Paramètres
        Dim R1 As Long = R.Dim1
        Dim lDummy As Double, lMaxa As Double, lPivot As Double

        'Vérifications
        If mDim1 <> mDim2 Or mDim1 <> R1 Or R.Dim2 <> 1 Then
            Err.Raise(Number:=513, Source:="MATRICE", Description:="MATRICE NON CARRE")
        Else
            'Recherche pivot
            For k As Long = 1 To mDim1
                lPivot = k
                lMaxa = Math.Abs(A(k, k))
                For i As Long = k + 1 To mDim1
                    lDummy = Math.Abs(A(i, k))
                    If lDummy > lMaxa Then
                        lMaxa = lDummy
                        lPivot = i
                    End If
                Next i
                If lPivot <> k Then
                    For j As Long = k To mDim1
                        lDummy = A(lPivot, j)
                        A(lPivot, j) = A(k, j)
                        A(k, j) = lDummy
                    Next j
                    lDummy = R(lPivot, 1)
                    R(lPivot, 1) = R(k, 1)
                    R(k, 1) = lDummy
                End If
            Next k
        End If
        Return R
    End Function

    Function GAUSS_SEIDEL(ByVal R As Matrice, Optional ByVal Lambda As Long = 1, _
    Optional ByVal Maxit As Long = 150, Optional ByVal eps As Double = 0.0000000000001) As Matrice
        'RESOLUTION DE GAUSS-SEIDEL

        'Paramètres
        Dim R1 As Long = R.Dim1
        Dim i As Long, j As Long
        Dim Dummy As Double, Sum As Double
        Dim epsa As Double, iter As Long, Old As Double
        Dim x As New Matrice(mDim1, 1)

        'Vérifications
        If mDim1 <> mDim2 Or mDim1 <> R1 Or R.Dim2 <> 1 Then
            Err.Raise(Number:=513, Source:="MATRICE", Description:="MATRICE NON CARRE")
        ElseIf Lambda < 0 Or Lambda > 2 Then
            Err.Raise(Number:=513, Source:="MATRICE", Description:="RELAXATION ERONNEE")
        Else
            'Algo
            epsa = 1.1 * eps
            For i = 1 To mDim1
                Dummy = A(i, i)
                For j = 1 To mDim1
                    A(i, j) = A(i, j) / Dummy
                Next j
                R(i, 1) = R(i, 1) / Dummy
            Next i

            iter = 0

            Do While (iter < Maxit And epsa > eps)
                iter = iter + 1
                For i = 1 To mDim1
                    Old = x(i, 1)
                    Sum = R(i, 1)
                    For j = 1 To mDim1
                        If i <> j Then Sum = Sum - A(i, j) * x(j, 1)
                    Next j
                    x(i, 1) = Lambda * Sum + (1 - Lambda) * Old
                    If x(i, 1) <> 0 Then epsa = Math.Abs((x(i, 1) - Old) / x(i, 1)) * 100
                Next i
            Loop
        End If

        'renvoi
        Return x
    End Function

    Function Gauss(ByVal R As Matrice) As Matrice
        'ELIMINATION DE GAUSS

        'Paramètres
        Dim R1 As Long = R.Dim1
        Dim i As Long, j As Long, k As Double
        Dim lFactor As Double, lSum As Double
        Dim x As New Matrice(mDim1, 1)

        'Vérifications
        If mDim1 <> mDim2 Or mDim1 <> R1 Or R.Dim2 <> 1 Then
            Err.Raise(Number:=513, Source:="MATRICE", Description:="MATRICE NON CARRE")
        Else
            'Algo élimination
            For k = 1 To mDim1 - 1
                For i = k + 1 To mDim1
                    R = Me.Pivot(R) ' modifie la matrice ME (A) et renvoie R(?)
                    If A(k, k) = 0 Then
                        Err.Raise(Number:=513, Source:="MATRICE", Description:="MATRICE NON INVERSIBLE")
                    End If
                    lFactor = A(i, k) / A(k, k)
                    For j = k + 1 To mDim1
                        A(i, j) = A(i, j) - lFactor * A(k, j)
                    Next j
                    R(i, 1) = R(i, 1) - lFactor * R(k, 1)
                Next i
            Next k

            'Algo remontée
            If A(mDim1, mDim1) = 0 Then
                Err.Raise(Number:=513, Source:="MATRICE", Description:="MATRICE SANS BASE PRINCIPALE / SOLUTIONS MULTIPLES")
            Else
                x(mDim1, 1) = R(mDim1, 1) / A(mDim1, mDim1)
                For i = mDim1 - 1 To 1 Step -1
                    lSum = 0
                    For j = i + 1 To mDim1
                        lSum = lSum + A(i, j) * x(j, 1)
                    Next j
                    x(i, 1) = (R(i, 1) - lSum) / A(i, i)
                Next i
            End If
        End If

        'renvoi
        Gauss = x

    End Function

#End Region

End Structure


'############# test, dans un MODULE SEPARE #############
    'taille
    Public Const t As Long = 3

    Public Sub DebugPrintMatrice(ByVal A As Matrice, Optional ByVal label As String = "")
        Dim str As String
        For i As Long = 1 To A.Dim1
            str = ""
            For j As Long = 1 To A.Dim2
                str &= "," & A.Item(i, j)
            Next j
            Debug.Print(label & i & "(" & Mid(str, 2) & ")")
        Next i
    End Sub

    Public Function InitMatrice() As Matrice
        Dim A As New Matrice(t, t)
        'remplis deux matrices
        A.Item(1, 1) = 1 : A.Item(1, 2) = 1 : A.Item(1, 3) = -1
        A.Item(2, 1) = 2 : A.Item(2, 2) = 1 : A.Item(2, 3) = 1
        A.Item(3, 1) = 4 : A.Item(3, 2) = -4 : A.Item(3, 3) = 2
        Return A
    End Function
    Public Sub Main()

        'définitions des MATRICES
        Dim A As New Matrice(t, t), B As New Matrice(t, t), R As New Matrice(t, 1)
        A = InitMatrice()

        Debug.Print("MATRICE A")
        DebugPrintMatrice(A, "A")

        B = A * 1
        Debug.Print("MATRICE B = MATRICE A")
        DebugPrintMatrice(B, "B")

        'Calculs
        Debug.Print("DETERMINANT A=" & A.Determinant)
        Debug.Print("DETERMINANT B=" & B.Determinant)

        Debug.Print("TRACE A=" & A.Trace)
        Debug.Print("TRACE B=" & B.Trace)

        A = A ^ 2
        Debug.Print("MATRICE A^(-2)")
        DebugPrintMatrice(A, "A")

        'résolution SYSTEME EQUATION AX=C
        Dim x As Matrice, c As New Matrice(t, 1)
        A = InitMatrice()
        c.Item(1, 1) = 3 : c.Item(2, 1) = 5 : c.Item(3, 1) = 1

        x = A.Gauss(c)
        Debug.Print("RESOLUTION EQUATION PAR GAUSS")
        DebugPrintMatrice(x, "X")

        A = InitMatrice()
        c.Item(1, 1) = 3 : c.Item(2, 1) = 5 : c.Item(3, 1) = 1

        Debug.Print("RESOLUTION EQUATION PAR GAUSS-SEIDEL")
        x = A.GAUSS_SEIDEL(c, 0.5)
        DebugPrintMatrice(x, "X")

    End Sub