1. Użyty algorytm (macierze jak w języku C tzn. wierszami)

for i := 1 to n do
  for k := 1 to n do begin
    A_tmp := A[i][k];
    for j := 1 to n do 
      C[i][j] := C[i][j] + A_tmp*B[k][j];
  end;

Zmienna A_tmp jest zmienną rejestrową.

Załączam pełne programy: algorytm odniesienia ijk, algorytm ikj oraz
algorytm ikj z rozwinięta pętlą wewnętrzną (4-krotnie).

2. Kompilator:
gcc version egcs-2.91.66 19990314/Linux (egcs-1.1.2 release)

Dla algorytmu odniesienia ijk była wykonana zwykła kompilacja i kompilacja
z optymalizacją (opcja -O2).

Dla obu algorytmów typu ikj kompilacja była wykonana z wyłączoną
optymalizacją (opcja -O0).

3. Czasy wykonania algorytmów. 
Macierze rzeczywiste pojedyńczej precyzji (float), globalne.
Rozmiar typu float: 4 bajty (sizeof(float) = 4).

rozmiar    ijk   ijk -O2  ikj -O0  ikj4 -O0
macierzy   [s]     [s]      [s]       [s]

256        16.37    9.68     6.31     5.87
512       130.37   78.08    50.15    47.59
1024     1098.40  636.35   417.62   383.65

rozmiar   ikj/ijk ikj/ijkO2 ikj4/ijk ikj4/ijkO2
macierzy    [%]     [%]       [%]       [%] 

512        38.47   64.24     36.50     60.95
1024       38.02   65.63     34.93     60.29 

4. Przyspieszenie (najlepszy algorytm vs. algorytm wzorcowy z
optymalizacją -O2) dla n = 1024

636.350785/383.652896 * 100% = 166 %

5. Parametry maszyny testowej
pamięć: 40960k RAM
procesor: Intel Pentium 75
zegar: 133271115 Hz = 133 MHz
wydajność: 53.04 BogoMIPS

6. System operacyjny
Linux version 2.2.5-15 (Redhat 6.0)

7. Dyskusja wyników:
Algorytm ikj ma tę przewagę nad algorytmem ijk, że dostęp do danych
następuje sekwencyjnie. Jeszcze większą wydajność można uzyskać używając
kody ikj na maszynach wektorowych.

Ulepszenie algorytmu przez rozwinięcie pętli wewnętrznej nie wpływa
znacząco na osiągi.


UWAGA 1: Poprzednie wyniki były osiągnięte na innej maszynie!

UWAGA 2: Pomiar czasu był dokonany za pomocą funkcji gettimeofday(...),
która podaje czas z dokładnością do mikrosekund. Podaje ona jednak czas
systemowy. Na systemach wielowątkowych nie musi być to dokładnie czas jaki
zajęło wykonanie zadania. Można by użyć funkcji getrusage(...) aby znaleźć
czas jaki zajęło wykonanie zadania. Funkcja times(...) która powinna to
także robić (t.j. podawac czas użytkownika, a nie systemowy) nie ma
odpowiedniej rozdzielczości (niezgodnie z opisem w info i manualach). 

UWAGA 3: Warto obejrzeć wyniki czas(rozmiar) na wykresie podwójnie
logarytmicznym. Np. za pomocą gnuplota:
gnuplot> set logscale xy
gnuplot> plot 'ijk.dat', 'ikj.dat'
Na życzenie służę danymi dla n = 32..1024 dla ijk, ijkO2, ikjO0, ikj4_O0. 

-- 
Jakub Narębski
  MISMaP UW
  
  
/******************************************
Algorytm odniesienia
*******************************************/

#include <stdio.h>
#include <time.h>
#include <sys/times.h>
#include <sys/time.h>
#include <unistd.h>

#define MAXN 1024


float A[MAXN][MAXN],B[MAXN][MAXN],C[MAXN][MAXN];


void init(int n);
int timeval_subtract( struct timeval *result, 
		      struct timeval *x, 
		      struct timeval *y );



int main(void)
{
  int i,j,k;
  int n;

  FILE *f;

  struct timeval t_beg, t_end, t_diff;
  struct tms c_beg, c_end;
  clock_t c_diff;


  if((f = fopen("matmul_ijk.dat","w")) == NULL) {
    perror("Blad otwarcia pliku\n");
    exit(1);
  }
  
  fprintf(f,"floats: %d\n",sizeof(A[0][0]));
  printf("zaczynamy: IJK %d-%d\n",16,MAXN);

  for(n = 16; n <= MAXN; n *= 2) {
   
    printf("%4i init ",n);
    
    init(n);

    printf("before ");

    gettimeofday(&t_beg, NULL);

    printf("%li.%li:",t_beg.tv_sec,t_beg.tv_usec);
    fflush(stdout);

    /* ---------------------------------------- */

    times(&c_beg);
    gettimeofday(&t_beg, NULL);
      
    for(i = 0; i < n; i++)
      for(j = 0; j < n; j++)
	for(k = 0; k < n; k++)
	  C[i][j] += A[i][k]*B[k][j];
      
    gettimeofday(&t_end, NULL);
    times(&c_end);
    
    /* ---------------------------------------- */

    timeval_subtract(&t_diff, &t_end, &t_beg);
    c_diff = c_end.tms_utime - c_beg.tms_utime;

    printf("= %li.%li\n",t_diff.tv_sec,t_diff.tv_usec);

    fprintf(f,"%i %li.%li %li %.20g\n",
	    n, t_diff.tv_sec, t_diff.tv_usec, c_diff,
	    ((double)c_diff)/CLOCKS_PER_SEC);
    fflush(f);
  
  }

  fclose(f);

  printf("koniec.\n");

  return 0;
}

void init( int n )
{
  int i,j;
  
   for(i = 0; i < n; i++)
    for(j = 0; j < n; j++) {
      A[i][j] = i%2 + 1;
      B[i][j] = j%3 + 1;
      C[i][j] = (i+j)%4;
    }
}

/* Subtract the `struct timeval' values X and Y,
   storing the result in RESULT.
   Return 1 if the difference is negative, otherwise 0.  */

int
timeval_subtract (result, x, y)
  struct timeval *result, *x, *y;
{
  /* Perform the carry for the later subtraction by updating Y. */
  if (x->tv_usec < y->tv_usec) {
    int nsec = (y->tv_usec - x->tv_usec) / 1000000 + 1;
    y->tv_usec -= 1000000 * nsec;    
    y->tv_sec += nsec;
  }
  if (x->tv_usec - y->tv_usec > 1000000) {
    int nsec = (x->tv_usec - y->tv_usec) / 1000000;
    y->tv_usec += 1000000 * nsec;
    y->tv_sec -= nsec;
  }
  
  /* Compute the time remaining to wait.
     `tv_usec' is certainly positive. */
  result->tv_sec = x->tv_sec - y->tv_sec;
  result->tv_usec = x->tv_usec - y->tv_usec;
    
  /* Return 1 if result is negative. */
  return x->tv_sec < y->tv_sec;
}


/******************************************
Algorytm ikj
*******************************************/

#include <stdio.h>
#include <time.h>
#include <sys/times.h>
#include <sys/time.h>
#include <unistd.h>

#define MAXN 1024

float A[MAXN][MAXN],B[MAXN][MAXN],C[MAXN][MAXN];


void init(int n);
int timeval_subtract( struct timeval *result, 
		      struct timeval *x, 
		      struct timeval *y );



int main(void)
{
  register float A_tmp;

  int i,j,k;
  int n;

  FILE *f;

  struct timeval t_beg, t_end, t_diff;
  struct tms c_beg, c_end;
  clock_t c_diff;


  if((f = fopen("matmul_ikj.dat","w")) == NULL) {
    perror("Blad otwarcia pliku\n");
    exit(1);
  }

  fprintf(f,"floats: %d\n",sizeof(A[0][0]));  
  printf("zaczynamy: IKJ %d-%d\n",16,MAXN);

  for(n = 16; n <= MAXN; n *= 2) {
   
    printf("%4i init ",n);
    
    init(n);

    printf("before ");

    gettimeofday(&t_beg, NULL);

    printf("%li.%li:",t_beg.tv_sec,t_beg.tv_usec);
    fflush(stdout);

    /* ---------------------------------------- */

    times(&c_beg);
    gettimeofday(&t_beg, NULL);

    for(i = 0; i < n; i++)
      for(k = 0; k < n; k++) {
	A_tmp = A[i][k];
	for(j = 0; j < n; j++)
	  C[i][j] += A_tmp*B[k][j];
      }

    gettimeofday(&t_end, NULL);
    times(&c_end);
    
    /* ---------------------------------------- */

    timeval_subtract(&t_diff, &t_end, &t_beg);
    c_diff = c_end.tms_utime - c_beg.tms_utime;

    printf("= %li.%li\n", t_diff.tv_sec, t_diff.tv_usec);

    fprintf(f,"%i %li.%li %li %.20g\n",
	    n, t_diff.tv_sec, t_diff.tv_usec, c_diff,
	    ((double)c_diff)/CLOCKS_PER_SEC);
    fflush(f);
  
  }

  fclose(f);

  printf("koniec.\n");

  return 0;
}

void init( int n )
{
  int i,j;
  
   for(i = 0; i < n; i++)
    for(j = 0; j < n; j++) {
      A[i][j] = i%2 + 1;
      B[i][j] = j%3 + 1;
      C[i][j] = (i+j)%4;
    }
}

/* Subtract the `struct timeval' values X and Y,
   storing the result in RESULT.
   Return 1 if the difference is negative, otherwise 0.  */

int
timeval_subtract (result, x, y)
  struct timeval *result, *x, *y;
{
  /* Perform the carry for the later subtraction by updating Y. */
  if (x->tv_usec < y->tv_usec) {
    int nsec = (y->tv_usec - x->tv_usec) / 1000000 + 1;
    y->tv_usec -= 1000000 * nsec;    
    y->tv_sec += nsec;
  }
  if (x->tv_usec - y->tv_usec > 1000000) {
    int nsec = (x->tv_usec - y->tv_usec) / 1000000;
    y->tv_usec += 1000000 * nsec;
    y->tv_sec -= nsec;
  }
  
  /* Compute the time remaining to wait.
     `tv_usec' is certainly positive. */
  result->tv_sec = x->tv_sec - y->tv_sec;
  result->tv_usec = x->tv_usec - y->tv_usec;
    
  /* Return 1 if result is negative. */
  return x->tv_sec < y->tv_sec;
}


/******************************************
Algorytm ikj4 - wykorzystujacy technike loop unrolling
*******************************************/

#include <stdio.h>
#include <time.h>
#include <sys/times.h>
#include <sys/time.h>
#include <unistd.h>

#define MAXN 1024

float A[MAXN][MAXN],B[MAXN][MAXN],C[MAXN][MAXN];


void init(int n);
int timeval_subtract( struct timeval *result, 
		      struct timeval *x, 
		      struct timeval *y );


int main(void)
{
  register float A_tmp;

  int i,j,k, jj;
  int n, n4;

  FILE *f;

  struct timeval t_beg, t_end, t_diff;
  struct tms c_beg, c_end;
  clock_t c_diff;


  if((f = fopen("matmul_ikj4.dat","w")) == NULL) {
    perror("Blad otwarcia pliku\n");
    exit(1);
  }

  fprintf(f,"floats: %d\n",sizeof(A[0][0]));  
  printf("zaczynamy IKJ + loop unrolling(4): %d-%d\n",16,MAXN);

  for(n = 16; n <= MAXN; n *= 2) {
   
    n4 = n/4;
    printf("%4i init ",n);
    
    init(n);

    printf("before ");

    gettimeofday(&t_beg, NULL);

    printf("%li.%li:",t_beg.tv_sec,t_beg.tv_usec);
    fflush(stdout);

    /* ---------------------------------------- */

    times(&c_beg);
    gettimeofday(&t_beg, NULL);

    for(i = 0; i < n; i++)
      for(k = 0; k < n; k++) {
	A_tmp = A[i][k];
	for(j = 0; j < n; j += 4 ) {
	  C[i][j  ] += A_tmp*B[k][j  ];
	  C[i][j+1] += A_tmp*B[k][j+1];
	  C[i][j+2] += A_tmp*B[k][j+2];
	  C[i][j+3] += A_tmp*B[k][j+3];
	}
      }

    gettimeofday(&t_end, NULL);
    times(&c_end);
    
    /* ---------------------------------------- */

    timeval_subtract(&t_diff, &t_end, &t_beg);
    c_diff = c_end.tms_utime - c_beg.tms_utime;

    printf("= %li.%li\n", t_diff.tv_sec, t_diff.tv_usec);

    fprintf(f,"%i %li.%li %li %.20g\n",
	    n, t_diff.tv_sec, t_diff.tv_usec, c_diff,
	    ((double)c_diff)/CLOCKS_PER_SEC);
    fflush(f);
  
  }

  fclose(f);

  printf("koniec.\n");

  return 0;
}

void init( int n )
{
  int i,j;
  
   for(i = 0; i < n; i++)
    for(j = 0; j < n; j++) {
      A[i][j] = i%2 + 1;
      B[i][j] = j%3 + 1;
      C[i][j] = (i+j)%4;
    }
}

/* Subtract the `struct timeval' values X and Y,
   storing the result in RESULT.
   Return 1 if the difference is negative, otherwise 0.  */

int
timeval_subtract (result, x, y)
  struct timeval *result, *x, *y;
{
  /* Perform the carry for the later subtraction by updating Y. */
  if (x->tv_usec < y->tv_usec) {
    int nsec = (y->tv_usec - x->tv_usec) / 1000000 + 1;
    y->tv_usec -= 1000000 * nsec;    
    y->tv_sec += nsec;
  }
  if (x->tv_usec - y->tv_usec > 1000000) {
    int nsec = (x->tv_usec - y->tv_usec) / 1000000;
    y->tv_usec += 1000000 * nsec;
    y->tv_sec -= nsec;
  }
  
  /* Compute the time remaining to wait.
     `tv_usec' is certainly positive. */
  result->tv_sec = x->tv_sec - y->tv_sec;
  result->tv_usec = x->tv_usec - y->tv_usec;
    
  /* Return 1 if result is negative. */
  return x->tv_sec < y->tv_sec;
}