ProfQuiz 03 – Pour les profs...

Pour les détails quant à cette démarche et pour des liens vers les autres questions : index.html#profquiz.

Pour un rappel des principes : index.html#profquiz_principes

Pour cette question, les solutionnaires suivants sont disponibles :

Question 03

Cette question fut proposée le 30 mai 2013.

Un(e) étudiant(e) suit le cours d'algèbre matricielle et souhaite écrire une classe Matrice représentant une matrice de même que la fonction permettant de multiplier deux matrices pour obtenir la matrice résultante. Pouvez-vous l'aider? On vise ici une solution élégante, mais pas nécessairement une version archi optimisée (pour ce type d'application, on peut aller très loin dans l'optimisation, mais notre étudiant(e) veut surtout essayer de bien comprendre son cours de maths!). Implémentez le nécessaire pour en arriver à (a) une classe Matrice et (b) une fonction pour multiplier deux matrices et obtenir le résultat. Ne lui conseillez pas du code tout fait par un tiers : elle/ il souhaite programmer pour mieux comprendre!

Solution possible avec C#

Un gros merci à Nicolas Chourot qui a proposé l'essentiel de ce qui suit, et m'a ainsi économisé du temps. Dans ses propres mots :

« Voici une solution produite sur un coin de table. Elle n'est pas décorée de validation des bornes.

Il y a certainement place à l'amélioration.

J'ai choisi un tableau dynamique de tableau dynamique interne ([ ][ ]) plutôt qu'un tableau dynamique bidimensionnel ([ , ]) pour m'éviter d'avoir à maintenir les propriétés NbRangées et NbColonnes. Le hic pédagogique est que l'étudiant devra vraiment comprendre la différence entre matrice[x][y] et matrice[x,y] (surchargé dans la classe Matrice). Il faudra aussi expliquer le mécanisme de passage de paramètres en nombre variable. »

Sa proposition suit :

// Classe Matrice
using System;

namespace MatriceSimple
{
    class Matrice
    {
        private double[][] matrice;
        public Matrice(int nbRangées, int nbColonnes)
        {
            matrice = new double[nbRangées][];
            for (int ran = 0; ran < nbRangées; ran++)
            {
                matrice[ran] = new double[nbColonnes];
            }
       }
        public int NbRangées
        {
            get { return matrice.Length; }
        }
        public int NbColonnes
        {
            get { return matrice[0].Length; }
        }
        public double this[int ran, int col]
        {
            set { matrice[ran][col] = value; }
            get { return matrice[ran][col]; }
        }
        public void Set(int numRan, params double[] listVal)
        {
            for (int col = 0; col < listVal.Length; col++)
            {
                matrice[numRan][col] = listVal[col];
            }
        }
        public void Affiche()
        {
            for (int ran = 0; ran < NbRangées; ran++)
            {
                Console.Write("|");
                for (int col = 0; col < NbColonnes; col++)
                {
                    Console.Write("\t {0}", matrice[ran][col]);
                }
                Console.WriteLine("\t|");
            }
        }
        public static Matrice operator * (Matrice matriceX, Matrice matriceY)
        {
            if (matriceX.NbColonnes != matriceY.NbRangées)
            {
                throw new Exception("matrices incompatibles pour la multiplication!");
            }
            Matrice matriceXY = new Matrice(matriceX.NbRangées, matriceY.NbColonnes);
            for (int ran = 0; ran < matriceX.NbRangées; ran++)
            {
                for (int col = 0; col < matriceY.NbColonnes; col++)
                {
                    double prodScal = 0;
                    for (int pos = 0; pos < matriceY.NbRangées; pos++)
                    {
                        prodScal += matriceX[ran,pos] * matriceY[pos,col];
                    }
                    matriceXY[ran,col] = prodScal;
                }
            }
            return matriceXY;
        }
    }
    class Program
    {
        static void Main(string[] args)
        {
            Matrice matriceX = new Matrice(4, 3);
            matriceX.Set(0, 1, 1, 1);
            matriceX.Set(1, 2, 2, 2);
            matriceX.Set(2, 3, 3, 3);
            matriceX.Set(3, 4, 4, 4);
            Matrice matriceY = new Matrice(3, 4);
            matriceY.Set(0, 1, 2, 3, 4);
            matriceY.Set(1, 1, 2, 3, 4);
            matriceY.Set(2, 1, 2, 3, 4);
            Matrice matriceXY = matriceX * matriceY;
            Console.WriteLine("La matrice suivante");
            matriceX.Affiche();
            Console.WriteLine("multipliée par celle ci");
            matriceY.Affiche();
            Console.WriteLine("donne");
            matriceXY.Affiche();
        }
    }
}

Le support par C# de la surcharge d'opérateurs permet d'en arriver à une écriture naturelle pour la multiplication de matrices dans le code client.

À l'exécution, ce programme donnera (prudence : je ne sais pas comment les divers fureteurs traiteront l'affichage des tabulations, en soi correctes dans un programme console, mais l'affichage d'origine est tout à fait correct) :

La matrice suivante
|        1       1       1      |
|        2       2       2      |
|        3       3       3      |
|        4       4       4      |
multipliée par celle-ci
|        1       2       3       4      |
|        1       2       3       4      |
|        1       2       3       4      |
donne
|        3       6       9       12     |
|        6       12      18      24     |
|        9       18      27      36     |
|        12      24      36      48     |

Quelques remarques supplémentaires :

Nicolas a aussi proposé une autre solution pour laquelle la signature de l'accès aux éléments est plus homogène, ce qui peut réduire les obstacles à la compréhension :

« Voici une autre solution, dérivée de ma première. Elle est selon moi plus élégante.

En découvrant l'existence la méthode Array.GetLength(int dimNumber), l'approche « tableau de tableau » n'était plus requise pour évaluer NbRangées et NbColonnes. Aussi le constructeur dans cette solution me semble plus intuitif pour l'utilisateur. Plus besoin de paramètres pour l'initialisation et plus besoin d'expliquer la différence entre [ ][ ] et [ , ]. »

Le code suit :

using System;

namespace MatriceSimple
{
    class Matrice
    {
        private double[,] matrice;
        public Matrice(int nbRangées, int nbColonnes)
        {
            matrice = new double[nbRangées, nbColonnes];
        }
        public Matrice(double [,] tableau_2D)
        {
            matrice = tableau_2D;
        }
        public int NbRangées
        {
           get { return matrice.GetLength(0); }
        }
        public int NbColonnes
        {
           get { return matrice.GetLength(1); }
        }
        public double this[int ran, int col]
        {
            set { matrice[ran,col] = value; }
            get { return matrice[ran,col]; }
        }
        public void Affiche()
        {
            for (int ran = 0; ran < NbRangées; ran++)
            {
                Console.Write("|");
                for (int col = 0; col < NbColonnes; col++)
                {
                    Console.Write("\t {0}", matrice[ran,col]);
                }
                Console.WriteLine("\t|");
            }
        }
        public static Matrice operator *(Matrice matriceX, Matrice matriceY)
        {
            if (matriceX.NbColonnes != matriceY.NbRangées)
            {
                throw new Exception("matrices incompatibles pour la multiplication!");
            }
            Matrice matriceXY = new Matrice(matriceX.NbRangées, matriceY.NbColonnes);
            for (int ran = 0; ran < matriceX.NbRangées; ran++)
            {
                for (int col = 0; col < matriceY.NbColonnes; col++)
                {
                    double prodScal = 0;
                    for (int pos = 0; pos < matriceY.NbRangées; pos++)
                    {
                        prodScal += matriceX[ran, pos] * matriceY[pos, col];
                    }
                    matriceXY[ran, col] = prodScal;
                }
            }
            return matriceXY;
        }
    }
    class Program
    {
        static void Main(string[] args)
        {
            double [,] X = new double [,] {
               {1, 1, 1},
               {2, 2, 2},
               {3, 3, 3},
               {4, 4, 4}
            };
            Matrice matriceX = new Matrice(X);
            double [,] Y = new double [,] {
               {1, 2, 3, 4},
               {1, 2, 3, 4},
               {1, 2, 3, 4}
            };
            Matrice matriceY = new Matrice(Y);
            Matrice matriceXY = matriceX * matriceY;
            Console.WriteLine("la matrice suivante");
            matriceX.Affiche();
            Console.WriteLine("multipliée par celle ci");
            matriceY.Affiche();
            Console.WriteLine("donne");
            matriceXY.Affiche();
        }
    }
}

Nicolas complétait cette proposition comme suit :

« Il serait possible mais complexe d'enseigner l'implantation de l'interface IEnumerable pour simplifier encore plus la vie de l'utilisateur avec la construction suivante d'un objet de classe Matrice :

Matrice matriceX = new Matrice(){
   {1, 1, 1},
   {2, 2, 2},
   {3, 3, 3},
   {4, 4, 4}
};

... plutôt que le détour suivant (qui, néanmoins, est plus simple et intuitif que dans ma solution version 1) :

double [,] X = new double [,] {
   {1, 1, 1},
   {2, 2, 2},
   {3, 3, 3},
   {4, 4, 4}
};
Matrice matriceX = new Matrice(X);

... ou encore :

Matrice matriceX = new Matrice(
   new double [,] {
      {1, 1, 1},
      {2, 2, 2},
      {3, 3, 3},
      {4, 4, 4}
   }
);

Il s'agit de savoir ce que l'on cherche à enseigner. »

Comme Nicolas et moi en avions discuté ultérieurement, il serait bien sûr possible d'appliquer l'écriture suivante (synopsis seulement) :

// ...
   class Matrice
   {
      double[,] elems;
      public Matrice(double[,] zeElems)
      {
         elems = zeElems;
      }
      // ...
   }
   class Program
   {
      static void Main(string[] args)
      {
         Matrice mat = new Matrice(
            new double[,] {
               { 1.0, 2.0, 3.0 },
               { 4.0, 5.0, 6.0 },
               { 7.0, 8.0, 9.0 }
            }
         );
         // ...
      }
   }
// ...

Un certain nombre de problèmes demeurent, mais pour une implémentation naïve mais fonctionnelle, tel que demandé, c'est pas mal du tout.

Solution possible avec C++

Pour une proposition en C++, j'ai pensé à ceci. Notez que ce programme ne compilera pas avec Visual Studio 2012 (qui ne supporte pas std::initializer_list, tristement) mais compile sans problème avec Visual Studio 2013 Preview, de même qu'avec les versions à jour de Clang et de g++. J'ai inséré des annotations pour celles et ceux qui seraient moins familièes ou familiers avec C++ 11 :

Les bibliothèques sur lesquelles reposera cette proposition sont strictement des bibliothèques standards. Portez une attention particulière à <initializer_list>, une particularité de C++ 11 qui rend l'initialisation des primitifs et des objets plus homogène que par le passé; c'est grâce à cette fonctionnalité que nous pourrons ici régler l'irritant relevé par Nicolas dans l'implémentation C#.

#include <vector>
#include <cassert>
#include <istream>
#include <ostream>
#include <initializer_list>
#include <algorithm>
using namespace std;

J'ai placé ma classe Matrice dans un espace nommé (ici : vanille), pour être plus facilement en mesure de passer d'une implémentation à l'autre dans le code client si le coeur m'en dit... Dans la mesure bien sûr où les interfaces (au sens large de méthodes et types publics) de chacune des classes Matrice concordent, quel que soit leur espace nommé.

J'ai choisi d'exposer les types internes et publics size_type, pour décrire le type de la largeur et de la hauteur d'une Matrice, et value_type pour le type des valeurs dans chacune de ses cellules. Ce choix tient du fait que j'ai défini ici une Matrice contenant des double et utilisant des std::vector comme substrat, mais ce sont des choix qui pourraient changer – par exemple, je pourrais généraliser en offrant une Matrice<T> pour permettre des trucs tels que Matrice<std::complex> pour les gens qui en ont besoin; en exposant les alias appropriés pour ces types au code client, cela permet d'écrire des programmes qui survivront mieux aux changements d'implémentation.

Les méthodes nb_lignes() et nb_colonnes() permettent de consulter les dimensions de la Matrice sans pouvoir les modifier, comme le font les propriétés NbRangées et NbColonnes dans la proposition en C#.

Une différence dans les choix d'implémentation est que j'ai défini ici un opérateur [] prenant un numéro de ligne comme intrant, et retournant une ligne entière comme extrant. La ligne ayant son propre opérateur [] pour accéder à une cellule, il est possible d'uitliser la notation m[i][j] pour accéder à la cellule à la ligne i et à la cellule j de la Matrice nommée m. Notez toutefois que je n'ai offert aucune validation des bornes dans un cas comme dans l'autre – un débordement de capacité aura donc les conséquences habituelles en C++, menant à un vilain Undefined Behavior.

J'ai choisi d'offrir explictement deux constructeurs :

  • un constructeur paramétrique muni de deux paramètres, soit le nombre de lignes et le nombre de colonnes. La validation est quelque peu brutale : les valeurs de ces deux paramètres sont nécessairement supérieures ou égales à zéro (le type size_type choisi ici étant un type entier non-signé), et une dimension nulle mènera à un arrêt immédiat du programme (violation de précondition). Notez que le constructeur d'un vector<T> ayant un nombre non-nul d'éléments initalise chacun à la valeur par défaut du type T, et que cette valeur est le zéro du type dans le cas des primitifs. Conséquemment, ici, les éléments seront initialisés à 0.0 car value_type est double; et
  • un constructeur paramétrique acceptant une initializer_list contenant des initializer_list<double>, ce qui permettra d'utiliser la notation entre accolades à l'initialisation d'une matrice. Ici, il faut noter qu'une familiarité avec les itérateurs est requise, mais c'est le cas en général en C++ pour tout ce qui touche à la bibliothèque standard, aux conteneurs et aux algorithmes. Dans ce constructeur, l'utilisation de reserve() est une optimisation; enlever cette ligne ne change rien au comportement observable du programme.

Comme c'est toujours le cas en C++, la Sainte-Trinité est offerte par défaut sauf si nous la refusons explicitement. Conséquemment, le constructeur de copie, l'affectation et le destructeurs sont définis en termes de ces opérations sur les attributs de la classe, et il se trouve que c'est le cas pour un vecteur dont les éléments offrent correctement ces trois opérations (donc pour double, pour vector<double> et pour vector<vector<double>>). C'est un gros avantage de C++ dans un programme comme celui-ci, comme nous le verrons plus bas.

namespace vanille
{
   class Matrice
   {
   public:
      using value_type = double;
   private:
      vector<vector<value_type>> vals_;
   public:
      using size_type = vector<value_type>::size_type;
      Matrice(size_type lignes, size_type colonnes)
         : vals_(lignes, vector<value_type>(colonnes))
      {
         assert(lignes && colonnes);
      }
      Matrice(initializer_list<initializer_list<value_type>> init)
      {
         vals_.reserve(init.size());
         for (auto & v : init)
            vals_.push_back(vector<value_type>(begin(v), end(v)));
      }
      size_type nb_lignes() const
         { return vals_.size(); }
      size_type nb_colonnes() const
         { return vals_.front().size(); }
      vector<value_type>& operator [](size_type n)
         { return vals_[n]; }
      const vector<value_type>& operator [](size_type n) const
         { return vals_[n]; }
   };

Dans le même espace nommé que la classe Matrice, j'ai ajouté des fonctions globales hauteur() et largeur() toutes deux applicables à une instance de Matrice. C'est une pratique courante en C++ que d'élargir l'interface d'une classe par des fonctions globales; ceci permet de surcharger ces fonctions pour différents types dans l'optique d'offrir une API homogène et non-intrusive pour le code client, même si diverses classes représentant une Matrice sont utilisées dans le même programme; j'y reviendrai un peu plus bas en décrivant la fonction est_carree().

La fonction de projection d'une Matrice sur un flux en sortie est banale. Je ne l'ai pas adaptée pour avoir le même affichage que celui apparaissant dans le code C# plus haut, mais les changements pour avoir la même sortie dans les deux cas sont banals et cosmétiques alors je vous invite à les implémenter si le coeur vous en dit.

// ...
   Matrice::size_type hauteur(const Matrice &mat)
      { return mat.nb_lignes(); }
   Matrice::size_type largeur(const Matrice &mat)
      { return mat.nb_colonnes(); }
   ostream &operator<<(ostream &os, const Matrice &mat)
   {
      for (Matrice::size_type i = 0; i < hauteur(mat); ++i)
      {
         for (Matrice::size_type j = 0; j < largeur(mat); ++j)
            os << mat[i][j] << ' ';
         os << endl;
      }
      return os;
   }
}

Pour que le code client qui suit fonctionne en utilisant la classe Matrice (et ses outils accessoires) décrite ci-dessus, il suffit d'ajouter une instruction using. Notez que ceci rend aussi implictement accessibles (de manière contextuelle) les fonctions qui utilisent une vanille::Matrice, comme vanille::largeur() par exemple.

// ...
using vanille::Matrice;

L'opération de multiplication de deux instances de Matrice, m0 et m1, est une expression classique de l'algorithme naïf en ce sens. Sa complexité est si m0 est une matrice et si m1 est une matrice ; conséquemment, la complexité pour la mutliplication de deux matrices carrées est .

La matrice résultante est nommée resultat et est retournée par valeur... Notez qu'une optimisation importante serait possible ici en introduisant dans Matrice le support à la sémantique de mouvement, mais nous avions convenu de viser une solution simple plutôt qu'une solution à haute performance, alors je me suis abstenu.

J'ai offert une fonction produit_matriciel(m0,m1) qui délègue à la multiplication, mais je dois vous avouer que je ne suis plus certain des termes (les termes anglais pour les opérations sur les matrices tendent à différer des termes français, alors mes excuses si j'ai pris le mauvais nom... ça m'arrive souvent!).

// ...
Matrice operator*(const Matrice &m0, const Matrice &m1)
{
   assert(largeur(m0) == hauteur(m1));
   Matrice resultat(hauteur(m0), largeur(m1));
   for (Matrice::size_type i = 0; i < hauteur(resultat); ++i)
      for (Matrice::size_type j = 0; j < largeur(resultat); ++j)
         for (Matrice::size_type k = 0; k < largeur(m0); ++k)
            resultat[i][j] += m0[i][k] * m1[k][j];
   return resultat;
}
Matrice produit_matriciel(const Matrice &m0, const Matrice &m1)
   { return m0 * m1; }

Là où la multiplication de deux matrices n'est pas commutative, la multiplication d'une matrice par un scalaire l'est. Conséquemment, j'ai implémenté cette opération des deux manières, en réalisant la multiplication d'un scalaire par une matrice à travers une délégation à une multiplication d'une matrice par un scalaire; si la vitesse d'exécution vous préoccupe, sachez qu'il est banal pour le compilateur d'appliquer le Function Inlining au cas qui se limite à une délégation, et donc de réaliser cette simplification sans qu'il n'y ait de coût à l'exécution.

Comme mentionné plus haut, le retour par valeur pourrait être remplacé par un retour par mouvement avec très peu d'effort, mais tel n'est pas l'objectif ici.

// ...
Matrice operator*(Matrice mat, Matrice::value_type val)
{
   for (Matrice::size_type i = 0; i < hauteur(mat); ++i)
      for (Matrice::size_type j = 0; j < largeur(mat); ++j)
         mat[i][j] *= val;
   return mat;
}
Matrice operator*(Matrice::value_type val, Matrice mat)
   { return mat * val; }

Pour le plaisir, j'ai implémenté un prédicat unaire est_carree() permettant de déterminer si une Matrice donnée est carrée ou non, de même qu'une fonction de fabrication d'une matrice identité d'une certaine taille.

Le prédicat a été écrit pour montrer comment il est possible d'enrichir l'interface de base d'une classe de manière non-intrusive par l'ajout de fonctions globales se limitant à exploiter des services plus primitifs de cete interface. Ceci permet de garder les classes plus simples, plus circonscrites, sans perte de flexibilité. Notez que C# permet d'utiliser ce qu'on y nomme des méthodes d'extension, enrichissant de manière analogue l'interface d'une classe par des services externes, mais utilisables par le code client comme s'il s'agissait de méthodes. Les moeurs des deux langages diffèrent mais les résultats sont les mêmes.

// ...
bool est_carree(const Matrice &mat)
   { return hauteur(mat) == largeur(mat); }
Matrice creer_matrice_identite(Matrice::size_type n)
{
   using size_type = Matrice::size_type;
   Matrice mat(n, n);
   for (size_type i = 0; i < n; ++i)
      for (size_type j = 0; j < n; ++j)
         mat[i][j] = i == j ? 1.0 : 0.0;
   return mat;
}

Le code client proposé est celui à droite. Quelques remarques :

  • l'initialisation d'une Matrice se fait à l'aide de la même syntaxe que celle applicable aux tableaux, ceci grâce à initializer_list;
  • une Matrice s'affiche comme on afficherait un entier, respectant la maxime de Scott Meyers : « Whenever possible, do as the ints do », ou en français « Chaque fois que cela s'avère raisonnable, faites en sorte que vos classes s'utilisent comme des int ». Mieux vaut viser la simplicité;
  • aucune allocation dynamique de mémoire n'est visible, toutes les opérations de ce genre étant encapsulées dans le substrat (std::vector) de Matrice;
  • aucune fuite de mémoire ou de ressources, car le destructeur (implicite) de Matrice assure la destruction (implicite) de ses attributs, tous responsables de leurs propres ressources.
// ...
#include <iostream>
int main()
{
   Matrice mat = {
      { 1.0, 1.0, 1.0, 1.0 },
      { 1.0, 1.0, 1.0, 1.0 },
      { 1.0, 1.0, 1.0, 1.0 }
   };
   cout << mat << endl;
   mat = mat * 2.5;
   cout << mat << endl;
   auto mid4x4 = creer_matrice_identite(4);
   cout << mid4x4 << endl;
   //
   // boum!
   // cout << mid4x4 * mat << endl;
   //
   cout << mat * mid4x4 << endl;
   Matrice m2x2a = {
      { 1.0, 2.0 },
      { 3.0, 4.0 }
   };
   Matrice m2x2b = {
      { 1.0, 1.0 },
      { 2.0, 2.0 }
   };
   cout << m2x2a << '*' << endl << m2x2b << '=' << endl << m2x2a * m2x2b << endl
        << m2x2b << '*' << endl << m2x2a << '=' << endl << m2x2b * m2x2a << endl;
}

À l'exécution de ce programme, on obtient :

1 1 1 1
1 1 1 1
1 1 1 1

2.5 2.5 2.5 2.5
2.5 2.5 2.5 2.5
2.5 2.5 2.5 2.5

1 0 0 0
0 1 0 0
0 0 1 0
0 0 0 1

2.5 2.5 2.5 2.5
2.5 2.5 2.5 2.5
2.5 2.5 2.5 2.5

1 2
3 4
*
1 1
2 2
=
5 5
11 11

1 1
2 2
*
1 2
3 4
=
4 6
8 12

Ceci donne un aperçu de ce qu'il est possible de faire dans un langage comme dans l'autre avec un effort minimal.

Considérations de design

Passer d'une implémentation naïve mais utilisable, tel que demandé ici, à une implémentation de qualité plus près de celle d'un produit utilisable en pratique demande un effort beaucoup plus important... Effort que nous n'offririons pas, faute d'espace et de temps. Toutefois, quelques pistes d'exploration sont possibles, même avec des étudiant(e)s de niveau collégial, et peuvent servir de complément de formation intéressants.

Robustesse et respect des invariants

Il serait possible, dans un cas comme dans l'autre, d'assurer un meilleur respect des invariants des classes Matrice proposées. Pour ce faire, il faudrait :

Ces étapes accroîtraient la robustesse de manière superficielle, mais moins qu'on ne pourrait le croire de prime abord.

Immuabilité et aliasing

Un risque très concret avec une Matrice tient à l'aliasing, ou du partage indésiré de données. La solution au problème de l'initialisation d'une Matrice telle que celle que j'avais proposée à Nicolas est un cas patent d'un tel risque. Pour les besoins de l'illustration, en voici une version très légèrement modifiée :

// ...
   class Matrice
   {
      double[,] elems;
      public Matrice(double[,] zeElems)
      {
         elems = zeElems;
      }
      // ...
   }
   class Program
   {
      static void Main(string[] args)
      {
         var données = new double[,] {
            { 1.0, 2.0, 3.0 },
            { 4.0, 5.0, 6.0 },
            { 7.0, 8.0, 9.0 }
         };
         Matrice mat = new Matrice(données);
         // ...
      }
   }
// ...

Après l'initialisation de la Matrice nommée mat ici, le programme principal possède une variable données qui pointe vers un double [ , ] alors que mat.elems pointe sur le même double [ , ]. Ceci signifie qu'un changement à un élément dans l'un impactera l'autre à son insu, ce qui pourrait en soi poser problème si mat, par exemple, avait des invariants à faire respecter sur ses éléments (ce n'est pas le cas ici) : même si mat avait assuré le respect de ses invariants à la construction, il demeurerait vulnérable à un changement hors de son contrôle par la suite.

Un problème plus grave est le risque d'accès concurrents si le programme utilise plus d'un thread, ce qui sera par exemple le cas si une interface personne/ machine graphique utilise une Matrice comme source de données affichables. Une telle situation est susceptible d'entraîner des bogues très vilains, et très ardus à dépister.

L'aliasing est un irritant incontournable des langages qui n'utilisent les objets que de manière indirecte, ce qui inclut à ma connaissance tous les langages à collecte d'ordures (Java et les langages .NET en tête), et est susceptible de frapper tout programme C++ utilisant des pointeurs bruts ou des std::shared_ptr de manière imprudente. Il est possible de ne pas s'en préoccuper pour des programmes destinés à usage interne seulement; il est dangereux de ne pas s'en préoccuper dans le cas d'une API destinée à être utilisée par des tiers ou dans le cas de l'interface d'une classe susceptible d'être utilisée dans plusieurs projets.

En fait, c'est inexact; vous remarquerez peut-être que l'implémentation vanille::Matrice retourne, de par ses opérateurs [ ] en version const et non-const, une référence sur un vecteur plutôt qu'une copie. Dans la version non-const, c'est la chose à faire, mais dans la version const on parle d'une bêtre optimisation, qui accélère beaucoup l'exécution du programme mais pourrait éventuellement nous mordre...

En C++, privilégier la copie (et le mouvement) règle de facto ce problème. En C#, par exemple, une classe a deux grands axes à privilégier pour éviter des bogues d'aliasing :

Les interfaces immuables sont l'une des clés de la saine programmation dans les cas multiprogrammés. La programmation défensive est coûteuse (on remplace la copie d'une référence par une allocation de mémoire suivie de répétitives réalisant des copies, donc par des opérations plus coûteuses par quelques ordres de grandeur) mais est nécessaire dès qu'on fait plus que des programmes jouets. Il est, en pratique, très difficile de sécuriser en Java ou en C# une classe dont l'iinterface expose des objets qui ne sont pas immuables, du fait qu'en passant systématiquement par des références sur des objets plutôt que par des objets à part entière, le mode d'opération par défaut est le partage des pointés. Ainsi, pratiques recommandables :

Validation statique des dimensions

Bien que la syntaxe puisse être rébarbative pour qui s'y frotte une première fois, il est possible d'exprimer le code de Matrice de manière générique, flexible et à certains égards plus simple en utilisant un type Matrice dont les dimensions M et N sont des constantes connues à la compilation. évidemment, ce ne sont pas tous les cas d'application qui se prêtent à cette pratique, mais plusieurs programmes de nature scientifique opèrent sur des matrices dont les dimensions sont connues a priori.

Ce qui suit est donc une illustration, pensée pour votre divertissement et, qui sait, peut-être soulever quelques questions. Cette illustration ne permettra pas de mettre en valeur le recours à des espaces nommés, cependant, l'interface étant différente entre une classe générique et une qui ne l'est pas – on pourrait aplanir ces différences, évidemment.

Les ajustements aux bibliothèques incluses sont mineurs : j'utiliserai dans ce cas-ci des std::array<T,N> plutôt que des std::vector<T> à titre de substrat, le nombre d'éléments étant chaque fois connu a priori.

Le code de vanille::Matrice, étant logé dans un namespace, n'interférera pas avec celui de notre nouvelle classe.

#include <vector>
#include <cassert>
#include <istream>
#include <ostream>
#include <initializer_list>
#include <algorithm>
#include <array>
using namespace std;

namespace vanille
{
   // ...
}

La nouvelle classe sera statique::Matrice, un type générique sur la base d'un entier M (la hauteur), d'un entier N (la largeur) et d'un type T (le type des valeurs). Nous utiliserons double par défaut à titre de type T par conformité avec les exemples précédents.

Le constructeur par défaut remplit la Matrice avec le zéro du type. pour y arriver, cet exemple crée une ligne avec des zéros, puis copie cette ligne dans chacune des lignes de l'objet en cours de construction (ce serait insuffisant si nous utilisions un langage où, par défaut, les objets sont partagés plutôt que copiés).

Le constructeur avec une initializer_list<initializer_list<value_type>> prend soin de valider la hauteur de même que la largeur de chaque ligne, ce qui explique sa complexité.

Par contre, les accesseurs nb_lignes() et nb_colonnes() sont très simples... Ils pourraient même être qualifiés constexpr, ce qui serait considérablement avantageux.

namespace statique
{
   template <int M, int N, class T = double>
      class Matrice
      {
      public:
         using value_type = T;
      private:
         array<array<value_type, N>, M> vals_;
      public:
         using size_type = typename
            array<value_type, N>::size_type;
         Matrice()
         {
            array<value_type, N> temp;
            fill(begin(temp), end(temp), value_type());
            fill(begin(vals_), end(vals_), temp);
         }
         Matrice(initializer_list<initializer_list<value_type>> init)
         {
            assert(init.size() == M);
            size_type i = 0;
            for (auto & ligne : init)
            {
               assert(ligne.size() == N);
               size_type j = 0;
               for (auto & colonne : ligne)
               {
                  vals_[i][j] = colonne;
                  ++j;
               }
               ++i;
            }
         }
         size_type nb_lignes() const
            { return M; }
         size_type nb_colonnes() const
            { return N; }
         array<value_type, N> &operator [](size_type n)
            { return vals_[n]; }
         array<value_type, N> operator [](size_type n) const
            { return vals_[n]; }
      };

Les fonctions globales accompagnant cette classe Matrice sont aussi en partie simplifiées, tirant leur valeur de retour de la signature des types impliqués, donc de données connues dès la compilation.

// ...
      template <int M, int N, class T>
         int hauteur(const Matrice<M, N, T> &mat)
            { return M; }
      template<int M, int N, class T>
         int largeur(const Matrice<M, N, T> &mat)
            { return N; }
      template <int M, int N, class T>
         ostream &operator<<(ostream &os, const Matrice<M, N, T> &mat)
         {
            for (int i = 0; i < M; ++i)
            {
               for (int j = 0; j < N; ++j)
                  os << mat[i][j] << ' ';
               os << endl;
            }
            return os;
         }
}

Tel que mentionné plus haut, une instruction using nous donnera accès à cette nouvelle classe de même qu'à ses outils, même si les avantages de cette pratique ne transparaissent pas beaucoup ici.

// ...
using statique::Matrice;

La multiplication de matrices est légèrement simplifiée, les dimensions étant inscrites à même le type plutôt que décrites dynamiquement. Conséquemment, si le nombre de colonnes de la première ne correspond pas au nombre de lignes de la seconde, le code ne compilera tout simplement pas.

// ...
template <int M, int N, int K, class T>
   Matrice<M,N,T> operator*(const Matrice<M,K,T> &m0, const Matrice<K,N,T> &m1)
   {
      Matrice<M, N, T> resultat;
      for (int i = 0; i < M; ++i)
         for (int j = 0; j < N; ++j)
            for (int k = 0; k < K; ++k)
               resultat[i][j] += m0[i][k] * m1[k][j];
      return resultat;
   }
template <int M, int N, int K, class T>
   Matrice<M, N, T> produit_matriciel(const Matrice<M, K, T> &m0, const Matrice<K, N, T> &m1)
      { return m0 * m1; }

La multiplication d'une matrice par un scalaire est légèrement simplifiée, sans plus.

// ...
template <int M, int N, class T>
   Matrice<M,N,T> operator*(Matrice<M,N,T> mat, const T &val)
   {
      for (int i = 0; i < M; ++i)
         for (int j = 0; j < N; ++j)
            mat[i][j] *= val;
      return mat;
   }
template <int M, int N, class T>
   Matrice<M,N,T> operator*(const T &val, Matrice<M,N,T> mat)
      { return mat * val; }

Le prédicat est_carree() devient une constante, susceptible d'être constexpr, du fait que l'état d'être carrée ou non est inscrit à même le type de la matrice.

// ...
template <int M, int N, class T>
   bool est_carree(const Matrice<M,N,T> &mat)
      { return false; }
template <int N, class T>
   bool est_carree(const Matrice<N, N, T> &mat)
      { return true; }

Une matrice identité d'une certaine taille est maintenant construite à partir d'une taille statique. Le code en est encore une fois allégé.

// ...
template <int N, class T = double>
   Matrice<N,N,T> creer_matrice_identite()
   {
      Matrice<N,N,T> mat;
      for (int i = 0; i < N; ++i)
         mat[i][i] = 1.0;
      return mat;
   }

Enfin, le code client n'a que peu changé, outre la manière (statique plutôt que dynamique) d'instancier les matrices.

// ...
#include <iostream>
int main()
{
   Matrice<3,4> mat = {
      { 1.0, 1.0, 1.0, 1.0 },
      { 1.0, 1.0, 1.0, 1.0 },
      { 1.0, 1.0, 1.0, 1.0 }
   };
   cout << mat << endl;
   mat = mat * 2.5;
   cout << mat << endl;
   auto mid4x4 = creer_matrice_identite<4>();
   cout << mid4x4 << endl;
   //
   // boum!
   // cout << mid4x4 * mat << endl;
   //
   cout << mat * mid4x4 << endl;
   Matrice<2,2> m2x2a = {
      { 1.0, 2.0 },
      { 3.0, 4.0 }
   };
   Matrice<2,2> m2x2b = {
      { 1.0, 1.0 },
      { 2.0, 2.0 }
   };
   cout << m2x2a << '*' << endl << m2x2b << '=' << endl << m2x2a * m2x2b << endl
        << m2x2b << '*' << endl << m2x2a << '=' << endl << m2x2b * m2x2a << endl;
}

En espérant que cet exemple, qui ne fait qu'effleurer la surface de ce qui est possible, vous donne des idées. Vous remarquerez que dans le cas du code client, il n'y a que peu de changements, et ceux qui y apparaissent ne sont pas (à mon humble avis) déraisonnables.

Organisation des données et optimisation

Enfin, jetons un coup d'oeil à des questions d'organisation des données et d'optimisation.

Accès à la Cache

La version naïve de l'opérateur de multiplication de deux matrices (version C++, dimensions statiques pour fins de simplicité) est répétée à droite. Cette implémentation est terriblement inefficace, pour des raisons qui vous auront peut-être échappé, et qui ont trait au matériel contemporain. En effet, supposant (pour simplifier le propos) que T soit le type double :

  • la matrice m0 est un tableau de double. Ainsi, dans m0, les données de chaque ligne vont de à inclusivement, et sont contiguës en mémoire pour chaque ligne;
  • la matrice m1 est un tableau de double. Ainsi, dans m1, les données de chaque ligne vont de à inclusivement, et sont contiguës en mémoire pour chaque ligne;
  • la matrice resultat est un tableau de double. Ainsi, dans resultat, les données de chaque ligne vont de à inclusivement, et sont contiguës en mémoire pour chaque ligne;
  • dans un ordinateur contemporain, les accès à la mémoire dominent le temps d'exécution des programmes. Il est donc essentiel de choisir des structures de données qui favorisent un accès efficace aux données. Concrètement, cela signifie des structures de données contiguës en mémoire (tableaux, vecteurs);
  • il est aussi essentiel d'utiliser des algorithmes qui profitent au maximum de cette proximité. C'est là que nous pouvons faire mieux.
// ...
template <int M, int N, int K, class T>
   Matrice<M,N,T> operator*(const Matrice<M,K,T> &m0, const Matrice <K,N,T>&m1)
{
   Matrice<M,N,T> resultat;
   for (int i = 0; i < M; ++i)
      for (int j = 0; j < N; ++j)
         for (int k = 0; k < K; ++k)
            resultat[i][j] += m0[i][k] * m1[k][j];
   return resultat;
}

En effet, dans l'expression resultat[i][j] += m0[i][k] * m1[k][j], les données indicées par j varient pour deux des trois matrices sur une plage contiguë en mémoire, alors que la répétitive interne (la Tight Loop) fait varier k, pas j, ce qui fait que dans deux des trois matrices, les accès successifs ne sont pas faits sur des éléments adjacents en mémoire. Cette faible localité des données a un effet adverse sur la qualité de la gestion de la Cache, et fait dégénérer la vitesse d'exécution.

Le simple fait d'intervertir les deux boucles les plus internes de l'algorithme, sans changer quoi que ce soit aux résultats du calcul, accroît très fortement la localité des accès dans la boucle la plus interne. Ce faisant, l'algorithme profite nettement plus de la Cache, et la vitesse d'exécution du programme résultant est elle aussi accrue... De manière significative!

// ...
template <int M, int N, int K, class T>
   Matrice<M,N,T> operator*(const Matrice<M,K,T> &m0, const Matrice <K,N,T>&m1)
{
   Matrice<M,N,T> resultat;
   for (int i = 0; i < M; ++i)
      for (int k = 0; k < K; ++k)
         for (int j = 0; j < N; ++j)
            resultat[i][j] += m0[i][k] * m1[k][j];
   return resultat;
}

À titre d'exemple, comparons les deux algorithmes à partir du programme de test suivant :

// ...
#include <chrono>
using namespace std::chrono;
int main()
{
   enum { N = 100, NTESTS = 1000 };
   Matrice<N, N> m0, m1;
   auto avant = system_clock::now();
   for (int i = 0; i < NTESTS; ++i)
      m0 = m0 * m1;
   auto apres = system_clock::now();
   cout << "mat " << N << 'x' << N << " * mat " << N << 'x' << N << ", " << NTESTS << " fois en "
        << duration_cast<milliseconds>(apres - avant).count() << " ms." << endl;
}

Un algorithme de complexité comme la multplication matricielle peut prendre un temps significatif à exécuter, même pour une paire de matrices . Avec ce même programme de test, en utilisant l'algoirithme naïf (boucles sur i, j et k dans l'ordre) et l'algorithme restructuré pour profiter de la Cache (boucles sur i, k et j dans l'ordre), j'obtiens sur un même ordinateur les résultats suivants :

Algorithme naïf (i, j, k)Algorithme Cache-Friendly (i, k, j)
mat 100x100 * mat 100x100, 1000 fois en 38002 ms.
mat 100x100 * mat 100x100, 1000 fois en 1166 ms.

Vous ne rêvez pas : l'algorithme Cache-Friendly prend environ (autour de ) du temps d'exécution de l'algorithme naïf pour les mêmes opérations sur les mêmes données.

De manière amusante, avec l'implémentation C# proposée par Nicolas, le compilateur réorganise (en mode Release) les boucles dans chaque cas pour offrir des performances à peu près identiques dans chaque cas, soit environ huit fois le temps de l'algorithme Cache-Friendly en C++ (8212 ms). pour une version cumulant dans une variable temporaire, comme le fait habilement Nicolas dans sa proposition, et environ dix fois le temps de la version Cache-Friendly (10271 ms.) si le cumul est fait directement dans la matrice résultante du calcul avec l'accesseur [ , ] – en C#, les propriétés sont sympathiques sur le plan syntaxique, mais elles coûtent cher! Il y a plus que l'accès à la propriété ici, cela dit : opérer sur une donnée située sur la pile permet une meilleure localité des accès en mémoire que ne le fait opérer sur une donnée ailleurs en mémoire, et cette meilleure localité entraîne un gain à chaque écriture. C'est significatif.

Le modèle de compilation et d'exécution de C# est justement pensé pour des programmes où l'exécution passe fréquemment par les mêmes endroits dans le code, compilant les fonctions en code machine dans une antémémoire et offrant des performances très décentes dans un cas comme celui-ci.

Profiter de la sémantique de mouvement

Dans le cas de l'implémentation par vecteurs en C++, quelques-unes des fonctions (dont l'opérateur de multiplication) tendent à retourner des variables par copie. En implémentant la sémantique de mouvement (chose très simple dans ce cas-ci), les copies de vecteurs de éléments sont remplacés par des déplacements de pointeurs, une opération en temps constant (complexité ). Contrairement aux agrégats simples que sont les tableaux bruts ou les std::array, un std::vector alloue de la mémoire dynamiquement, ce qui fait dd sa copie une opération potentiellement coûteuse.

Les retouches à faire dans la classe vanille::Matrice sont les suivantes :

  • ajouter un constructeur de mouvement, qui sera utilisé par le compilateur lorsque l'objet sources est manifestement destiné à être éliminé de toute manière; et
  • ajouter un opérateur d'affectation par mouvement, qui sera utilisé lui aussi lorsque la source est un objet jetable.

Dans notre programme, les performances à l'exécution sont bien moindres que ce qu'elles pourraient être avec une implémentation naïve (qui copie plutôt que déplacer) puisque nous avons favorisé l'opérateur * binaire (à deux opérandes), qui crée chaque fois une variable temporaire pour contenir le fruit de la multiplication. En fait, si nous en restons à une implémentation naïve, le code C# sera plus rapide que le code C++, du fait qu'on y manipule les objets par des indirections, évitant par le fait-même les copies.

// ...
       Matrice(Matrice && m)
         : vals_(std::move(m.vals_))
      {
      }
      Matrice & operator=(Matrice && m)
      {
         vals_ = std::move(m.vals_);
         return *this;
      }
// ...

Nous testerons notre implémentation avec le programme suivant :

// ...
#include <chrono>
using namespace std::chrono;
int main()
{
   enum { N = 100, NTESTS = 1000 };
   Matrice m0{ N, N }, m1{ N, N };
   auto avant = system_clock::now();
   for (int i = 0; i < NTESTS; ++i)
      m0 = m0 * m1;
   auto apres = system_clock::now();
   cout << "mat " << N << 'x' << N << " * mat " << N << 'x' << N << ", " << NTESTS << " fois en "
      << duration_cast<milliseconds>(apres - avant).count() << " ms." << endl;
}

Le tableau suivant montre un comparatif des vitesses avec un algorithme naïf ou Cache-Friendly, avec ou sans support pour la sémantique de mouvement, de même qu'un test utilisant des variables temporaires pour l'évaluation des dimensions des matrices (dans le but d'offrir une meilleure localité) :

Algorithme naïf (i, j, k), sans mouvementAlgorithme Cache-Friendly (i, k, j), sans mouvement
mat 100x100 * mat 100x100, 1000 fois en 2642 ms.
mat 100x100 * mat 100x100, 1000 fois en 1420 ms.
Algorithme naïf (i, j, k), avec mouvementAlgorithme Cache-Friendly (i, k, j), avec mouvement
mat 100x100 * mat 100x100, 1000 fois en 2623 ms.
mat 100x100 * mat 100x100, 1000 fois en 1353 ms.
Algorithme naïf (i, j, k), avec mouvement et temporairesAlgorithme Cache-Friendly (i, k, j), avec mouvement et temporaires
mat 100x100 * mat 100x100, 1000 fois en 1267 ms.
mat 100x100 * mat 100x100, 1000 fois en 671 ms.

On comprendra que de la version la plus sophistiquée du lot consomme environ () de la version la plus naïve d'entre elles. Le jeu en vaut la chandelle.


Valid XHTML 1.0 Transitional

CSS Valide !