Wierszami czy kolumnami? cz III. (Mnożenie macierzy)

W niniejszej notce omówię algorytm mnożenia macierzy. Jeśli czytelnik z mnożeniem macierzy spotyka się po raz pierwszy to odsyłam do wikipedii (http://pl.wikipedia.org/wiki/Mno%C5%BCenie_macierzy). Operacja mnożenia macierzy jest tam bardzo ładnie wytłumaczona.
Przed napisaniem tej notki wyszukałem w googlach hasło "mnożenie macierzy C++". Pobieżne przejrzenie wyników potwierdziło moje przypuszczenia - zdecydowana większość programistów pisze pętle w złej kolejności. Przedstawmy najpierw typową implementację mnożenia macierzy:

const int N = 1000;
static int T[N][N];
static int A[N][N];
static int B[N][N];

for (int i = 0; i < N; i++)
  for (int j = 0; j < N; j++)
  {
    T[i][j] = 0;
    for (int k = 0; k < N; k++)
      T[i][j] += A[i][k] * B[k][j];
  }

Gdzie leży błąd można na pierwszy rzut oka nie zauważyć. Wytłumaczę to zatem, aby wszyscy mogli to zauważyć, a w razie innego algorytmu - powtórzyć rozumowanie.
Z tablicą T nie ma problemów. Wewnątrz pętli odwołujemy się do elementu T[i][j]. Przy czym przechodzimy ją po wierszach. To znaczy, że pętla iterująca po j-tach znajduje się wewnątrz pętli iterującej po zmiennej i. Podobna sytuacja jest z tablicą A. Odwołujemy się do elementu A[i][k], czyli przechodzimy ją po wierszach (pętla iterująca po zmiennej k znajduje się wewnątrz pętli iterującej po zmiennej i). Problem jest z tablicą B.
W przypadku tablicy B odwołujemy się do elementu B[k][j]. Przechodzimy ją zatem po kolumnach - pętla iterująca po j-tach znajduje się na zewnątrz pętli iterującej po zmiennej k. Procesor zatem nie będzie wczytywał optymalnie tablicy B do pamięci cache.
Możemy to poprawić przestawiając miejscami pętle ze zmienną k z pętlą ze zmienną j. Musimy też stworzyć osobne pętle zerujące tablicę.

const int N = 1000;
static int T[N][N];
static int A[N][N];
static int B[N][N];

for (int i = 0; i < N; i++)
  for (int j = 0; j < N; j++)
    T[i][j] = 0;

for (int i = 0; i < N; i++)
  for (int k = 0; k < N; k++)
    for (int j = 0; j < N; j++)
      T[i][j] += A[i][k] * B[k][j];

Tutaj już wszystko powinno działać w porządku. Wszystkie tablice przechodzimy po wierszach.
Przykład ten jest często wykorzystywany przy omawianiu przechodzenia tablicy po wierszach i kolumnach. Oba algorytmy mają nawet swoje nazwy potoczne. Pierwszy z nich nazywamy algorytmem ijk, natomiast drugi (jak nietrudno się domyślić) ikj.
Czas na przetestowanie obu rozwiązań. Różnica jest odczuwalna:
algorytm ijk - 0:13.63
algorytm ikj - 0:10.45
Motto tego trzy częściowego kursu jest następujące: "Iterowanie po wierszach przyśpiesza program". Czasem o kilka sekund, czasem nawet ośmiokrotnie. Warto o tym pamiętać bo przestawianie dwóch pętli nas nic nie kosztuje.