Voici les codes des algorithmes. Le fichier complet est MATRICE.CPP.
Les codes ont été vérifiés avec Purify de Rational.

#define TYPE double // Type des éléments des matrices
#define N 720 // Taille des matrices
#define M 720
#define P 720

TYPE m1[N][N]; // Matrices: r = m1 . m2
TYPE m2[N][N];
TYPE r[N][N];

TYPE a[N][M];
TYPE b[M][P];
TYPE p[N][P];

//-------------------------------------------------------------------------------------------------
// r = 0
// p = 0

void clear(void)
{
 int
  l,
  c;

 for(l = 0; l < N; l++)
  for(c = 0; c < N; c++)
   r[l][c] = 0.0;
 for(l = 0; l < N; l++)
  for(c = 0; c < P; c++)
   p[l][c] = 0.0;
}

//-------------------------------------------------------------------------------------------------
// Vérifie que r = m1 . m2
int testr(void)
{
 int
  l,
  c,
  k,
  t = 1;
 double
  v;

 for(l = 0; l < N; l++)
  for(c = 0; c < N; c++)
  {
   v = 0.0;
   for(k = 0; k < N; k++)
    v += m1[l][k] * m2[k][c];
   t &= r[l][c] == v;
  }

 if(t)
  cout << "Ok" << " - ";
 else
  cout << "Erreur" << " - ";
 return t;
}

//-------------------------------------------------------------------------------------------------
// Vérifie que p = a . b
int testp(void)
{
 int
  l,
  c,
  k,
  t = 1;
 double
  v;

 for(l = 0; l < N; l++)
  for(c = 0; c < P; c++)
  {
   v = 0.0;
   for(k = 0; k < M; k++)
    v += a[l][k] * b[k][c];
   t &= p[l][c] == v;
  }

 if(t)
  cout << "Ok" << " - ";
 else
  cout << "Erreur" << " - ";
 return t;
}

//-------------------------------------------------------------------------------------------------
// Algorithme de base

void mul_1(void)
{
 int
  l,
  c,
  i;

 for(l = 0; l < N; l++)
  for(c = 0; c < N; c++)
  {
   r[l][c] = 0.0;
   for(i = 0; i < N; i++)
    r[l][c] += m1[l][i] * m2[i][c];
  }
}

//-------------------------------------------------------------------------------------------------
// Algorithme de base avec variable intermédiaire

void mul_2(void)
{
 int
  l,
  c,
  i;
 TYPE s;

 for(l = 0; l < N; l++)
  for(c = 0; c < N; c++)
  {
   s = 0.0;
   for(i = 0; i < N; i++)
    s += m1[l][i] * m2[i][c];
   r[l][c] = s;
  }
}

//-------------------------------------------------------------------------------------------------
// Algorithme de base avec initialisation

void mul_3(void)
{
 int
  l,
  c,
  i;

 for(l = 0; l < N; l++)
  for(c = 0; c < N; c++)
  {
   r[l][c] = m1[l][0] * m2[0][c];
   for(i = 1; i < N; i++)
    r[l][c] += m1[l][i] * m2[i][c];
  }
}

//-------------------------------------------------------------------------------------------------
// Algorithme de base avec variable intermédiaire et initialisation

void mul_4(void)
{
 int
  l,
  c,
  i;
 TYPE s;

 for(l = 0; l < N; l++)
  for(c = 0; c < N; c++)
  {
   s = m1[l][0] * m2[0][c];
   for(i = 1; i < N; i++)
    s += m1[l][i] * m2[i][c];
   r[l][c] = s;
  }
}

//-------------------------------------------------------------------------------------------------
// Maeno Toshinori
// Tokyo Institute of Technology
// N mod L = 0

#define L 30 // Taille des matrices placées en cache
#define L2 30

void mul_5(void)
{
 int
  i,
  j,
  k,
  kk,
  i2,
  kt,
  it;
 double
  t0,
  t1,
  t2,
  t3,
  s,
  t4,
  t5,
  t6,
  t7;

 for(i = 0; i < N; i++)
  for(j = 0; j < N; j++)
   r[i][j] = 0.0;
 for(i2 = 0; i2 < N; i2 += L2)
  for(kk = 0; kk < N; kk += L)
  {
   it = i2 + L;
   kt = kk + L;
   for(j = 0; j < N; j += 4)
    for(i = i2; i < it; i += 2)
    {
     t0 = t1 = t2 = t3 = 0.0;
     t4 = t5 = t6 = t7 = 0.0;
     for(k = kk; k < kt; k++)
     {
      s = m1[i][k];
      t0 += s * m2[k][j];
      t1 += s * m2[k][j + 1];
      t2 += s * m2[k][j + 2];
      t3 += s * m2[k][j + 3];
      s = m1[i + 1][k];
      t4 += s * m2[k][j];
      t5 += s * m2[k][j + 1];
      t6 += s * m2[k][j + 2];
      t7 += s * m2[k][j + 3];
     }
     r[i][j] += t0;
     r[i][j + 1] += t1;
     r[i][j + 2] += t2;
     r[i][j + 3] += t3;
     r[i + 1][j] += t4;
     r[i + 1][j + 1] += t5;
     r[i + 1][j + 2] += t6;
     r[i + 1][j + 3] += t7;
    }
  }
}

//-------------------------------------------------------------------------------------------------
// Dan Warner
// Dept. of Mathematics, Clemson University

#define NB 40

void mul_6(void)
{
 int
  i,
  j,
  k,
  ii,
  jj,
  kk;
 double
  s00,
  s01,
  s10,
  s11;

 for(ii = 0; ii < N; ii += NB)
  for(jj = 0; jj < N; jj += NB)
  {
   for(i = ii; i < ii + NB; i++)
    for(j = jj; j < jj + NB; j++)
     r[i][j] = 0.0;
   for(kk = 0; kk < N; kk += NB)
    for(i = ii; i < ii + NB; i += 2)
     for(j = jj; j < jj + NB; j += 2)
     {
      s00 = r[i][j];
      s01 = r[i][j + 1];
      s10 = r[i + 1][j];
      s11 = r[i + 1][j + 1];
      for(k = kk; k < kk + NB; k++)
      {
       s00 += m1[i][k] * m2[k][j];
       s01 += m1[i][k] * m2[k][j + 1];
       s10 += m1[i + 1][k] * m2[k][j];
       s11 += m1[i + 1][k] * m2[k][j + 1];
      }
      r[i][j] = s00;
      r[i][j + 1] = s01;
      r[i + 1][j] = s10;
      r[i + 1][j + 1] = s11;
     }
  }
}

//-------------------------------------------------------------------------------------------------
// Inversion des boucles les plus intérieures

void mul_7(void)
{
 int
  i,
  j,
  k;

 for(i = 0; i < N; i++)
  for(j = 0; j < N; j++)
   r[i][j] = 0.0;
 for(i = 0; i < N; i++)
  for(k = 0; k < N; k++)
   for(j = 0; j < N; j++)
    r[i][j] += m1[i][k] * m2[k][j];
}

//-------------------------------------------------------------------------------------------------
// Inversion des boucles les plus intérieures avec variable intermédiaire

void mul_8(void)
{
 int
  i,
  j,
  k;
 double
  temp;

 for(i = 0; i < N; i++)
  for(j = 0; j < N; j++)
   r[i][j] = 0.0;
 for(i = 0; i < N; i++)
  for(k = 0; k < N; k++)
  {
   temp = m1[i][k];
   for(j = 0; j < N; j++)
    r[i][j] += temp * m2[k][j];
  }
}

//-------------------------------------------------------------------------------------------------
// Pavage de longueur step

#define min(a, b) ((a) < (b) ? (a) : (b))

void mul_9(void)
{
 int
  i,
  j,
  k,
  jj,
  kk,
  step = 8; // 4 or 8 or 16 or 32 or 64 or 96
 double
  temp;

 for(i = 0; i < N; i++)
  for(j = 0; j < N; j++)
   r[i][j] = 0.0;
 for(kk = 0; kk < N; kk += step)
  for(jj = 0; jj < N; jj += step)
   for(i = 0; i < N; i++)
    for(k = kk; k < min(kk + step, N); k++)
    {
     temp = m1[i][k];
     for(j = jj; j < min(jj + step, N); j++)
      r[i][j] += temp * m2[k][j];
    }

}

//-------------------------------------------------------------------------------------------------
// Boucle raccourcie de 4
// N >= 4

void mul_10(void)
{
 int
  i,
  j,
  k;
 double
  temp;

 for(i = 0; i < N; i++)
  for(j = 0; j < N; j++)
  {
   temp = 0.0;
   for(k = 0; k < N - 3; k += 4)
   {
    temp += m1[i][k] * m2[k][j];
    temp += m1[i][k + 1] * m2[k + 1][j];
    temp += m1[i][k + 2] * m2[k + 2][j];
    temp += m1[i][k + 3] * m2[k + 3][j];
   }
   for(; k < N; k++)
    temp += m1[i][k] * m2[k][j];
   r[i][j] = temp;
  }
}

//-------------------------------------------------------------------------------------------------
// Boucle raccourcie de 8
// N >= 8

void mul_11(void)
{
 int
  i,
  j,
  k;
 double
  temp;

 for(i = 0; i < N; i++)
  for(j = 0; j < N; j++)
  {
   temp = 0.0;
   for(k = 0; k < N - 7; k += 8)
   {
    temp += m1[i][k] * m2[k][j];
    temp += m1[i][k + 1] * m2[k + 1][j];
    temp += m1[i][k + 2] * m2[k + 2][j];
    temp += m1[i][k + 3] * m2[k + 3][j];
    temp += m1[i][k + 4] * m2[k + 4][j];
    temp += m1[i][k + 5] * m2[k + 5][j];
    temp += m1[i][k + 6] * m2[k + 6][j];
    temp += m1[i][k + 7] * m2[k + 7][j];
   }
   for(; k < N; k++)
    temp += m1[i][k] * m2[k][j];
   r[i][j] = temp;
  }
}

//-------------------------------------------------------------------------------------------------
// Boucle raccourcie de 16
// N >= 16

void mul_12(void)
{
 int
  i,
  j,
  k;
 double
  temp;

 for(i = 0; i < N; i++)
  for(j = 0; j < N; j++)
  {
   temp = 0.0;
   for(k = 0; k < N - 15; k += 16)
   {
    temp += m1[i][k] * m2[k][j];
    temp += m1[i][k + 1] * m2[k + 1][j];
    temp += m1[i][k + 2] * m2[k + 2][j];
    temp += m1[i][k + 3] * m2[k + 3][j];
    temp += m1[i][k + 4] * m2[k + 4][j];
    temp += m1[i][k + 5] * m2[k + 5][j];
    temp += m1[i][k + 6] * m2[k + 6][j];
    temp += m1[i][k + 7] * m2[k + 7][j];
    temp += m1[i][k + 8] * m2[k + 8][j];
    temp += m1[i][k + 9] * m2[k + 9][j];
    temp += m1[i][k + 10] * m2[k + 10][j];
    temp += m1[i][k + 11] * m2[k + 11][j];
    temp += m1[i][k + 12] * m2[k + 12][j];
    temp += m1[i][k + 13] * m2[k + 13][j];
    temp += m1[i][k + 14] * m2[k + 14][j];
    temp += m1[i][k + 15] * m2[k + 15][j];
   }
   for(; k < N; k++)
    temp += m1[i][k] * m2[k][j];
   r[i][j] = temp;
  }
}

//-------------------------------------------------------------------------------------------------
// Utilisation de pointeurs
// N mod 8 = 0

void mul_13(void)
{
 int
  i,
  j,
  k;
 double
  t,
  *aa,
  *bb;

 for(i = 0; i < N; i++)
  for(j = 0; j < N; j++)
  {
   t = 0;
   aa = m1[i];
   bb = &m2[0][j];
   for(k = 0; k < N / 8; k++)
   {
    t += *aa++ * bb[0];
    t += *aa++ * bb[N];
    t += *aa++ * bb[2 * N];
    t += *aa++ * bb[3 * N];
    t += *aa++ * bb[4 * N];
    t += *aa++ * bb[5 * N];
    t += *aa++ * bb[6 * N];
    t += *aa++ * bb[7 * N];
    bb += 8 * N;
   }
   r[i][j] = t;
  }
}

//-------------------------------------------------------------------------------------------------
// Vincent LESAUVAGE
// N mod 2 = 0
// Carrés de 2

void mul_14(void)
{
 int
  c,
  c2,
  l;
 double
  t0,
  t1,
  t2,
  t3;

 for(l = 0; l < N; l++)
  for(c = 0; c < N; c++)
   r[l][c] = 0.0;
 for(l = 0; l < N; l += 2)
  for(c = 0; c < N; c += 2)
  {
   t0 = t1 = t2 = t3 = 0.0;
   for(c2 = 0; c2 < N; c2 += 2)
   {
    t0 += m1[l][c2] * m2[c2][c] + m1[l][c2 + 1] * m2[c2 + 1][c];
    t1 += m1[l][c2] * m2[c2][c + 1] + m1[l][c2 + 1] * m2[c2 + 1][c + 1];
    t2 += m1[l + 1][c2] * m2[c2][c] + m1[l + 1][c2 + 1] * m2[c2 + 1][c];
    t3 += m1[l + 1][c2] * m2[c2][c + 1] + m1[l + 1][c2 + 1] * m2[c2 + 1][c + 1];
   }
   r[l][c] = t0;
   r[l][c + 1] = t1;
   r[l + 1][c] = t2;
   r[l + 1][c + 1] = t3;
  }
}

//-------------------------------------------------------------------------------------------------
// Vincent LESAUVAGE
// N mod 3 = 0
// Carrés de 3

void mul_15(void)
{
 int
  c,
  c2,
  l;
 double
  t0,
  t1,
  t2,
  t3,
  t4,
  t5,
  t6,
  t7,
  t8;

 for(l = 0; l < N; l++)
  for(c = 0; c < N; c++)
   r[l][c] = 0.0;
 for(l = 0; l < N; l += 3)
  for(c = 0; c < N; c += 3)
  {
   t0 = t1 = t2 = 0.0;
   t3 = t4 = t5 = 0.0;
   t6 = t7 = t8 = 0.0;
   for(c2 = 0; c2 < N; c2 += 3)
   {
    t0 += m1[l][c2] * m2[c2][c] + m1[l][c2 + 1] * m2[c2 + 1][c] +
     m1[l][c2 + 2] * m2[c2 + 2][c];
    t1 += m1[l][c2] * m2[c2][c + 1] + m1[l][c2 + 1] * m2[c2 + 1][c + 1] +
     m1[l][c2 + 2] * m2[c2 + 2][c + 1];
    t2 += m1[l][c2] * m2[c2][c + 2] + m1[l][c2 + 1] * m2[c2 + 1][c + 2] +
     m1[l][c2 + 2] * m2[c2 + 2][c + 2];
 
    t3 += m1[l + 1][c2] * m2[c2][c] + m1[l + 1][c2 + 1] * m2[c2 + 1][c] +
     m1[l + 1][c2 + 2] * m2[c2 + 2][c];
    t4 += m1[l + 1][c2] * m2[c2][c + 1] + m1[l + 1][c2 + 1] * m2[c2 + 1][c + 1] +
     m1[l + 1][c2 + 2] * m2[c2 + 2][c + 1];
    t5 += m1[l + 1][c2] * m2[c2][c + 2] + m1[l + 1][c2 + 1] * m2[c2 + 1][c + 2] +
     m1[l + 1][c2 + 2] * m2[c2 + 2][c + 2];

    t6 += m1[l + 2][c2] * m2[c2][c] + m1[l + 2][c2 + 1] * m2[c2 + 1][c] +
     m1[l + 2][c2 + 2] * m2[c2 + 2][c];
    t7 += m1[l + 2][c2] * m2[c2][c + 1] + m1[l + 2][c2 + 1] * m2[c2 + 1][c + 1] +
     m1[l + 2][c2 + 2] * m2[c2 + 2][c + 1];
    t8 += m1[l + 2][c2] * m2[c2][c + 2] + m1[l + 2][c2 + 1] * m2[c2 + 1][c + 2] +
     m1[l + 2][c2 + 2] * m2[c2 + 2][c + 2];
   }
   r[l][c] = t0;
   r[l][c + 1] = t1;
   r[l][c + 2] = t2;
   r[l + 1][c] = t3;
   r[l + 1][c + 1] = t4;
   r[l + 1][c + 2] = t5;
   r[l + 2][c] = t6;
   r[l + 2][c + 1] = t7;
   r[l + 2][c + 2] = t8;
  }
}

//-------------------------------------------------------------------------------------------------
// Vincent LESAUVAGE
// a (NxM), b (MxP), p (NxP)
// N, M, P mod 2 = 0
// Carrés de 2

void mul_16(void)
{
 int
  c,
  c2,
  l;
 double
  t0,
  t1,
  t2,
  t3;

 for(l = 0; l < N; l++)
  for(c = 0; c < P; c++)
   p[l][c] = 0.0;
 for(l = 0; l < N; l += 2)
  for(c = 0; c < P; c += 2)
  {
   t0 = t1 = t2 = t3 = 0.0;
   for(c2 = 0; c2 < M; c2 += 2)
   {
    t0 += a[l][c2] * b[c2][c] + a[l][c2 + 1] * b[c2 + 1][c];
    t1 += a[l][c2] * b[c2][c + 1] + a[l][c2 + 1] * b[c2 + 1][c + 1];
    t2 += a[l + 1][c2] * b[c2][c] + a[l + 1][c2 + 1] * b[c2 + 1][c];
    t3 += a[l + 1][c2] * b[c2][c + 1] + a[l + 1][c2 + 1] * b[c2 + 1][c + 1];
   }
   p[l][c] = t0;
   p[l][c + 1] = t1;
   p[l + 1][c] = t2;
   p[l + 1][c + 1] = t3;
  }
}